LAMMPSでCO2を計算する EPM, EPM2, TraPPE

LAMMPSで剛体のCO₂分子を計算します。

パラメータ
計算入力ファイル
計算実行＆結果

パラメータ

EPM、EPM2、TraPPEの3種類の剛体モデルで計算をします。表の数値はLAMMPSのunits realでそのまま使用できます。

	$$ \epsilon_{C-C} $$[kcal/mol]	$$ \sigma_{C-C} $$[Å]	$$ \epsilon_{O-O} $$[kcal/mol]	$$ \sigma_{O-O} $$[Å]	$$ q_{C} $$[e]	$$ q_{O} $$[e]
EPM	0.057627	2.785	0.164932	3.064	+0.6645	-0.33225
EPM2	0.055898	2.757	0.159984	3.033	+0.6512	-0.3256
TraPPE	0.053655	2.800	0.156989	3.050	+0.700	-0.350

EPM, EPM2 https://doi.org/10.1021/j100031a034
TraPPE https://doi.org/10.1002/aic.690470719

それぞれLAMMPSのpair_coeffの記載は以下のようになります。

#EPM
pair_style lj/cut/coul/long 10.0 8.0
pair_modify mix arithmetic
pair_coeff 1 1 0.057627 2.785
pair_coeff 2 2 0.164932 3.064

#EPM2
pair_style lj/cut/coul/long 10.0 8.0
pair_modify mix arithmetic
pair_coeff 1 1 0.055898 2.757
pair_coeff 2 2 0.159984 3.033

#TraPPE
pair_style lj/cut/coul/long 10.0 8.0
pair_modify mix arithmetic
pair_coeff 1 1 0.053655 2.800 
pair_coeff 2 2 0.156989 3.050

また、分子をmoleculeコマンドでテキストファイルを読み込みます。テキストファイルには、分子中の原子座標、電荷、BondsとAnglesの定義が書かれています。

EPM.txt

# CO2 molecule. EPM geometry
# header section:
3 atoms
2 bonds
1 angles

# body section:
Coords

1    0.00000   0.00000   0.00000
2    1.16100   0.00000   0.00000
3   -1.16100   0.00000   0.00000

Types

1        1   # C
2        2   # O
3        2   # O

Charges

1        0.6645
2       -0.33225
3       -0.33225

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3

EPM2.txt

# CO2 molecule. EPM2 geometry
# header section:
3 atoms
2 bonds
1 angles

# body section:
Coords

1    0.00000   0.00000   0.00000
2    1.14900   0.00000   0.00000
3   -1.14900   0.00000   0.00000

Types

1        1   # C
2        2   # O
3        2   # O

Charges

1        0.6512
2       -0.3256
3       -0.3256

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3

TraPPE.txt

# CO2 molecule. TraPPE geometry
# header section:
3 atoms
2 bonds
1 angles

# body section:
Coords

1    0.00000   0.00000   0.00000
2    1.16000   0.00000   0.00000
3   -1.16000   0.00000   0.00000

Types

1        1   # C
2        2   # O
3        2   # O

Charges

1        0.70
2       -0.35
3       -0.35

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3

計算入力ファイル

LAMMPSのinファイルは以下のような形式で書きました。（EPM2の場合）

EPM2_0.in

units real
atom_style full
boundary p p p

variable length equal 100.0
region box block 0.0 ${length} 0.0 ${length} 0.0 ${length}
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 12.011
mass 2 15.9994

pair_style lj/cut/coul/long 10.0 8.0
pair_modify mix arithmetic
pair_coeff 1 1 0.055898 2.757
pair_coeff 2 2 0.159984 3.033

bond_style zero
bond_coeff 1 1.149

angle_style zero
angle_coeff 1 180.0

kspace_style ewald 1.0e-3

molecule co2 EPM2.txt
create_atoms 0 random 500 34564 NULL mol co2 25367 overlap 5.0

timestep 1.0


group co2_mols type 1 2
velocity all create ${temp} 5463576
fix fixgiridnpt co2_mols rigid/npt molecule temp ${temp} ${temp} 100.0 iso ${press} ${press} 1000.0

thermo_style custom step temp press etotal density pe ke
thermo 1000
run 100000 upto

write_data EPM2_1.data nocoeff

大きめのボックスを作り、EPM2のパラメータを定義し、分子を500個生成してNPTアンサンブルの計算を行います。温度と圧力は変数にしてあるので、実行時に以下のように指定します。

# 45℃ 20MPaの場合
lmp -var temp 318.15 -var press 197.385 -in EPM2_0.in

圧力の単位はatmです。197.385 atm ＝ 20 MPa

このEPM2_0.inの計算の最後にwrite_dataでEPM2_1.dataを作成します。これは密度緩和が完了した構造データです。このデータの続きから、次のinファイルを使って再計算をします。

EPM2_1.in

units real
atom_style full
boundary p p p

read_data EPM2_1.data

pair_style lj/cut/coul/long 10.0 8.0
pair_modify mix arithmetic
pair_coeff 1 1 0.055898 2.757
pair_coeff 2 2 0.159984 3.033

bond_style zero
bond_coeff 1 1.149

angle_style zero
angle_coeff 1 180.0

kspace_style ewald 1.0e-3

timestep 1.0

group co2_mols type 1 2
fix fixgiridnpt co2_mols rigid/npt molecule temp ${temp} ${temp} 100.0 iso ${press} ${press} 1000.0

thermo_style custom step temp press etotal density pe ke
thermo 1000
run 100000 upto

write_data EPM2_2.data nocoeff

計算実行＆結果

EPM2_0.inとEPM2_1.inの2段階の計算を実行します。

lmp -var temp 318.15 -var press 197.385 -in EPM2_0.in |tee 0.log
lmp -var temp 318.15 -var press 197.385 -in EPM2_1.in |tee 1.log

EPM2だけでなく、EPMとTraPPEでも同様にファイルを用意して計算を実行します。

そして第2段階の計算ログから密度の平均値を求めます。

Model	Density [g/cm³]
EPM	0.6857
EPM2	0.6541
TraPPE	0.6934

45℃ 20MPaの超臨界CO₂の密度は0.8 g/cm³くらいですので、今回の計算結果はいずれも少し低いです。大きく外れたわけではないですが、計算条件を見直すことで改善はしそうです。とりあえずちゃんと動いたからヨシ！