LAMMPSのインプットファイルの構造【バルクの水のMD計算】
更新サボっていて申し訳ありませんでした。久々の更新になります。今回は前回ダウンロードしたLAMMPSを使ってバルク水の分子動力学計算を行ってみます。
これはLAMMPSのインプットファイルの構造を学ぶだけなので、インプットファイルをどうやって作成するかというところは今のところは無視して構いません。
バルク水のMD計算のインプットファイル
units real atom_style full boundary p p p pair_style lj/cut/coul/long 10. 10. kspace_style pppm 1e-5 bond_style harmonic angle_style harmonic dihedral_style charmm read_data bulkwater.data neighbor 2.0 bin neigh_modify every 1 delay 0 check yes # 出力 dump 1 all custom 100 bulkwater.dump id type xs ys zs ix iy iz dump 2 all xtc 100 bulkwater.xtc # 速度の割り当て velocity all create 300.0 12345 run 0 velocity all scale 300.0 # 温度制御 fix 1 all shake 1e-4 1000 0 b 1 b 2 a 1 fix 2 all nvt temp 300.0 300.0 100. tchain 3 thermo_style multi thermo 100 # MDメイン計算 timestep 1.0 run 500000 write_restart bulkwater.restart
各コマンドの意味について以下で解説していきます。
units
シミュレーションに使用する単位の設定。例えばreal
なら質量はg/mol、距離はÅを使用する。他の単位については詳しくは以下のリンクを参照してください。
real
以外によく使用するコマンドはmetal
とlj
があるが、MD計算に使用するのは基本的にreal
とmetal
です。ソフトマターはreal
で金属はmetal
という考え方でいいと思います。
ここでは水なのでreal
を使用しますが、いずれはmetal
の場合も例示したいと思います。
http://lammps.sandia.gov/doc/units.html
atom_style
シミュレーションに使用する原子の型を決める。これも様々な設定方法がある(詳細は下記リンク参照)が、基本的に使用するのはunits real
の場合はfull
、units metal
の場合はatomic
で良いと思います。
大は小を兼ねるなので、あまり制限されたものを使用するとあとで大変なので大きい型(full
)を使うのがよいはず。
http://lammps.sandia.gov/doc/atom_style.html
boundary
境界条件の決定。boundary x y z
で各方向の境界条件を定める。p
は周期境界であり基本はこれを使用することが多いです。他にはf
やm
などがあり、前者は壁に到達すると原子を消去し、後者は原子が遠く離れてもそこを境界に再設定(Shrink-wrapping)する手法。
http://lammps.sandia.gov/doc/boundary.html
pair_style
原子間のポテンシャルを記述するのだが、ポテンシャルはMD計算の要なので、ここはかなり多数のコマンドが用意されているのでこれも詳しくは以下のマニュアルを参照してください。
http://lammps.sandia.gov/doc/pair_style.html
今回はlj/cut/coul/long
を使用します。この意味は短距離相互作用にはLennard-Jonesポテンシャルを、長距離相互作用のクーロン力にはEwald法、もしくはPPPM法を使用するという意味です。
数字はカットオフ値でlj/cut/coul/long 10. 10.
の場合、LJポテンシャルを10Åまで計算、クーロン力は10Åまで直接計算した後、それ以外のところでは上記手法などで計算することを意味しています。
原子のMDの場合はLJポテンシャルは10Å程度、クーロン力はEwald(PPPM)法を用いるのが最も多いです。
kspace_style
前項で述べた長距離相互作用の計算はk空間で行うので、ここはk空間の設定をします。pppm
はparticle-particle-particle-mesh Ewald法を意味しており、原子数が大きな系では計算時間がEwald法よりも速くなるという利点があります。Ewald
は通常のEwald法を使用します。その次の数字は許容される計算誤差で小さい方が精度が高く計算できます。
出先でパソコンの充電がやばいのでまたあとで追記します。
充電復活しましたので続きです。
bond_style
原子内のボンド間のポテンシャルを設定する。これも様々なコマンドが存在するので詳しくは以下のリンクを参照してください。
http://lammps.sandia.gov/doc/bond_style.html
使用頻度が高いのはharmonic
でE=K(r-r0)で表される調和振動子ポテンシャルです。
その他に比較的よく使われるのはMorseポテンシャルmorse
でE=D[1-exp{-α(r-r0)}]と表されます。
angle_style
原子内の結合角ポテンシャルの設定。こちらもbond_style
と同様にharmonic
を使用することが多い。表式はE=K(θ-θ0)となる。
他にも様々な記述方法があるのでリンク参照。
http://lammps.sandia.gov/doc/angle_style.html
dihedral_style
二面角ポテンシャルの設定。二面角の場合はcharmm
を使用することが多い。式はE=K{1-cos(nφ-d)}と表される。
harmonic
も存在するが、あまり使われない。(表式としてはcharmm
に似ている)
http://lammps.sandia.gov/doc/dihedral_style.html
read_data
Geometricalファイルの読み込み。Geomファイルの作り方は後日。
neighbor
粒子からneighbor listを構成するまでの距離と方法の指定を行う。
これをいじったらどう変化するのかよく分からん。色んな例を見ているとどれもneighbor 2.0 bin
になっているます。
http://lammps.sandia.gov/doc/neighbor.html
neigh_modify
neighbor listをどのタイミングで更新するかの設定。
例えばevery 1 delay 0 check yesなら0ステップ後に1ステップ毎に更新。(ただし実際はneighborで設定した距離の半分を移動したら更新する。)
http://lammps.sandia.gov/doc/neigh_modify.html
dump
dump
はアウトプットのフォーマットを指定します。Syntaxはdump ID group-ID style N file args
であり簡単な例はdump 2 all xtc 100 bulkwater.xtc
です。
all
は全原子に対してのアウトプットを意味しており、xtc
はアウトプットのファイルフォーマット、100
はアウトプット間隔、bulkwater.xtc
はアウトプットファイル名です。
ファイル名の後に、アウトプットする属性を指定することもできる。上記例ではid type xs ys zs ix iy iz
となっているが、id
は原子番号、type
は原子種番号、xs ys zs
は原子位置、ix iy iz
は周期境界条件のどこの仮想ボックスに位置しているかを出力します。
それ以外にも例えばvx vy vz
は原子の速度、fx fy fz
は原子に働く力、q
は原子電荷を出力します。それ以外のパラメータに関しては下記リンク参照。
http://lammps.sandia.gov/doc/dump.html
やったことはないが画像のアウトプットもできるらしい。構文はdump image
です。
http://lammps.sandia.gov/doc/dump_image.html
velocity
速度の割り当て。グループallの原子に対して300.0Kの温度を乱数12345で割り当てることを意味しています。
ここでshake法(後述)を使用する際は、速度を上手く割り当てられないことに注意してください。それを回避するには、
velocity all create 300.0 12345 run 0 velocity all scale 300.0
と、一度速度を割り当ててから計算を0ステップ実行して、その後300Kに再スケールを行うという方法があります。
http://lammps.sandia.gov/doc/velocity.html
fix
fix
コマンドはかなり多くのことができるので個別に記事を書きたいと思います。ここではSHAKE法の設定と、温度制御を行っています。
SHAKE法は水素原子を含む結合を固定する方法で、今回は水の計算にSPC/E力場を用いているので、shake法を使用します。fix 1 all shake 1e-4 1000 0 b 1 b 2 a 1
はall分子に対して1e-4の精度でSHAKE法を実行することを意味している。1000はSHAKE方の反復計算回数。0はSHAKE法を使用している結合と角度の情報の出力を何ステップごとに行うかということです。(0は出力しない。)
b
は結合長の制限、a
は角度の制限です。ここでは結合1と結合2を固定、角度1を固定しています。
http://lammps.sandia.gov/doc/fix_shake.html
また、fix 2 all nvt temp 300.0 300.0 100. tchain 3
は熱浴の設定です。今回はNVTアンサンブルを使用していますが、NPTアンサンブルnpt
やNPHアンサンブルnph
も使用できます。
temp
は温度制御の宣言。NPTの場合は圧力も制御するのでiso
で指定します。
300.0 300.0
は初期温度と最終温度で100.
は温度緩和の時間設定です。NVTアンサンブルの場合LAMMPSはNose-Hoover法を使用しているので温度緩和は100にすると良いとのこと。
tchain
は熱浴における振動制御のパラメータでNose-Hoover法の場合デフォルトは3。
http://lammps.sandia.gov/doc/fix_nh.html
thermo_style
出力する熱力学量の設定。multi
にするとetotal ke temp pe ebond eangle edihed eimp evdwl ecoul elong press
が出力される(全エネルギー、運動エネルギー、温度、全ポテンシャルエネルギー、結合エネルギー、結合角のエネルギー、二面角のエネルギー、improperな二面角のエネルギー、ファンデルワールスエネルギー、クーロンエネルギー、k空間のエネルギー、圧力)。
http://lammps.sandia.gov/doc/thermo_style.html
themo
アウトプット間隔の設定。
timestep
計算の時間刻みの設定。units
で指定した単位系で設定する。real
の場合はfs単位。水素原子が存在するので1fs刻みにする。
run
計算時間の設定。real
の場合はfs単位。
write_restart
リスタートファイルの作成。