lammpsで混合液系を作成

## packmolのインストール ##
- packmolのgithubリポジトリからソースコードをダウンロード
- Usageを見ながらコンパイル
  
https://github.com/m3g/packmol/tree/v20.14.3

- 使用できるか確認
```
$packmol

################################################################################

 PACKMOL - Packing optimization for the automated generation of
 starting configurations for molecular dynamics simulations.
 
                                                              Version 20.010 

################################################################################

  Packmol must be run with: packmol < inputfile.inp 

  Userguide at: http://m3g.iqm.unicamp.br/packmol 

  Reading input file... (Control-C aborts)

  
```
## moltemplateのインストール ##
- moltemplateをgithubからダウンロード

https://github.com/jewettaij/moltemplate/blob/master/README.md#installation-example


- 以下のサイトを参考にbashファイルを書き換えてmoltemplateのPATHを通す　


https://github.com/jewettaij/moltemplate/blob/master/INSTALL.md

- 使用できるか確認
```
$moltemplate.sh

moltemplate.sh v2.20.20 2023-4-20

lttree_check.py v0.81.2 2021-5-24
########################################################
##            WARNING: atom_style unspecified         ##
## --> "Data Atoms" column data has an unknown format ##
##              Assuming atom_style = "full"          ##
########################################################

Error: This program requires at least one argument
       the name of a file containing ttree template commands
```
## packmolで混合液系の系を作成 ##
- ディレクトリを作成
以降、ファイルはこの場所に保存する。
```
$mkdir test
$cd test
```


- 水のPDBファイルを次のようにする。
```
#water.pdb
ATOM      1  O   UNK   1       5.901   7.384   1.103       
ATOM      2  H1   UNK   1       6.047   8.238   0.581        
ATOM      3  H2   UNK   1       6.188   7.533   2.057 
```
![water.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/2eb952de-a891-650f-8209-ca3d2f3139e1.png)
##

- ligpargenで有機分子のPDBファイルを作成する。
今回はグリセロールを作成した。
```
#glycerol.pdb
REMARK LIGPARGEN GENERATED PDB FILE
ATOM      1  C00 UNK     1       1.000   1.000   0.000
ATOM      2  O01 UNK     1      -0.427   1.000   0.000
ATOM      3  C02 UNK     1       1.530   1.000  -1.432
ATOM      4  O03 UNK     1       0.924  -0.109  -2.117
ATOM      5  C04 UNK     1       1.202   2.299  -2.175
ATOM      6  O05 UNK     1       1.733   2.220  -3.505
ATOM      7  H06 UNK     1       1.324   0.098   0.528
ATOM      8  H07 UNK     1       1.350   1.876   0.558
ATOM      9  H08 UNK     1      -0.673   0.205  -0.531
ATOM     10  H09 UNK     1       2.610   0.835  -1.432
ATOM     11  H0A UNK     1       0.684   0.212  -3.011
ATOM     12  H0B UNK     1       1.641   3.174  -1.678
ATOM     13  H0C UNK     1       0.127   2.441  -2.262
ATOM     14  H0D UNK     1       2.702   2.118  -3.439
TER 
CONECT    1    2 
CONECT    1    3 
CONECT    3    4 
CONECT    3    5 
CONECT    5    6 
CONECT    1    7 
CONECT    1    8 
CONECT    2    9 
CONECT    3   10 
CONECT    4   11 
CONECT    5   12 
CONECT    5   13 
CONECT    6   14 
END
```
![glycerol.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/519c64c5-0161-b1ee-c052-c085aedfff12.png)
##

- packmolのスクリプトファイルを作成する
```
#water_glycerol.inp
tolerance 2.0                #隣り合う分子の最小距離
filetype pdb                 #インプットファイルの形式
output out.pdb               #出力ファイル名
 
structure glycerol.pdb       #混合液成分①のファイル名
  number 500                 #用意する分子数
  inside box 0 0 0 50 50 50  #分子を配置する枠を定義
end structure

structure water.pdb          #混合液成分②のファイル名
  number 5000                #
 inside box 0 0 0 50 50 50   #
end structure
```

