有機分子の分子動力学計算
有機分子の描画とlammpsファイルの取得
今回使用する有機分子は1pentanolを使用する
まずはhttps://www.nite.go.jp/chem/chrip/chrip_search/srhInput
で1-pentanolを検索し、以下のSMILESを取得する。
OCCCCC
このSMILESをhttp://zarbi.chem.yale.edu/ligpargen/
になげて投げてOPLS力場のパラメータを得る。
すると以下のようなファイルが得られる。
LAMMPS data file Created by LigParGen - (Written by Leela S. Dodda)
18 atoms
17 bonds
31 angles
39 dihedrals
10 impropers
18 atom types
17 bond types
31 angle types
39 dihedral types
10 improper types
-4.863660 45.136340 xlo xhi
0.107370 50.107370 ylo yhi
-0.930060 49.069940 zlo zhi
Masses
1 15.999
2 12.011
3 12.011
4 12.011
5 12.011
6 12.011
7 1.008
8 1.008
9 1.008
10 1.008
11 1.008
12 1.008
13 1.008
14 1.008
15 1.008
16 1.008
17 1.008
18 1.008
Pair Coeffs
1 0.170 3.1200000
2 0.066 3.5000000
3 0.066 3.5000000
4 0.066 3.5000000
5 0.066 3.5000000
6 0.066 3.5000000
7 0.000 0.0000000
8 0.030 2.5000000
9 0.030 2.5000000
10 0.030 2.5000000
11 0.030 2.5000000
12 0.030 2.5000000
13 0.030 2.5000000
14 0.030 2.5000000
15 0.030 2.5000000
16 0.030 2.5000000
17 0.030 2.5000000
18 0.030 2.5000000
Bond Coeffs
1 320.0000 1.4100
2 268.0000 1.5290
3 268.0000 1.5290
4 268.0000 1.5290
5 268.0000 1.5290
6 553.0000 0.9450
7 340.0000 1.0900
8 340.0000 1.0900
9 340.0000 1.0900
10 340.0000 1.0900
11 340.0000 1.0900
12 340.0000 1.0900
13 340.0000 1.0900
14 340.0000 1.0900
15 340.0000 1.0900
16 340.0000 1.0900
17 340.0000 1.0900
Angle Coeffs
1 50.000 109.500
2 58.350 112.700
3 58.350 112.700
4 58.350 112.700
5 55.000 108.500
6 35.000 109.500
7 35.000 109.500
8 37.500 110.700
9 37.500 110.700
10 37.500 110.700
11 37.500 110.700
12 37.500 110.700
13 37.500 110.700
14 37.500 110.700
15 37.500 110.700
16 37.500 110.700
17 37.500 110.700
18 37.500 110.700
19 33.000 107.800
20 37.500 110.700
21 33.000 107.800
22 33.000 107.800
23 37.500 110.700
24 37.500 110.700
25 37.500 110.700
26 37.500 110.700
27 33.000 107.800
28 33.000 107.800
29 33.000 107.800
30 37.500 110.700
31 33.000 107.800
Dihedral Coeffs
1 -1.552 0.000 0.000 0.000
2 1.300 -0.200 0.200 0.000
3 1.300 -0.200 0.200 0.000
4 -0.356 -0.174 0.492 0.000
5 0.000 0.000 0.300 0.000
6 0.000 0.000 0.300 0.000
7 0.000 0.000 0.300 0.000
8 0.000 0.000 0.300 0.000
9 0.000 0.000 0.468 0.000
10 0.000 0.000 0.300 0.000
11 0.000 0.000 0.300 0.000
12 0.000 0.000 0.300 0.000
13 0.000 0.000 0.300 0.000
14 0.000 0.000 0.300 0.000
15 0.000 0.000 0.300 0.000
16 0.000 0.000 0.300 0.000
17 0.000 0.000 0.300 0.000
18 0.000 0.000 0.300 0.000
19 0.000 0.000 0.300 0.000
20 0.000 0.000 0.300 0.000
21 0.000 0.000 0.300 0.000
22 0.000 0.000 0.300 0.000
23 0.000 0.000 0.300 0.000
24 0.000 0.000 0.300 0.000
25 0.000 0.000 0.300 0.000
26 0.000 0.000 0.300 0.000
27 0.000 0.000 0.300 0.000
28 0.000 0.000 0.300 0.000
29 0.000 0.000 0.352 0.000
30 0.000 0.000 0.300 0.000
31 0.000 0.000 0.300 0.000
32 0.000 0.000 0.352 0.000
33 0.000 0.000 0.300 0.000
34 0.000 0.000 0.468 0.000
35 0.000 0.000 0.300 0.000
36 0.000 0.000 0.300 0.000
37 0.000 0.000 0.300 0.000
38 0.000 0.000 0.300 0.000
39 0.000 0.000 0.300 0.000
Improper Coeffs
1 0.000 -1 2
2 0.000 -1 2
3 0.000 -1 2
4 0.000 -1 2
5 0.000 -1 2
6 0.000 -1 2
7 0.000 -1 2
8 0.000 -1 2
9 0.000 -1 2
10 0.000 -1 2
Atoms
1 1 1 -0.59270000 1.000 1.00000 0.00000
2 1 2 0.01280000 -0.421 1.00000 0.00000
3 1 3 -0.18580000 -0.913 1.00000 1.44051
4 1 4 -0.17550000 -2.438 0.99820 1.51208
5 1 5 -0.17760000 -2.928 1.00001 2.95981
6 1 6 -0.23910000 -4.447 0.99644 3.03271
7 1 7 0.40790000 1.284 0.99882 -0.93006
8 1 8 0.06810000 -0.768 0.10780 -0.53077
9 1 9 0.06810000 -0.768 1.89212 -0.53073
10 1 10 0.10490000 -0.509 1.87681 1.96163
11 1 11 0.10490000 -0.506 0.12459 1.96234
12 1 12 0.08900000 -2.829 0.11408 0.99353
13 1 13 0.08900000 -2.831 1.87968 0.99054
14 1 14 0.09080000 -2.545 1.88523 3.48082
15 1 15 0.09080000 -2.540 0.11893 3.48444
16 1 16 0.08140000 -4.864 1.88201 2.54268
17 1 17 0.08140000 -4.776 0.99909 4.07643
18 1 18 0.08140000 -4.859 0.10737 2.54537
Bonds
1 1 2 1
2 2 3 2
3 3 4 3
4 4 5 4
5 5 6 5
6 6 7 1
7 7 8 2
8 8 9 2
9 9 10 3
10 10 11 3
11 11 12 4
12 12 13 4
13 13 14 5
14 14 15 5
15 15 16 6
16 16 17 6
17 17 18 6
Angles
1 1 1 2 3
2 2 2 3 4
3 3 3 4 5
4 4 4 5 6
5 5 2 1 7
6 6 1 2 8
7 7 1 2 9
8 8 2 3 10
9 9 2 3 11
10 10 3 4 12
11 11 3 4 13
12 12 4 5 14
13 13 4 5 15
14 14 5 6 16
15 15 5 6 17
16 16 5 6 18
17 17 4 3 11
18 18 5 4 12
19 19 17 6 18
20 20 6 5 14
21 21 16 6 17
22 22 14 5 15
23 23 5 4 13
24 24 4 3 10
25 25 3 2 8
26 26 6 5 15
27 27 8 2 9
28 28 12 4 13
29 29 10 3 11
30 30 3 2 9
31 31 16 6 18
Dihedrals
1 1 4 3 2 1
2 2 5 4 3 2
3 3 6 5 4 3
4 4 7 1 2 3
5 5 16 6 5 4
6 6 10 3 2 9
7 7 11 3 2 8
8 8 8 2 3 4
9 9 11 3 2 1
10 10 18 6 5 14
11 11 17 6 5 15
12 12 18 6 5 4
13 13 16 6 5 14
14 14 12 4 3 10
15 15 13 4 3 10
16 16 11 3 2 9
17 17 11 3 4 5
18 18 15 5 4 12
19 19 16 6 5 15
20 20 17 6 5 14
21 21 14 5 4 13
22 22 9 2 3 4
23 23 17 6 5 4
24 24 10 3 2 8
25 25 13 4 3 11
26 26 18 6 5 15
27 27 15 5 4 3
28 28 12 4 5 6
29 29 9 2 1 7
30 30 15 5 4 13
31 31 14 5 4 12
32 32 8 2 1 7
33 33 13 4 5 6
34 34 10 3 2 1
35 35 12 4 3 11
36 36 13 4 3 2
37 37 10 3 4 5
38 38 12 4 3 2
39 39 14 5 4 3
Impropers
1 1 2 1 3 8
2 2 2 9 3 1
3 3 3 10 4 2
4 4 3 2 11 4
5 5 4 3 12 5
6 6 4 5 3 13
7 7 5 4 14 6
8 8 5 4 6 15
9 9 6 17 5 16
10 10 6 18 5 16
これをOVITO等で描画すると以下のような分子が得られる。
lammpsファイルからmolfile形式に書き換え
得られたlammpsファイルをもとに以下ののようにmolfileに書き換える。
#1-pentanol molfile
18 atoms
17 bonds
31 angles
39 dihedrals
10 impropers
Coords
1 1 1 0
2 -0.421 1 0
3 -0.913 1 1.44051
4 -2.438 0.9982 1.51208
5 -2.928 1.00001 2.95981
6 -4.447 0.99644 3.03271
7 1.284 0.99882 -0.93006
8 -0.768 0.1078 -0.53077
9 -0.768 1.89212 -0.53073
10 -0.509 1.87681 1.96163
11 -0.506 0.12459 1.96234
12 -2.829 0.11408 0.99353
13 -2.831 1.87968 0.99054
14 -2.545 1.88523 3.48082
15 -2.54 0.11893 3.48444
16 -4.864 1.88201 2.54268
17 -4.776 0.99909 4.07643
18 -4.859 0.10737 2.54537
Types
1 1
2 2
3 3
4 4
5 5
6 6
7 7
8 8
9 9
10 10
11 11
12 12
13 13
14 14
15 15
16 16
17 17
18 18
Charges
1 -0.5927
2 0.0128
3 -0.1858
4 -0.1755
5 -0.1776
6 -0.2391
7 0.4079
8 0.0681
9 0.0681
10 0.1049
11 0.1049
12 0.089
13 0.089
14 0.0908
15 0.0908
16 0.0814
17 0.0814
18 0.0814
Bonds
1 1 2 1
2 2 3 2
3 3 4 3
4 4 5 4
5 5 6 5
6 6 7 1
7 7 8 2
8 8 9 2
9 9 10 3
10 10 11 3
11 11 12 4
12 12 13 4
13 13 14 5
14 14 15 5
15 15 16 6
16 16 17 6
17 17 18 6
Angles
1 1 1 2 3
2 2 2 3 4
3 3 3 4 5
4 4 4 5 6
5 5 2 1 7
6 6 1 2 8
7 7 1 2 9
8 8 2 3 10
9 9 2 3 11
10 10 3 4 12
11 11 3 4 13
12 12 4 5 14
13 13 4 5 15
14 14 5 6 16
15 15 5 6 17
16 16 5 6 18
17 17 4 3 11
18 18 5 4 12
19 19 17 6 18
20 20 6 5 14
21 21 16 6 17
22 22 14 5 15
23 23 5 4 13
24 24 4 3 10
25 25 3 2 8
26 26 6 5 15
27 27 8 2 9
28 28 12 4 13
29 29 10 3 11
30 30 3 2 9
31 31 16 6 18
Dihedrals
1 1 4 3 2 1
2 2 5 4 3 2
3 3 6 5 4 3
4 4 7 1 2 3
5 5 16 6 5 4
6 6 10 3 2 9
7 7 11 3 2 8
8 8 8 2 3 4
9 9 11 3 2 1
10 10 18 6 5 14
11 11 17 6 5 15
12 12 18 6 5 4
13 13 16 6 5 14
14 14 12 4 3 10
15 15 13 4 3 10
16 16 11 3 2 9
17 17 11 3 4 5
18 18 15 5 4 12
19 19 16 6 5 15
20 20 17 6 5 14
21 21 14 5 4 13
22 22 9 2 3 4
23 23 17 6 5 4
24 24 10 3 2 8
25 25 13 4 3 11
26 26 18 6 5 15
27 27 15 5 4 3
28 28 12 4 5 6
29 29 9 2 1 7
30 30 15 5 4 13
31 31 14 5 4 12
32 32 8 2 1 7
33 33 13 4 5 6
34 34 10 3 2 1
35 35 12 4 3 11
36 36 13 4 3 2
37 37 10 3 4 5
38 38 12 4 3 2
39 39 14 5 4 3
Impropers
1 1 2 1 3 8
2 2 2 9 3 1
3 3 3 10 4 2
4 4 3 2 11 4
5 5 4 3 12 5
6 6 4 5 3 13
7 7 5 4 14 6
8 8 5 4 6 15
9 9 6 17 5 16
10 10 6 18 5 16
分子を格子状に配置するlammps実行ファイルを作成する
units real
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style cvff
pair_style lj/cut/coul/cut 12 #12
pair_modify mix geometric
special_bonds lj/coul 0.0 0.0 0.5
#kspace_style pppm 0.0001
variable N_A index 6.022140857e+23
variable mass_C index 12.001
variable mass_O index 15.9994
variable mass_H index 1.008
variable dens_1pentanol index 0.05 # g/cm^3 #https://ja.wikipedia.org/wiki/%E3%83%A1%E3%82%BF%E3%83%8E%E3%83%BC%E3%83%AB
variable mass_1pentanol index $(v_mass_C+v_mass_O+4*v_mass_H)
variable ndens_1pentanol index $(v_dens_1pentanol/v_mass_1pentanol*v_N_A*1e6*1e-30) # 1/A^3
variable a_1pentanol index $((1/v_ndens_1pentanol)^(1/3))
variable lx index 100
variable ly index 100
variable lz index 100
variable lx_1pentanol index $(floor(v_lx/v_a_1pentanol)*v_a_1pentanol)
variable ly_1pentanol index $(floor(v_ly/v_a_1pentanol)*v_a_1pentanol)
print "size: ${lx} x ${ly} x ${lz}"
print "size: ${ly_1pentanol} x ${ly_1pentanol} x ${lz}"
region box block 0 $(v_lx_1pentanol) 0 $(v_ly_1pentanol) 10 $(v_lz+10)
create_box 18 box bond/types 17 angle/types 31 dihedral/types 39 improper/types 10 extra/bond/per/atom 17 extra/angle/per/atom 31 extra/dihedral/per/atom 39 extra/improper/per/atom 10 extra/special/per/atom 15
include "pentanol.setting"
molecule 1pentanol 1pentanol.molfile #offset 6 5 7 3 2 # Toff 6 boff 5 aoff 7 doff 3 ioff 2
lattice sc ${a_1pentanol}
create_atoms 0 box mol 1pentanol 101
#change_box all z final 0 $(v_lz+20) units box
############################################
write_data data.init
ここで、[pentanol.setting]は以下のOPLS力場パラメータを記述している。
pair_coeff 1 1 0.17 3.12
pair_coeff 2 2 0.066 3.5
pair_coeff 3 3 0.066 3.5
pair_coeff 4 4 0.066 3.5
pair_coeff 5 5 0.066 3.5
pair_coeff 6 6 0.066 3.5
pair_coeff 7 7 0 0
pair_coeff 8 8 0.03 2.5
pair_coeff 9 9 0.03 2.5
pair_coeff 10 10 0.03 2.5
pair_coeff 11 11 0.03 2.5
pair_coeff 12 12 0.03 2.5
pair_coeff 13 13 0.03 2.5
pair_coeff 14 14 0.03 2.5
pair_coeff 15 15 0.03 2.5
pair_coeff 16 16 0.03 2.5
pair_coeff 17 17 0.03 2.5
pair_coeff 18 18 0.03 2.5
#Bond Coeffs
bond_coeff 1 320 1.41
bond_coeff 2 268 1.529
bond_coeff 3 268 1.529
bond_coeff 4 268 1.529
bond_coeff 5 268 1.529
bond_coeff 6 553 0.945
bond_coeff 7 340 1.09
bond_coeff 8 340 1.09
bond_coeff 9 340 1.09
bond_coeff 10 340 1.09
bond_coeff 11 340 1.09
bond_coeff 12 340 1.09
bond_coeff 13 340 1.09
bond_coeff 14 340 1.09
bond_coeff 15 340 1.09
bond_coeff 16 340 1.09
bond_coeff 17 340 1.09
#Angle Coeffs
angle_coeff 1 50 109.5
angle_coeff 2 58.35 112.7
angle_coeff 3 58.35 112.7
angle_coeff 4 58.35 112.7
angle_coeff 5 55 108.5
angle_coeff 6 35 109.5
angle_coeff 7 35 109.5
angle_coeff 8 37.5 110.7
angle_coeff 9 37.5 110.7
angle_coeff 10 37.5 110.7
angle_coeff 11 37.5 110.7
angle_coeff 12 37.5 110.7
angle_coeff 13 37.5 110.7
angle_coeff 14 37.5 110.7
angle_coeff 15 37.5 110.7
angle_coeff 16 37.5 110.7
angle_coeff 17 37.5 110.7
angle_coeff 18 37.5 110.7
angle_coeff 19 33 107.8
angle_coeff 20 37.5 110.7
angle_coeff 21 33 107.8
angle_coeff 22 33 107.8
angle_coeff 23 37.5 110.7
angle_coeff 24 37.5 110.7
angle_coeff 25 37.5 110.7
angle_coeff 26 37.5 110.7
angle_coeff 27 33 107.8
angle_coeff 28 33 107.8
angle_coeff 29 33 107.8
angle_coeff 30 37.5 110.7
angle_coeff 31 33 107.8
#Dihedral Coeffs
dihedral_coeff 1 -1.552 0 0 0
dihedral_coeff 2 1.3 -0.2 0.2 0
dihedral_coeff 3 1.3 -0.2 0.2 0
dihedral_coeff 4 -0.356 -0.174 0.492 0
dihedral_coeff 5 0 0 0.3 0
dihedral_coeff 6 0 0 0.3 0
dihedral_coeff 7 0 0 0.3 0
dihedral_coeff 8 0 0 0.3 0
dihedral_coeff 9 0 0 0.468 0
dihedral_coeff 10 0 0 0.3 0
dihedral_coeff 11 0 0 0.3 0
dihedral_coeff 12 0 0 0.3 0
dihedral_coeff 13 0 0 0.3 0
dihedral_coeff 14 0 0 0.3 0
dihedral_coeff 15 0 0 0.3 0
dihedral_coeff 16 0 0 0.3 0
dihedral_coeff 17 0 0 0.3 0
dihedral_coeff 18 0 0 0.3 0
dihedral_coeff 19 0 0 0.3 0
dihedral_coeff 20 0 0 0.3 0
dihedral_coeff 21 0 0 0.3 0
dihedral_coeff 22 0 0 0.3 0
dihedral_coeff 23 0 0 0.3 0
dihedral_coeff 24 0 0 0.3 0
dihedral_coeff 25 0 0 0.3 0
dihedral_coeff 26 0 0 0.3 0
dihedral_coeff 27 0 0 0.3 0
dihedral_coeff 28 0 0 0.3 0
dihedral_coeff 29 0 0 0.352 0
dihedral_coeff 30 0 0 0.3 0
dihedral_coeff 31 0 0 0.3 0
dihedral_coeff 32 0 0 0.352 0
dihedral_coeff 33 0 0 0.3 0
dihedral_coeff 34 0 0 0.468 0
dihedral_coeff 35 0 0 0.3 0
dihedral_coeff 36 0 0 0.3 0
dihedral_coeff 37 0 0 0.3 0
dihedral_coeff 38 0 0 0.3 0
dihedral_coeff 39 0 0 0.3 0
#Improper Coeffs
improper_coeff 1 0 -1 2
improper_coeff 2 0 -1 2
improper_coeff 3 0 -1 2
improper_coeff 4 0 -1 2
improper_coeff 5 0 -1 2
improper_coeff 6 0 -1 2
improper_coeff 7 0 -1 2
improper_coeff 8 0 -1 2
improper_coeff 9 0 -1 2
improper_coeff 10 0 -1 2
mass 1 15.999
mass 2 12.011
mass 3 12.011
mass 4 12.011
mass 5 12.011
mass 6 12.011
mass 7 1.008
mass 8 1.008
mass 9 1.008
mass 10 1.008
mass 11 1.008
mass 12 1.008
mass 13 1.008
mass 14 1.008
mass 15 1.008
mass 16 1.008
mass 17 1.008
mass 18 1.008
以上のコードを実行すると
以下のように1-pentanolが格子状に配置される
NPT計算結果
トラジェクトリ
密度の時間経過
