2
0

lammpsで混合液系を作成

Last updated at Posted at 2024-03-03

packmolのインストール

  • packmolのgithubリポジトリからソースコードをダウンロード
  • Usageを見ながらコンパイル

  • 使用できるか確認
$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からダウンロード

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

  • 使用できるか確認
$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

  • 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

  • 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

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 の3つであり、SPC力場パラメータを割り当てたいので、Oには@atom:76、H1,H2は@atom:77を使用する。

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

image.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

※packmolでぎゅうぎゅうに詰め込みすぎるとNVT計算した途端に分子が吹っ飛んでしまうので注意
(私は吹っ飛びました)

参考資料・サイト

2
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
2
0