- スクリプトを実行
```
$packmol < water_glycerol.inp

################################################################################

 PACKMOL - Packing optimization for the automated generation of
 starting configurations for molecular dynamics simulations.
 
                                                              Version 20.010 

################################################################################

  Packmol must be run with: packmol < inputfile.inp 

  Userguide at: http://m3g.iqm.unicamp.br/packmol 

  Reading input file... (Control-C aborts)
  Seed for random number generator:      1234567
  Output file: out.pdb
  Reading coordinate file: glycerol.pdb
  Reading coordinate file: water.pdb
  Number of independent structures:            2
  The structures are: 
  Structure            1 :glycerol.pdb(          14  atoms)
  Structure            2 :water.pdb(           3  atoms)
  Maximum number of GENCAN loops for all molecule packing:          400
  Total number of restrictions:            2
  Distance tolerance:    2.0000000000000000     
  Warning: Type of residue numbering not set for structure            1
  Residue numbering set for structure            1 :           0
  Swap chains of molecules of structure            1 : F
  Warning: Type of residue numbering not set for structure            2
  Residue numbering set for structure            2 :           0
  Swap chains of molecules of structure            2 : F
  Number of molecules of type            1 :          500
  Number of molecules of type            2 :         5000
  Total number of atoms:        22000
  Total number of molecules:         5500
  Number of fixed molecules:            0
  Number of free molecules:         5500
  Number of variables:        33000
  Total number of fixed atoms:            0
  Maximum internal distance of type            1 :    4.8483665290487270     
  Maximum internal distance of type            2 :    1.6417923132966601     
  All atoms must be within these coordinates: 
   x: [   -997.94225253270474      ,    1002.0577474672953       ] 
   y: [   -997.47854533676423      ,    1002.5214546632358       ] 
   z: [   -996.79999999999995      ,    1003.2000000000000       ] 
  If the system is larger than this, increase the sidemax parameter. 

################################################################################

~~~
~~~
~~~

  Starting GENCAN loop:           22
  Scaling radii by:    1.0000000000000000     

  Packing:|0                                                             100%|
          |******************************************************************|
          |*********************

  Function value from last GENCAN loop: f = .19953E-01
  Best function value before: f = .60780E-01
  Improvement from best function value:    67.17 %
  Improvement from last loop:    67.17 %
  Maximum violation of target distance:     0.008874
  Maximum violation of the constraints: .47777E-02

################################################################################

                                 Success! 
              Final objective function value: .19953E-01
              Maximum violation of target distance:   0.008874
              Maximum violation of the constraints: .47777E-02

--------------------------------------------------------------------------------

              Please cite this work if Packmol was useful: 

           L. Martinez, R. Andrade, E. G. Birgin, J. M. Martinez, 
         PACKMOL: A package for building initial configurations for
                   molecular dynamics simulations. 
          Journal of Computational Chemistry, 30:2157-2164,2009.

################################################################################

  Solution written to file: out.pdb

--------------------------------------------------------------------------------

   Running time:    84.0848160      seconds. 

--------------------------------------------------------------------------------
```
![water_glycerol.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/aabc4464-afd6-636a-153e-09cbd85ebda1.png)


## moltemplateによってスクリプトファイルを作成 ##
### OPLS力場における各原子のタイプを確認 ###
water.pdbおよびglycerol.pdbで定義した分子についてmoltemplate/moltemplate/force_fields/oplsaa.lt  
を確認しながら、各原子のタイプを確認していく。

