LAMMPS SPC/E水モデルの計算

LAMMPS
2023.11.04

古くからよく使われている剛体水分子モデルであるSPC/Eを使って水の密度を計算してみたいと思います。

スポンサーリンク

SPC/E水分子モデルについて

水が持つ様々な性質を精度よく再現するために、様々な種類の水モデルが考案されています。
その一つがSPC/E (extended simple point charge model, 拡張単純点電荷モデル)です。
元の論文は
Berendsen, Herman JC, J. Raul Grigera, and Tjerk P. Straatsma. “The missing term in effective pair potentials.” Journal of Physical Chemistry 91.24 (1987): 6269-6271.
他の3サイト剛体水分子モデルなども含め、パラメータ等はWikipediaに詳しく載っています。

水モデル - Wikipedia
ja.wikipedia.org
スポンサーリンク

計算内容

LAMMPS公式ドキュメントのサンプルコード

LAMMPSの公式ドキュメントページにSPC/Eを使った計算の例が載っています。

8.4.6. SPC water model — LAMMPS documentation
docs.lammps.org

今回はこちらのページの内容をベースにして計算していきます。

ファイルの準備

spce.molファイルを用意します。

# Water molecule. SPC/E geometry

3 atoms
2 bonds
1 angles

Coords

1    0.00000  -0.06461   0.00000
2    0.81649   0.51275   0.00000
3   -0.81649   0.51275   0.00000

Types

1        1   # O
2        2   # H
3        2   # H

Charges

1       -0.8476
2        0.4238
3        0.4238

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3

Shake Flags

1 1
2 1
3 1

Shake Atoms

1 1 2 3
2 1 2 3
3 1 2 3

Shake Bond Types

1 1 1 1
2 1 1 1
3 1 1 1

Special Bond Counts

1 2 0 0
2 1 1 0
3 1 1 0

Special Bonds

1 2 3
2 1 3
3 1 2

このspce.molは、計算入力ファイルから呼び出す形になります。
今回は水分子を500個用意しようと思います。
密度を1.0g/cm3とすると、水分子1個当たりの体積は29.89Å3です。
500倍した場合、立方体セルは一辺24.6Åくらいなので、少し余裕を持たせた適当な大きさ(-12.7 – 12.7)で作ります。
また、温度はコマンドライン引数で指定できるように変数${temperature}としています。
これは初速度を与えるvelocity allの行と、計算中の温度制御のfix npt行に書きました。
pair_styleはEwald法を使いたいため、lj/cut/coul/cut からlj/cut/coul/long に変更し、kspace_styleを指定しています。
Minimizeは省略、shakeのログも出力しない設定で、以下のようなファイルになりました。

units real
atom_style full
region box block -12.7 12.7 -12.7 12.7 -12.7 12.7
create_box 2 box  bond/types 1 angle/types 1 &
                extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2

mass 1 15.9994
mass 2 1.008

pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 1 1 0.1553 3.166
pair_coeff 1 2 0.0    1.0
pair_coeff 2 2 0.0    1.0

kspace_style ewald 1.0e-4

bond_style zero
bond_coeff 1 1.0

angle_style zero
angle_coeff 1 109.47

molecule water spce.mol
create_atoms 0 random 500 34564 NULL mol water 25367 overlap 1.33

timestep 1.0
fix rigid     all shake 0.0001 20 0 b 1 a 1
velocity all create ${temperature} 5463576
fix integrate all npt temp ${temperature} ${temperature} 100.0 iso 1.0 1.0 1000.0

thermo_style custom step temp press etotal density pe ke
thermo 1000
run 100000 upto
write_data spce.data nocoeff

ファイル名をin.spceとして、準備完了です。

追記
上の例は原子の定義が1：酸素、2：水素の場合ですが、1：水素、2：酸素の場合はpair_coeffのセクションが以下のようになります。

pair_style lj/cut/coul/long 10.0 10.0
pair_coeff 1 1 0.0    1.0
pair_coeff 1 2 0.0    1.0
pair_coeff 2 2 0.1553 3.166
スポンサーリンク

計算

