有機分子の分子動力学シミュレーション方法

# 有機分子の分子動力学計算 #

## 有機分子の描画と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等で描画すると以下のような分子が得られる。


![1-pentanol.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/91ed6d27-4a72-c929-f63b-4e8f305175d2.png)

## 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が格子状に配置される
![mksys_pentanol.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/c3876686-d7cd-c4d3-099c-5f1526a00f2e.png)

## NPT計算結果 ##
トラジェクトリ
![npt_pentanol.gif](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/35857ebc-3f0f-2714-5c02-9eda85f8f005.gif)

密度の時間経過
![density_pentanol.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/271d62eb-074c-8ba0-293b-a407dca1bb33.png)