```
#oplsaa.lt
~~~
~~~
    set type @atom:70 charge 0.0  # "TIP4F Water O"
    set type @atom:71 charge 0.511  # "TIP4F Water H"
    set type @atom:72 charge -1.022  # "TIP4F Water M"
    set type @atom:73 charge 0.0  # "TIP5P Water O"
    set type @atom:74 charge 0.241  # "TIP5P Water H"
    set type @atom:75 charge -0.241  # "TIP5P Water LP"
    set type @atom:76 charge -0.82  # "SPC Water O"
    set type @atom:77 charge 0.41  # "SPC Water H"
    set type @atom:78 charge -1.02  # "Ammonia NH3"
    set type @atom:79 charge 0.34  # "Ammonia NH3"
    set type @atom:80 charge -0.18  # "Alkane CH3-"
    set type @atom:81 charge -0.12  # "Alkane -CH2-"
    set type @atom:82 charge -0.06  # "Alkane >CH-"
    set type @atom:83 charge -0.24  # "Methane CH4"
    set type @atom:84 charge 0.0  # "Alkane >C<"
    set type @atom:85 charge 0.06  # "Alkane H-C"
    set type @atom:86 charge 0.0  # "Alkene R2-C="
    set type @atom:87 charge -0.115  # "Alkene RH-C="
    set type @atom:88 charge -0.23  # "Alkene H2-C="
    set type @atom:89 charge 0.115  # "Alkene H-C="
    set type @atom:90 charge -0.115  # "Aromatic C"
    set type @atom:91 charge 0.115  # "Aromatic H-C"
    set type @atom:92 charge 0.0  # "Naphthalene Fusion C"
    set type @atom:93 charge -0.065  # "Ethyl Benzene CH3-"
    set type @atom:94 charge -0.005  # "Ethyl Benzene -CH2-"
    set type @atom:95 charge -0.115  # "Diene =CH-CH="
    set type @atom:96 charge -0.683  # "Alcohol -OH"
    set type @atom:97 charge 0.418  # "Alcohol -OH"
    set type @atom:98 charge 0.04  # "Methanol CH3-"
    set type @atom:99 charge 0.145  # "Alcohol CH3OH & RCH2OH"
    set type @atom:100 charge 0.205  # "Alcohol R2CHOH"
~~~
~~~
```

- water.pdbで定義した原子IDはO, H1,H2 の３つであり、SPC力場パラメータを割り当てたいので、Oには@atom:76、H1,H2は@atom:77を使用する。  



- glycerol.pdbで定義した分子については以下の画像のように割り当てる
原子がどのような結合の一部かで原子のタイプを選択する

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/dd812d1f-92f3-733f-fb58-2833bf69e73a.png)


### moltemplateの形式で各分子を定義する ###
基本のltファイル形式

```
# Use the OPLS-AA force field for all species.
import /home/username/moltemplate/moltemplate/force_fields/oplsaa.lt


#MOLECULE1はnoltemplateを使う上での分子の名前を表す
MOLECULE1 inherits OPLSAA {
    #pdbファイル内で定義された原子を全て再定義する。
    write("Data Atoms"){ 
       #$atom:ID    $mol:ID     @atom:atomtype charge xcood ycood zcood 
        $atom:C00	$mol:GLY	@atom:99       0      1	    1	  0 
        ~~~
        ~~~
        ~~~
    }

    write("Data Bond List") {
       #pdbファイル内で定義されたCONNET部分から結合を定義する。
       #$bond:ID    $atom:ID    $atom:ID
        $bond:CO01	$atom:C00	$atom:O01
        ~~~
        ~~~
        ~~~
    }
}

```

##
水のltファイル（water.lt）
```
# Use the OPLS-AA force field for all species.
import /home/username/moltemplate/moltemplate/force_fields/oplsaa.lt

Water inherits OPLSAA {
  write("Data Atoms"){
    $atom:O  $mol:Water @atom:76 0.  5.901   7.384   1.103
    $atom:H1 $mol:Water @atom:77 0.  6.047   8.238   0.581
    $atom:H2 $mol:Water @atom:77 0. 6.188   7.533   2.057
  }
  write("Data Bond List"){
    $bond:OH1 $atom:O $atom:H1
    $bond:OH2 $atom:O $atom:H2
  }
}

```

