search

jewettaij / moltemplate

A general cross-platform tool for preparing simulations of molecules and complex molecular assemblies
http://www.moltemplate.org
MIT License
257 stars 99 forks source link

Force field in the GROMACS topology format to LAMMPS format #75

Closed yunessalman closed 2 years ago

yunessalman commented 2 years ago

Dear Andrew,

How do I convert this Ethylene glycol topology and parameter set in GROMACS format to LAMMPS format?

So, I have this OPLS-AA file and the ethlyene glycol molecule file from ATB.

The set was obtained from this paper https://doi.org/10.1021/jp109914s

#define _FF_OPLS
#define _FF_OPLSAA

[ defaults ]
; nbfunc    comb-rule   gen-pairs   fudgeLJ fudgeQQ
       1            3         yes       0.5     0.5

[ atomtypes ]
; name  bond_type      mass    charge   ptype          sigma      epsilon
    CG      CG  6  12.01100    0.1650       A    3.50000e-01  2.76144e-01
    HG      HG  1   1.00800    0.0600       A    2.50000e-01  1.25520e-01
    OG      OG  8  15.99940   -0.7200       A    3.00000e-01  7.11280e-01
    HO      HO  1   1.00800    0.4350       A    0.00000e+00  0.00000e+00

[ bondtypes ]
;  i     j   func         b0         kb
  CG    CG      1    0.15290   224262.0
  CG    HG      1    0.10900   284512.0
  CG    OG      1    0.14300   267776.0
  OG    HO      1    0.09450   462750.4

[ angletypes ]
;  i     j     k   func      th0      cth
  CG    CG    OG      1   108.00   418.40
  CG    CG    HG      1   110.70   313.80
  HG    CG    HG      1   107.80   276.14
  HG    CG    OG      1   109.50   292.88
  CG    OG    HO      1   108.50   460.24

[ dihedraltypes ]
;   i    j     k     l   func  coefficients
   OG    CG    CG    OG  3    9.852  11.989   4.987 -26.828   0.000   0.000
   CG    CG    OG    HO  3   -0.141   5.591   3.156  -8.606   0.000   0.000

; Fourier coefficients
; angle    F1      F2      F3
; OCCO   16.264  -4.987  13.414
; CCOH    1.727  -3.156   4.303

Thank you for your time and support.

Best regards, Yunes Salman

jewettaij commented 2 years ago

Hi Yunes.

I don't know of an automatic way to convert GROMACS/ITP files into moltemplate (LT) format.

EDIT (6/10/2022): For the specific case of files obtained from the ATB service, I suggest downloading them in directly in moltemplate (LT) format instead. (To do that, when you download the file from the web site, select "LAMMPS format". The resulting files you receive will actually be in moltemplate/LT format. This is easier and safer than downloading the ATB files in GROMACS/ITP format and converting them later.)

For ITP files like this one, it should be very easy to write a program to convert them to moltemplate (LT) format automatically. I actually started to write a converter program, but its not finished. Currently the program I wrote is not general enough to handle all of the different kinds of information stored in ITP files. I would need help from someone who is more familiar with the details of the GROMACS ITP file format in order to produce such a general conversion tool that we be useful to people. If you know about Gromacs ITP files and want to help me do this, let me know.

Alternatively, it may be possible to use OpenBabel to convert some ITP files onto a LAMMPS DATA file. And then you could try using the "ltemplify.py" program included with moltemplate to convert that DATA file into an LT file. But I suspect that won't work in this case because the Bonds, Angles, and Dihedrals lookup atoms according to their "bond type", and LAMMPS data files don't store this kind of information.

Alternatively, you can try doing the conversion manually. As an example, I attempted to convert your ITP file example into moltemplate format by hand. I did this quickly by hand so there may be mistakes. But I tried to be clear about where there is ambiguity and the additional information we need.

# I DON'T KNOW WHAT THIS MEANS SO I'M IGNORING IT.
##define _FF_OPLS
##define _FF_OPLSAA
#[ defaults ]
#; nbfunc    comb-rule   gen-pairs   fudgeLJ fudgeQQ
#       1            3         yes       0.5     0.5