以下のコマンドを実行します。
なお-varオプションにより、計算ファイル内の ${temperature} に 298.15 が代入されます。
このコマンドラインの数値を書き換えるだけで温度設定を変更できます。
バッチファイルやシェルスクリプトで温度を連番で指定して計算ができるので便利です。

lmp -in in.spce -var temperature 298.15

出力の一部を以下に貼っています。

molecule water spce.mol
Read molecule template water:
  1 molecules
  0 fragments
  3 atoms with max type 2
  2 bonds with max type 1
  1 angles with max type 1
  0 dihedrals with max type 0
  0 impropers with max type 0
create_atoms 0 random 500 34564 NULL mol water 25367 overlap 1.33
Created 1500 atoms
  using lattice units in orthogonal box = (-12.7 -12.7 -12.7) to (12.7 12.7 12.7)
  create_atoms CPU = 0.030 seconds

こちらは分子をランダム生成する部分です。
1500個の原子（500個の水分子）が生成されていることが確認できます。
セルが小さいとoverlapの設定により分子が削除されることがあり、その場合は分子数が減ります。

   Step          Temp          Press          TotEng        Density         PotEng         KinEng    
         0   298.15         127252.61      8295.0578      0.91277098     7407.218       887.83979    
      1000   278.50803     -1816.5944     -2536.7322      0.83690492    -3366.0815      829.34937    
      2000   297.62152      767.54138     -764.51772      0.85248793    -1650.7838      886.26605    
      3000   294.41703     -940.33075     -1482.4847      0.85726488    -2359.2084      876.72365    
      4000   303.02266     -1029.3307     -2280.3291      0.86752698    -3182.6788      902.34974    
      5000   309.55568     -1625.7826     -2763.2554      0.88416182    -3685.0594      921.80395  
       :
       :
    495000   299.14584      1735.5836     -1348.4354      0.98740201    -2239.2406      890.80522    
    496000   298.58532      3505.1897     -661.91297      0.9898038     -1551.0491      889.1361     
    497000   303.76877     -61.1881       -2650.2411      0.9864683     -3554.8126      904.57153    
    498000   303.2069      -508.79768     -2335.7307      0.98587143    -3238.6291      902.89837    
    499000   300.95014      1437.7024     -1446.994       0.98625999    -2343.1721      896.17812    
    500000   308.30472      1458.3828     -1400.4053      0.98735969    -2318.4841      918.07881    
Loop time of 770.237 on 4 procs for 500000 steps with 1500 atoms

Performance: 56.087 ns/day, 0.428 hours/ns, 649.151 timesteps/s, 973.726 katom-step/s
9.4% CPU use with 4 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 454.23     | 463.55     | 473.13     |  35.4 | 60.18
Bond    | 0.23932    | 0.24539    | 0.25098    |   0.8 |  0.03
Neigh   | 0.020919   | 0.021483   | 0.022327   |   0.4 |  0.00
Comm    | 55.791     | 55.9       | 56.019     |   1.1 |  7.26
Output  | 0.014818   | 0.015891   | 0.018985   |   1.4 |  0.00
Modify  | 232.24     | 241.71     | 251        |  48.8 | 31.38
Other   |            | 8.79       |            |       |  1.14

出力項目はthermo_style custom step temp press etotal density pe keにある通り、ステップ、温度、圧力、全エネルギー、密度、ポテンシャルエネルギー、運動エネルギーです。

密度をプロットしたところ以下のようになりました。初期密度からほぼ変わらず1近傍をキープしています。
計算後半部分の平均密度は0.9888 g/cm3でした。

データはPythonで処理しました。
使用したコードはこちらの記事で紹介しています。

LAMMPSのログをPython Matplotlibでグラフ化
LAMMPSの計算結果の温度やエネルギーなどが書かれたログをNumpyで読み込み、Matplotlibでグラフにします。非GUI環境でMatplotlibを使う方法も解説しています。
tech.gijukatsu.com
2023.11.04
スポンサーリンク

まとめ

LAMMPS公式で紹介されているSPC/E水モデルの計算ファイルに手を加えて、密度の計算をやってみました。
水専用に作られたパラメータなので密度が合うのは当然ですが、簡単に再現できました。

コメント

タイトルとURLをコピーしました