##
グリセロールのltファイル(glycerol.lt)
```

# Use the OPLS-AA force field for all species.
import /home/username/moltemplate/moltemplate/force_fields/oplsaa.lt



Glycerol inherits OPLSAA {
    write("Data Atoms"){
        $atom:C00	$mol:GLY	@atom:99	0	1	1	0
        $atom:O01	$mol:GLY	@atom:96	0	-0.427	1	0
        $atom:C02	$mol:GLY	@atom:100	0	1.53	1	-1.432
        $atom:O03	$mol:GLY	@atom:96	0	0.924	-0.109	-2.117
        $atom:C04	$mol:GLY	@atom:99	0	1.202	2.299	-2.175
        $atom:O05	$mol:GLY	@atom:96	0	1.733	2.22	-3.505
        $atom:H06	$mol:GLY	@atom:85	0	1.324	0.098	0.528
        $atom:H07	$mol:GLY	@atom:85	0	1.35	1.876	0.558
        $atom:H08	$mol:GLY	@atom:97	0	-0.673	0.205	-0.531
        $atom:H09	$mol:GLY	@atom:85	0	2.61	0.835	-1.432
        $atom:H0A	$mol:GLY	@atom:97	0	0.684	0.212	-3.011
        $atom:H0B	$mol:GLY	@atom:85	0	1.641	3.174	-1.678
        $atom:H0C	$mol:GLY	@atom:85	0	0.127	2.441	-2.262
        $atom:H0D	$mol:GLY	@atom:97	0	2.702	2.118	-3.439
    }

    write("Data Bond List") {
        $bond:CO01	$atom:C00	$atom:O01
        $bond:CC02	$atom:C00	$atom:C02
        $bond:CO03	$atom:C02	$atom:O03
        $bond:CC04	$atom:C02	$atom:C04
        $bond:CO05	$atom:C04	$atom:O05
        $bond:CH06	$atom:C00	$atom:H06
        $bond:CH07	$atom:C00	$atom:H07
        $bond:OH08	$atom:O01	$atom:H08
        $bond:CO09	$atom:C02	$atom:H09
        $bond:OH10	$atom:O03	$atom:H0A
        $bond:CH11	$atom:C04	$atom:H0B
        $bond:CH12	$atom:C04	$atom:H0C
        $bond:OH13	$atom:O05	$atom:H0D
    }
}
```
##
### packmolによるpdbファイルにOPLS力場を割当てる ###
割り当て用のスクリプトファイル（water_glycerol.lt）を作成する。

基本の形式
```
import "MOLECULE1.lt"
import "MOLECULE2.lt"

# The lipids and water must be listed instantiated in the same order
# they appear in the packmol_files/mix_lipids+water.inp file:

MOL1 = new MOLECULE1[500]   #ltファイルで定義した分子名、packmolで作成した分子数に合わせる

MOL2 = new MOLECULE2[5000]  #複数の分子にそれぞれ割り当てる場合、上の行からパラメータ割当てがなされるため、packmolスクリプトで書いた順番に合わせる。


write_once("Data Boundary") {  #packmolで定義した枠を踏襲する
  0 50.0 xlo xhi
  0 50.0 ylo yhi
  0 50.0 zlo zhi
}
```

- 水とグリセロールの混合系割当てスクリプト（water + glycerol.lt）

```
import "water.lt"
import "glycerol.lt"

# The lipids and water must be listed instantiated in the same order
# they appear in the packmol_files/mix_lipids+water.inp file:

GLYs = new Glycerol[500]

waters = new Water[5000]


write_once("Data Boundary") {
  0 50.0 xlo xhi
  0 50.0 ylo yhi
  0 50.0 zlo zhi
}
```

- スクリプトを実行

```
$moltemplate.sh -pdb out.pdb water_glycerol.lt
```

ここまででpackmolで作成したPDBファイルにOPLS力場を割当てたlammps.dataファイルが作成される。
![water_glycerol_data.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2371660/3fa8f94e-95b7-b7a4-48ed-3050f86fd8a7.png)

※packmolでぎゅうぎゅうに詰め込みすぎるとNVT計算した途端に分子が吹っ飛んでしまうので注意
（私は吹っ飛びました）
## 参考資料・サイト ##
- lammps documentにmoltemplateについて記載がある。まずはここを見るのがいいと思う。
https://docs.lammps.org/latest/Howto_moltemplate.html
- 以下のサイト内の「Coarse-grained lipid bilayer using the MARTINI force field and PACKMOL」という例題が参考になる。
https://www.moltemplate.org/visual_examples.html
- 作成したファイル群は以下からダウンロードできる　
https://github.com/naonao07/my_lammps_examples/tree/main/water_glycerol