アルゴンの分子動力学計算
やったことなかったので、アルゴンの分子動力学計算をやってみた。まぁ、ただのGromacsの使い方メモのような。使用したGromacsのバージョンは、4.6.1
原理的なことに関しては、
分子動力学シミュレーションにおける温度・圧力制御第1回=能勢の熱浴と能勢・フーバー熱浴
https://www.jstage.jst.go.jp/article/mssj1998/10/44/10_44_29/_article/-char/ja/
分子動力学シミュレーションにおける温度・圧力制御 第3回：速度スケーリング法，ガウス束縛法，ベレンゼン熱浴
https://www.jstage.jst.go.jp/article/mssj/11/2/11_2_2_43/_article/-char/ja/
あたりを見ると、大体書いてある。まぁ、MDは原理自体は難しいことは、あまりないと思う。
ポテンシャルは、Lennard-Jonesポテンシャルで、あまり関係ないけど
量子力学はLennard-Jonesポテンシャルを再現できるか
http://d.hatena.ne.jp/m-a-o/20121217#p2
も参照。きっと、2030年くらいになれば、怪しい近似ポテンシャルとか力場とか使わなくても、量子力学的に原子核間に働くヘルマン・ファインマン力を計算して、原子核を古典論的に時間発展させるくらいのことは、そこらのPCで出来るようになるかもしれないけど、現在のところ厳しそう
(1)トポロジーファイルというものを作る。今の場合、実質的には、Lennard-Jonesパラメータと原子数の指定くらいしかやってない
;; argon.top [ defaults ] ; nbfuncは非結合粒子間ポテンシャルで、LJポテンシャルかBuckinghamポテンシャルを指定する(LJなら、1) ; comb-ruleは、異原子間のLJパラメータとかの決め方(sigma/epsilonを指定する場合は、2か3)。今はどっちでも同じ ; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ 1 3 no 1.0 1.0 [ atomtypes ] ; name mass charge ptype sigma(nm) epsilon(kJ/mol) AR 39.948 0.0 A 0.34 0.99774 [ moleculetype ] ; name nrexcl Ar 1 [ atoms ] ; id type resnr residu atom cgnr charge 1 AR 1 Ar Ar 1 0 [ system ] Argon [ molecules ] AR 512
(2)MD計算の初期状態を指定する。
pdbファイルで指定できるけど、その場合は、初速度の指定はできない(後で、Gromacsが適当な速度を付加してくれる)
速度も指定したい場合は、.groファイルフォーマットというのを使えばいいっぽい
とりあえず、何も考えず、立方格子の頂点に原子を配置するpdbを生成するコードをpythonで作る
#-*- coding:utf-8 -*-
#CellSizeの単位は(Å)
def gen_cubic_pdb(CellSize , N):
L = CellSize
A = CellSize/N
print('HEADER argon')
print('CRYST1 %9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d' % (L,L,L,90.0,90.0,90.0,"P 1",1))
atom_format = "ATOM %(serial)5d %(name)-4s %(resName)-4s%(chainID)1s%(resSeq)4d%(iCode)1s %(x)8.3f%(y)8.3f%(z)8.3f%(occupancy)6.2f%(tempFactor)6.2f"
for nx in xrange(N):
for ny in xrange(N):
for nz in xrange(N):
vars = dict()
vars['serial'] = nx*N*N+ny*N+nz+1
vars['x'] = nx * A
vars['y'] = ny * A
vars['z'] = nz * A
vars['name'] = 'Ar'
vars['resName'] = 'Ar'
vars['resSeq'] = nx*N*N+ny*N+nz+1
vars['iCode'] = ''
vars['occupancy'] = 1.0
vars['chainID'] = 'A'
vars['tempFactor'] = 0.0
print(atom_format % vars)
if __name__=="__main__":
gen_cubic_pdb(250.0 , 8)python genpdb.py > init.pdb
とかで、初期状態を指定する。トポロジーファイルに原子数を指定してあるので、一致させる
(3)MDの計算条件の指定(タイムステップやステップ数や相互作用のカットオフ等)
指定できる項目は、凄く多い
http://www.gromacs.org/Documentation/File_Formats/.mdp_File
とりあえず、最小に近い記述は、以下の通り(温度を273Kに保つという以外は、技術的な詳細)
;; NVT.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.005 ;; timestep(ps) 5fs nsteps = 200000 ;; total time=1000ps=1ns ; Output frequency and precision for xtc file nstxtcout = 1000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 273.0 ; Temperature coupling: no(NVE) , nose-hoover , berendsen tcoupl = berendsen tc-grps = system tau-t = 0.2 ;; time constant(ps) ref-t = 273.0 ;;reference temperature(K)
(4)MD計算
以上の入力ファイルを、まとめて、gromppというのにかけると、単一の入力ファイルを作ってくれる
grompp -p argon.top -c init.pdb -f NVT.mdp -o md.tpr
md.tprができれば成功。で、あとはmdrun
mdrun -deffnm md
出力は
md.xtc : trajectory
md.edr : エネルギー
md.gro : 最終的な位置と速度(継続して、別の計算に投げる時とかに使う)
md.log : ログファイル
など
(5)計算結果の確認と解析
g_energy -f md.edr -o energy -xvg none
とかやって、
End your selection with an empty line or a zero. ------------------------------------------------------------------- 1 LJ-(SR) 2 Disper.-corr. 3 Coulomb-(SR) 4 Potential 5 Kinetic-En. 6 Total-Energy 7 Temperature 8 Pres.-DC 9 Pressure 10 Vir-XX 11 Vir-XY 12 Vir-XZ 13 Vir-YX 14 Vir-YY 15 Vir-YZ 16 Vir-ZX 17 Vir-ZY 18 Vir-ZZ 19 Pres-XX 20 Pres-XY 21 Pres-XZ 22 Pres-YX 23 Pres-YY 24 Pres-YZ 25 Pres-ZX 26 Pres-ZY 27 Pres-ZZ 28 #Surf*SurfTen 29 T-System 30 Lamb-System
などと出てきたら、適当に
4 5 6 7 9
と打つと、energy.xvgというファイルが作られて、一定ステップごとに、ポテンシャルエネルギー・運動エネルギー・全エネルギー・温度・圧力の時間変化が見れる。とりあえず、温度が大体273Kに保たれていることとかが見て取れる。あと
Energy Average Err.Est. RMSD Tot-Drift ------------------------------------------------------------------------------- Potential -4.0542 0.18 1.72151 -0.99296 (kJ/mol) Kinetic En. 1739.87 0.031 1.07103 -0.172205 (kJ/mol) Total Energy 1735.81 0.2 1.36278 -1.1652 (kJ/mol) Temperature 273.003 0.0049 0.168055 -0.027022 (K) Pressure 1.23123 0.00017 0.00798183 -0.000493924 (bar)
とかいう平均値が出る。NVTの値から、理想気体の状態方程式によって、圧力を見積もると
p = NkT/V = 512*(1.3806488*1e-23)*273/(250*1e-10)^3 = 123508.202...(Pa)=1.23508...(bar)
くらいなので、大体合ってそう。Err.Est.が凄く小さいのは、多分最初からほぼ平衡状態にあることを示していると思う
MD計算で分かる物理量は、拡散係数、粘性係数(ずり粘性率、体積粘性率)、熱伝導率、比熱、密度、沸点、融点とか色々あると思うけど、一番簡単なのは、自己拡散係数で、平均二乗変位から見積もる方法と、速度相関関数からGreen-Kubo公式で求める方法があって、前者は、Gromacsに付いてるg_msdで簡単に計算できて
g_msd -s md.tpr -f md.xtc -o msd -xvg none
とかやって、適当に"0"と打つと、D=1.086e+04 (+/- 2985) 1e-5 cm^2/sであるとか返ってくる。
アルゴンガスの自己拡散係数のデータは探しても、あまり出てこないけど
DIFFUSION COEFFICIENT@thermopedia
http://www.thermopedia.com/content/696/
によると、常温常圧で、0.157(cm^2/s)とある。今、MD計算の方は、圧力を常圧にしてない(NPTアンサンブルでの計算は(6)を参照)けど、オーダーは合ってそう
g_msd -s md.tpr -f md.xtc -o msd -xvg none -beginfit 100 -endfit 600
とかやると、100psから600psのデータを使って自己拡散係数を出してくれるッぽい(今の場合、デフォルトでは100~900ps)。計算値は
100~400ps : 8346.9463 (+/- 2933.8926) 1e-5 cm^2/s 100~600ps : 9714.9385 (+/- 3332.4248) 1e-5 cm^2/s 100~800ps : 1.057e+04 (+/- 3164) 1e-5 cm^2/s 100~900ps : 1.086e+04 (+/- 2985) 1e-5 cm^2/s
1ns程度だと、まだ若干収束してない感じがする。計算時間を2nsと4nsに増やすと、それぞれ、1.242e+04 (+/- 1964) 1e-5 cm^2/sと1.287e+04 (+/- 166.5) 1e-5 cm^2/sくらいになって、標準偏差も大分小さくなるので、0.13(cm^2/s)弱くらいで収束するっぽい
Green-Kubo公式による自己拡散係数の計算もGromacsに付いているg_velaccとg_analyzeというので、できるらしい。
Diffusion Constant
http://www.gromacs.org/Documentation/How-tos/Diffusion_Constant
g_velaccというので、自己速度相関関数を計算してくれるっぽい。けど、そのためには中間速度を出力させる必要があって、それには、NVT.mdpに次の3行を追加して計算する。こうすると、md.trrとかいうファイルができるはず
nstxout = 100 nstvout = 100 nstfout = 100
計算できたら、
g_velacc -nonormalize -mol -f md.trr -s md.tpr g_analyze -f vac.xvg -integrate
とかで
Calculating the integral using the trapezium rule Integral 1 34.11745 +/- 0.00000
というのが出てくる。単位は、(nm^2/ps)らしいので、(cm^2/s)に直すためには、10^(-2)をかけて、これは速度自己相関関数の積分値なので、拡散係数にするために、3で割ると、0.1137(cm^2/s)という値になる。
(6)沸点・融点の計算
今度は、NPT一定の条件で、ちょっとずつ温度を下げていって、気液相転移・固液相転移を見る。歴史的には、Alderが、2次元系で剛体球ポテンシャルで、固液相転移を確認したのが、分子動力学の始まりらしい。
topファイルとpdbファイルは、さっきと同じで、mdpを以下のように変える
;; anneal.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.001 ;; timestep(ps) 1fs nsteps = 40000000 ;; total time=40ns ; Output frequency and precision for xtc file nstenergy = 10000 nstlog = 10000 nstxtcout = 10000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 150.0 ;; temperature and pressure coupling tcoupl = berendsen tc-grps = system tau_t = 0.2 ;; time constant(ps) ref_t = 150.0 ;;reference temperature(K) pcoupl = berendsen pcoupltype = isotropic nstpcouple = 1 tau_p = 0.015 ;; time constant(ps) compressibility = 4.5e-5 ;; (1/bar) ref_p = 1.0 ;;reference pressure(bar) ; SIMULATED ANNEALING annealing = single ; Type of annealing for each temperature group (no/single/periodic) annealing_npoints = 11 ; Number of time points annealing_time = 0 2800 3000 4800 5000 14800 15000 24800 25000 34800 35000 ; List of times(ps) at the annealing points annealing_temp = 150 150 110 110 90 90 70 70 50 50 30 ; Temp. at each annealing point
このように書くと、0~2800psでは、熱浴の温度を150Kに保ち、2800~3000psの間に、150Kから110Kへ冷却し(熱浴の温度をリニアに下げていく)、3000~4800psでは、温度を110Kに保ち...という感じの計算をやってくれる。圧力も1barに保つように設定されている。で、以下のコマンドで計算
grompp -p argon.top -c init.pdb -f anneal.mdp -o md.tpr mdrun -deffnm md
計算が終わると、今度はg_energyコマンドで、密度や体積も出力できるようになっている。各温度で平衡状態に於ける密度を見ると
150K:3.1(kg/m^3) 70K:7.4(kg/m^3) 50K:1690(kg/m^3) 30K:1742(kg/m^3)
くらいになっていて、70~50Kの間で、密度の急激な増大が起きていることが確認できる。実験値は
気体アルゴン密度(0℃、1気圧):1.784(kg/m^3) 液体アルゴン密度(87K、液体):1380(kg/m^3) 固体アルゴン密度(40K、固体):1656(kg/m^3) 融点:84K 沸点:87K
なので、融点・沸点は全然ずれている。密度は、大体合っているっぽい。一応、最終配置から、30Kに保って、自己拡散係数を出してみると、液体でなく、固体になっているっぽい
逆に、固体から始めて温度をあげていく。以下のようなpythonコードで、fcc結晶のpdbを作成する。今度は、原子数は3x3x3x4=108個
#!/usr/bin/env python
def gen_fcc(a , Nx , Ny , Nz):
ret = set([])
for n1 in xrange(2*Nx):
for n2 in xrange(2*Ny):
for n3 in xrange(2*Nz):
if (n1+n2+n3)%2==1:continue
x = n1*0.5*a
y = n2*0.5*a
z = n3*0.5*a
yield (x,y,z)
if __name__=="__main__":
a = 5.26 #-- lattice constant
Nx,Ny,Nz = 3,3,3
print('HEADER argon')
print('CRYST1 %9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %-11s%4d' % (Nx*a,Ny*a,Nz*a,90.0,90.0,90.0,"P 1",1))
atom_format = "ATOM %(serial)5d %(name)-4s %(resName)-4s%(chainID)1s%(resSeq)4d%(iCode)1s %(x)8.3f%(y)8.3f%(z)8.3f%(occupancy)6.2f%(tempFactor)6.2f"
for n,(x,y,z) in enumerate(gen_fcc(a,Nx,Ny,Nz)):
vars = dict()
vars['serial'] = n+1
vars['x'] = x
vars['y'] = y
vars['z'] = z
vars['name'] = 'Ar'
vars['resName'] = 'Ar'
vars['resSeq'] = n+1
vars['iCode'] = ''
vars['occupancy'] = 1.0
vars['chainID'] = 'A'
vars['tempFactor'] = 0.0
print(atom_format % vars).mdpは以下のように適当に作った(.topファイルの原子数も108に変更しておくこと)
;; heatup.mdp ; RUN CONTROL PARAMETERS integrator = md dt = 0.001 ;; timestep(ps) 1fs nsteps = 40000000 ;; total time=40ns ; Output frequency and precision for xtc file nstenergy = 100000 nstlog = 100000 nstxtcout = 100000 xtc-precision = 1000 ; Neighbour search algorithm (simple or grid) ns-type = grid nstlist = 1 ;; nblist update frequency rlist = 0.5 ;;(nm) rcoulomb = 0.5 rvdw = 0.5 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = EnerPres ; generate velocities gen-vel = yes gen-temp = 50.0 ;; temperature and pressure coupling tcoupl = berendsen tc-grps = system tau_t = 0.2 ;; time constant(ps) ref_t = 50.0 ;;reference temperature(K) pcoupl = berendsen pcoupltype = isotropic nstpcouple = 1 tau_p = 0.015 ;; time constant(ps) compressibility = 4.5e-5 ;; (1/bar) ref_p = 1.0 ;;reference pressure(bar) ; SIMULATED ANNEALING annealing = single ; Type of annealing for each temperature group (no/single/periodic) annealing_npoints = 11 ; Number of time points annealing_time = 0 2800 3000 4800 5000 14800 15000 24800 25000 34800 35000 ; List of times(ps) at the annealing points annealing_temp = 50 50 60 60 80 80 100 100 120 120 160 ; Temp. at each annealing point
今度は、平衡状態に達した後で、各温度での密度が
50K:1750(kg/m^3) 80K:1680(kg/m^3) 100K:1620(kg/m^3) 120K:1290(kg/m^3) 160K:3.01(kg/m^3)
くらいになっているので、100~120Kで液体化して、120~160Kの間で気体化しているっぽい。厳密には、もっと細かく点を取るべきだと思うけど、冷却していった時と、加熱していった時で、沸点と融点が全然合ってそうにないのは、何故だろう？一応、これらの値は、原子数を500に増やしても、ほぼ同一だった
(7)疑問
心当たりのある方は、教えてください
・沸点や融点が実験値とかなりずれる理由は何？(単純に、ポテンシャルのfittingがよくない?)
・固体から加熱していった時と、気体から冷却していった時で、融点・沸点が大きくずれる理由は何？経験的に知られているように、固体の液化温度と液体の固化温度が概ね一致するための条件は？