# WILL TRANSLATE THIS SECTION BELOW INTO MOLTEMPLATE (LT) FORMAT.
# THE "name" IS USED IN YOUR "Data Atoms" SECTION.
# THE "bond_type" IS USED TO LOOKUP BONDED INTERACTIONS.
# (I DON'T KNOW WHAT TO DO WITH COLUMNS MARKED WITH "?")
#
#[ atomtypes ]
#; name  bond_type      mass    charge   ptype          sigma      epsilon
#   CG      CG  6  12.01100    0.1650       A    3.50000e-01  2.76144e-01
#   HG      HG  1   1.00800    0.0600       A    2.50000e-01  1.25520e-01
#   OG      OG  8  15.99940   -0.7200       A    3.00000e-01  7.11280e-01
#   HO      HO  1   1.00800    0.4350       A    0.00000e+00  0.00000e+00
#
#              /|\                         /|\
#               |                           |
#               ?                           ?

write_once("Data Masses") {
  @atom:CG 12.01100
  @atom:HG 1.00000
  @atom:OG 15.99940
  @atom:HO 1.00000
}

# In moltemplate, the "full name" of each type of atom often contains
# both the original atom type "name" and "_b" followed by the "bond_type"
# (which, confusingly, is the same as the atom type name in this example).
# But you can refer to the atom type by its short name (eg "@atom:CG")
# in the "Data Atoms" section of your molecules.
# (That's what the "replace" command does.)

replace{ @atom:CG @atom:CG_bCG }
replace{ @atom:HG @atom:HG_bHG }
replace{ @atom:OG @atom:OG_bOG }
replace{ @atom:HO @atom:HO_bHO }

write_once("In Settings") {
  # The nonbonded Lennard-Jones parameters go in "In Settings"
  # which will end up in the "system.in.settings" file.
  # (Note: Gromacs and LAMMPS use different units by default.
  #        Epsilon values must be converted from kJ/mol to kcal/mol,
  #        and distances must be converted from nm to Angstroms.)
  # https://docs.lammps.org/pair_lj_cut_coul.html
  #            AtomType1   AtomType2     epsilon        sigma
  pair_coeff @atom:CG_bCG @atom:CG_bCG   0.0644594   3.50000
  pair_coeff @atom:HG_bHG @atom:HG_bHG   0.0292997   2.50000
  pair_coeff @atom:OG_bOG @atom:OG_bOG   0.1660317   3.00000
  pair_coeff @atom:HO_bHO @atom:HO_bHO   0.0000000   0.00000
}

# If you want to associate each type of atom with a charge, you
# can do that this way.  (This means the charges will not appear in
# the DATA file.  You will have to use "include system.in.charges"
# in your LAMMPS input script.

write_once("In Charges") {
  set type @atom:CG charge 0.1650
  set type @atom:HG charge 0.0600
  set type @atom:OG charge -0.7200
  set type @atom:HO charge 0.4350
}

# NOW I WILL TRANSLATE THIS SECTION INTO MOLTEMPLATE FORMAT:
#[ bondtypes ]
#;  i     j   func         b0         kb
#  CG    CG      1    0.15290   224262.0
#  CG    HG      1    0.10900   284512.0
#  CG    OG      1    0.14300   267776.0
#  OG    HO      1    0.09450   462750.4

write_once("In Settings") {
  # Assuming that "b0" is a bond distance, then this must
  # be converted from nm to Angstroms.
  # Furthermore, the spring constant "kb" must be converted
  # from kJ/mol/nm^2 to kcal/mol/Angstroms^2.
  # I attempted to do that below.
  # https://docs.lammps.org/bond_harmonic.html
  # BondTypeName   kb            r0
  @bond:CG_CG     536.000    1.5290 
  @bond:CG_HG     536.000    1.0900
  @bond:CG_OG     536.000    1.4300
  @bond:CG_HO     536.000    0.9450
}

# Names like "@bond:CG_HG" are just names.
# Moltemplate does not understand their meaning.
# So we have to generate rules so that moltemplate automatically knows to
# lookup bonds between pairs of atoms according to the atoms' "bond_type".
# In principle many different types of atoms can share the same bond_type
# although in this example that does not happen.
# (For example, there could be two different atom types "CGa" and "CGb"
#  use the same bond parameters, but have different Lennard Jones or charges.
#  In that case we would name them "@atom:CGa_bCG" and and "@atom:CGb_bCG".)
# So we use the "bond_type" (bCG) to lookup bonded interaction parameters
# and ignore the atom names ("CGa" or "CGb") when looking up bond parameters.
# (The * is a wildcard.)

