トップ > 分子動力学計算（LAMMPS） > LAMMPSのインプットファイルの構造【バルクの水のMD計算】

LAMMPSのインプットファイルの構造【バルクの水のMD計算】

分子動力学計算（LAMMPS） LAMMPS

更新サボっていて申し訳ありませんでした。久々の更新になります。今回は前回ダウンロードした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

前項で述べた長距離相互作用の計算はｋ空間で行うので、ここはｋ空間の設定をします。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

リスタートファイルの作成。