きっかけ
以前、液体関係の古典分子動力学をやっていた。しかしずっと機械学習ポテンシャルを使ってみたかったので、手始めに機械学習ポテンシャルの餌である第一原理計算データの作成やってみた。
参考図書:動かして理解する第一原理電子状態計算
ハード・ソフト情報
ハード情報
- CPU: Ryzen7 3800X 8コア16スレッド
- メモリ: DDR4 32GB
ソフト情報
- OS: Ubuntu-22.04
- 第一原理計算: quantumEspresso(v.6.7MAx)
- 分子集団作成: packmol(20.15.0)
- データ解析: python(3.10.12),MDAnalysis(2.8.0), matplotlib(3.8.4)
- データ可視化: OVITO
液体分子集団の作成
今回は20個の水分子集団を作成する。
適当に水分子一個の情報を書き込んだxyzファイル water.xyzを作成する。
3
Water molecule
O 0.000 0.000 0.000
H 0.757 0.586 0.000
H -0.757 0.586 0.000
次に、packmol用のインプットファイル water.inpを作成
tolerance 2.0
filetype xyz
output water_20.xyz
structure water.xyz
number 20
inside box 0. 0. 0. 10. 10. 10.
end structure
packmolを実行
packmol < water.inp
第一原理計算の実行
擬ポテンシャルの用意
計算に使用した2種類の原子それぞれの擬ポテンシャルは以下のものを選んだ。
- 水素: H.pbe-rrkjus_psl.1.0.0.UPF
- 酸素: O.pbe-n-rrkjus_psl.1.0.0.UPF
quantumEspressoの実行
本当はcar-parinello法を使いたかったが、うまいこと行かなかったのでBorn-Oppenheimer法を使った。
以下、quantumEspresso用スクリプトをwater.inとして保存する。
&CONTROL
calculation='md'
dt=10
nstep=1000
pseudo_dir='pseudo_dir/'
prefix = 'prefix',
tprnfor = .true.,
tstress = .true.,
verbosity = 'high',
/
&SYSTEM
ntyp = 2
nat = 60
ibrav = 0
nosym=.true.
/
&ELECTRONS
/
&IONS
ion_temperature = 'rescaling',
tempw = 300.0,
/
ATOMIC_SPECIES
H 1.0079 H.pbe-rrkjus_psl.1.0.0.UPF
O 15.9994 O.pbe-n-rrkjus_psl.1.0.0.UPF
K_POINTS gamma
CELL_PARAMETERS angstrom
10.0000000000000 0.00000000000000 0.00000000000000
0.00000000000000 10.0000000000000 0.00000000000000
0.00000000000000 0.00000000000000 10.0000000000000
ATOMIC_POSITIONS angstrom
# packmolで作成した原子の種類と位置の情報をここにコピペする
Born-Oppenheimer法なので、pw.xというソルバーを使う
pw.x < water.in > water.out
計算には8コアのCPUで16時間を要した。
計算結果の解析
可視化
出力されたoutファイルからどうにかOVITOで読めるようにxyzファイルに変換するpyhtonスクリプトを作成した。
def extract_xyz_from_cpout(cpout_path, xyz_output_path):
trajectory = []
with open(cpout_path, "r") as f:
lines = f.readlines()
current_atoms = []
recording = False
for line in lines:
if "ATOMIC_POSITIONS" in line:
recording = True
if current_atoms:
trajectory.append(current_atoms)
current_atoms = []
continue
elif recording:
if line.strip() == "" or not line[0].isalpha():
recording = False
if current_atoms:
trajectory.append(current_atoms)
current_atoms = []
else:
current_atoms.append(line.strip())
if current_atoms:
trajectory.append(current_atoms)
if not trajectory:
raise ValueError("No atomic positions found in the file.")
num_atoms = len(trajectory[0])
with open(xyz_output_path, "w") as f:
for i, frame in enumerate(trajectory):
f.write(f"{num_atoms}\n")
f.write(f"Step {i+1}\n")
for atom_line in frame:
f.write(f"{atom_line}\n")
if __name__ == "__main__":
import sys
if len(sys.argv) != 3:
print("Usage: python extract_cpout_xyz.py name.out name.xyz")
else:
extract_xyz_from_cpout(sys.argv[1], sys.argv[2])
水素結合は、2.3Åの位置に来たOとHを結んでいる。かわいい動きをしている。
動径分布関数
動径分布関数はMDAnalysisの機能を使って表示する。
import MDAnalysis as mda
from MDAnalysis.analysis.rdf import InterRDF
import matplotlib.pyplot as plt
# ファイル名を指定
xyz_file = "name.xyz" # あなたがアップロードしたファイル名と一致させること
box = [10.0, 10.0, 10.0, 90, 90, 90] # 推定セルサイズ（必要に応じて調整）
# トラジェクトリ読み込み
u = mda.Universe(xyz_file, format="XYZ")
# 各フレームにボックスサイズを設定
for ts in u.trajectory:
ts.dimensions = box
# O-H間のRDFを計算
g_OH = InterRDF(u.select_atoms("name O"), u.select_atoms("name H"),nbins=100, range=(1, 10.0))
g_OH.run()
# プロット
plt.plot(g_OH.bins, g_OH.rdf)
plt.xlabel("r (Å)")
plt.ylabel("g_OH(r)")
plt.title("Radial Distribution Function (O-H)")
plt.grid()
plt.tight_layout()
plt.savefig("rdf_OH.png")
plt.show()
-
酸素ー水素の動径分布 $g_{OH}$
-
水素ー水素の動径分布 $g_{HH}$
-
酸素ー酸素の動径分布 $g_{OO}$
下に挙げた先行実験の結果から見ても、いい感じではなかろうか
- https://fukuoka-u.repo.nii.ac.jp/record/2194/files/S3801_0095.pdf
- https://rdreview.jaea.go.jp/review_jp/2022/j2022_9_3f9_6.html
DeepMD-kitにいれるデータ数について
いくつか文献を見ると、数百個のデータセットを使って学習するらしい。
単純計算で16×100時間=1600時間=66日
途方もない