write_once("Data Bonds By Type") {
  @bond:CG_CG  @atom:*_bCG   @atom:*_bCG
  @bond:CG_HG  @atom:*_bCG   @atom:*_bHG
  @bond:CG_OG  @atom:*_bCG   @atom:*_bOG
  @bond:CG_HO  @atom:*_bCG   @atom:*_bHO
}

# We do the same for the Angle interactions.

#[ angletypes ]
#;  i     j     k   func      th0      cth
#  CG    CG    OG      1   108.00   418.40
#  CG    CG    HG      1   110.70   313.80
#  HG    CG    HG      1   107.80   276.14
#  HG    CG    OG      1   109.50   292.88
#  CG    OG    HO      1   108.50   460.24
#
#                     /|\
#                      |
#                      ?

# I am guessing that the "cth" number is in units of kJ/mol.
# IF that's true, w need to convert it to units of kcal/mol.
write_once("In Settings")
  # https://docs.lammps.org/angle_harmonic.html
  # AngleTypeName    cth        th0
  @angle:CG_CG_OG    418.40    108.00
  @angle:CG_CG_HG    313.80    110.70 
  @angle:HG_CG_HG    276.14    107.80
  @angle:HG_CG_OG    292.88    109.50 
  @angle:CG_OG_HO    460.24    108.50
}

write_once("Data Angles By Type") {
  @angle:CG_CG_OG   @atom:*_bCG   @atom:*_bCG   @atom:*_bOG
  @angle:CG_CG_HG   @atom:*_bCG   @atom:*_bCG   @atom:*_bHG
  @angle:HG_CG_HG   @atom:*_bHG   @atom:*_bCG   @atom:*_bHG
  @angle:HG_CG_OG   @atom:*_bHG   @atom:*_bCG   @atom:*_bOG
  @angle:CG_OG_HO   @atom:*_bCG   @atom:*_bOG   @atom:*_bHO
}

# Let's try to do the same thing for the dihedrals.
# I provide only a sketch how to do this below, but I left some numbers ???.
# I don't know what equation they are using (for "func 3").
# You will need to lookup what "func 3" means, and 
# CONVERT ALL PARAMETERS WITH UNITS OF ENERGY FROM kj/mol TO kcal/mol.
# I left these parameters as ??? below.
#
# [ dihedraltypes ]
#;   i    j     k     l   func  coefficients
#   OG    CG    CG    OG  3    9.852  11.989   4.987 -26.828   0.000   0.000
#   CG    CG    OG    HO  3   -0.141   5.591   3.156  -8.606   0.000   0.000
#
#; Fourier coefficients
#; angle    F1      F2      F3
#; OCCO   16.2    64  -4.987  13.414
#; CCOH    1.727  -3.156   4.303

write_once("Data Dihedrals By Type") {
  @dihedral:OG_CG_CG_OG   @atom:*_bOG  @atom:*_bCG  @atom:*_bCG  @atom:*_bOG
  @dihedral:CG_CG_OG_HO   @atom:*_bCG  @atom:*_bCG  @atom:*_bOG  @atom:*_bHO
}

write_once("In Settings") {
  @dihedral:OG_CG_CG_OG    ???   ???   ???   ???   ???   ???   ???
  @dihedral:CG_CG_OG_HO    ???   ???   ???   ???   ???   ???   ???
}

# If you had improper interactions in your molecule, then it gets more
# complicated because different MD file formats disagree on the order
# of the atoms.  (For example: Should the central atom in an improper
# interaction in an ITP file appear first in the list of atoms?  Not always.
# In some file formats, I think it the central atom is listed 2nd or 3rd atom.)
# So someone will have to look up the details of the ITP file format to figure
# out how it handles this ambiguity.
jewettaij commented 2 years ago

To clarify: I am happy to write the force field converter, but I would need someone I could ask questions who knows about the ITP format.

If anyone is curious, the format used by moltemplate to store force field information is explained in more detail here, and also demonstrated in these examples. (You have to click on the "forcefield.lt" file and open it in a text editor.)

yunessalman commented 2 years ago

Hi Andrew,

Thank you so much for your time and effort, I appreciate it. Also, thank you for moltemplate it is a very helpful tool.

Regarding Gromacs ITP files, unfortunately, I've never worked with Gromacs before so I don't know anything about their file structure. I hope you can find someone who can help you with this.

jewettaij commented 2 years ago

I just re-read your original post. Are you aware that the ATB service can provide files directly in moltemplate format? (When you download the files, select "LAMMPS" format. This is actually moltemplate/LT format.)