[gmx-users] Questions about itp files

[Justin A. Lemkul]

Andrew DeYoung wrote:
> Hi,
> 
> If you have time, I am wondering if you can help me understand some parts
> about .itp files.  Before I ask my questions, though, I will paste below a
> few relevant files (or parts of files).  My system consists of 1000 SPC/E
> water molecules in a cubic box of edge length 4.71 nm.  I am using the
> OPLS-AA force field.  I am running Gromacs 4.5.4.
> 
> Here is my topology file (water.top):
> ---------------------------------------------
> ; Include forcefield parameters
> #include "oplsaa.ff/forcefield.itp"
> 
> ; Include water topology
> #include "oplsaa.ff/spce.itp"
> 
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ;  i funct       fcx        fcy        fcz
>    1    1       1000       1000       1000
> #endif
> 
> ; Include topology for ions
> #include "oplsaa.ff/ions.itp"
> 
> [ system ]
> ; Name
> water
> 
> [ molecules ]
> ; Compound        #mols
> SOL              1000
> ---------------------------------------------
> 
> Here is forcefield.itp (located in the folder oplsaa.ff), omitting the
> reference citations listed in comments:
> ---------------------------------------------
> #define _FF_OPLS
> #define _FF_OPLSAA
> 
> ; This force field uses a format that requires Gromacs 3.1.4 or later.
> ;
> ; References for the OPLS-AA force field:
> ; ...
> 
> [ defaults ]
> ; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
> 1               3               yes             0.5     0.5
> 
> #include "ffnonbonded.itp"
> #include "ffbonded.itp"
> #include "gbsa.itp"
> ---------------------------------------------
> 
> Here is spce.itp (located in the folder oplsaa.ff):
> ---------------------------------------------
> [ moleculetype ]
> ; molname       nrexcl
> SOL             2
> 
> [ atoms ]
> ;   nr   type  resnr residue  atom   cgnr     charge       mass
>      1  opls_116   1    SOL     OW      1      -0.8476
>      2  opls_117   1    SOL    HW1      1       0.4238
>      3  opls_117   1    SOL    HW2      1       0.4238
> 
> #ifndef FLEXIBLE
> [ settles ]
> ; OW    funct   doh     dhh
> 1       1       0.1     0.16330
> 
> [ exclusions ]
> 1       2       3
> 2       1       3
> 3       1       2
> #else
> [ bonds ]
> ; i     j       funct   length  force.c.
> 1       2       1       0.1     345000  0.1     345000
> 1       3       1       0.1     345000  0.1     345000
> 
> [ angles ]
> ; i     j       k       funct   angle   force.c.
> 2       1       3       1       109.47  383     109.47  383
> #endif
> ---------------------------------------------
> 
> Here is the part of ffnonbonded.itp (located in the folder oplsaa.ff) that
> specifies the atomtypes called opls_116 and opls_117):
> ---------------------------------------------
> [ atomtypes ]
> ; full atom descriptions are available in ffoplsaa.atp
> ; name  bond_type    mass    charge   ptype          sigma      epsilon
> ; ...
> #ifdef HEAVY_H
> ; ...
> #else
> ; ...
>  opls_116   OW  8     15.99940    -0.820       A    3.16557e-01  6.50194e-01
>  opls_117   HW  1      1.00800     0.410       A    0.00000e+00  0.00000e+00
> #endif
> ; ...
> ---------------------------------------------
> In case it is helpful, I have posted the entire ffnonbonded.itp file here:
> http://www.andrew.cmu.edu/user/adeyoung/july20/ffnonbonded.itp
> 

Not necessary.  Everyone already has all of these files in the standard Gromacs 
distribution, along with most of what you've already posted :)

> If you have time, I have a few questions:
> 
> 1.  My water.top file invokes both the OPLS-AA force field and the SPC/E
> water model, using the commands #include "oplsaa.ff/forcefield.itp" and
> #include "oplsaa.ff/spce.itp", respectively.  But, in case of competition,
> which .itp file "wins"?  For example, spce.itp specifies atomic charges
> (-0.8476 for OW and 0.4238 for each of HW1 and HW2).  But, spce.itp also
> specifies atomtypes; spce.itp says that OW is opls_116 and that HW1 and HW2
> are each opls_117.  Now, opls_116 and opls_117 are also (and more
> explicitly) defined in ffnonbonded.itp.  In ffnonbonded.itp, atomic charges,
> among other things, are specified (-0.820 for OW and 0.410 for each of HW1
> and HW2).  So, two different files are specifying charges for the oxygens
> and hydrogens, and while the respective charges are similar, they are not
> identical.  Thus, given my water.top, spce.itp, and ffnonbonded.itp files,
> what will be the exact charges: will they be -0.8476/0.4238 for OW/HW, or
> will they be -0.820/0.410 for OW/HW?  Is there any way that I can know for
> sure?
> 

The charges in ffnonbonded.itp are generic charges for various functional 
groups.  They are not used by grompp, so the charges in spce.itp are what are 
used.  You can always confirm your input by running gmxdump on your .tpr file.

> 2.  If SPC/E water is a rigid model, then why does it appear that bond and
> angle force constants are specified in spce.itp?  Does the fact that these
> force constants are couched within #ifndef FLEXIBLE and #endif mean that
> they will be used only if the model is explicitly specified as flexible
> (which would be done, I think, by typing -DFLEXIBLE in my .mdp file(s), as
> described here: http://manual.gromacs.org/current/online/mdp_opt.html#pp)?
> How can I be sure that my waters are rigid?  Is it just that if I do not
> have -DFLEXIBLE in my .mdp file(s), then my waters are certainly rigid and
> thus the bond and angle force constants will be ignored?
> 

Right, #define blocks are not invoked unless the appropriate -D(whatever) is 
used.  This is a C preprocessor trick.

> 3.  In forcefield.itp, the following parameters appear:
> nbfunc = 1
> comb-rule = 3
> gen-pairs = yes
> fudgeLJ = 0.5
> fudgeQQ = 0.5 
> 
> 3a.  Page 129 of section 5.7.1 in the Gromacs 4.5.4 Manual says that the
> default value for gen-pairs is "no", but here in this OPLS-AA forcefield.itp
> file, gen-pairs has been automatically set to "yes".  This makes me wonder
> if I am thus doing something very non-standard.  What does setting gen-pairs
> to "yes" do?  Page 129 says, "gen-pairs is for pair generation. The default
> is 'no', i.e. get 1-4 parameters from the pairtypes list. When parameters
> are not present in the list, stop with a fatal error. Setting 'yes'
> generates 1-4 parameters that are not present in the pair list from normal
> Lennard-Jones parameters using fudgeLJ."  But, I do not think that I
> understand.  In general, is using gen-pairs = yes a good idea?  
> 

Force fields that use "gen-pairs = no" list all [pairtypes] explicitly.  If 
there's something missing, grompp raises a fatal error, as stated in the manual. 
  There should be a sentence right around what you've quoted that describes the 
behavior of OPLS-AA; it generates pair interactions by scaling, thus it needs to 
generate pairs itself (i.e., "gen-pairs = yes").  There's nothing to worry about 
in the forcefield.itp files.

> 3b.  What are fudgeLJ and fudgeQQ?  Page 129 of the manual says, "fudgeLJ is
> the factor by which to multiply Lennard-Jones 1-4 interactions, default 1."
> and "fudgeQQ is the factor by which to multiply electrostatic 1-4
> interactions, default 1."  What does this mean?  It sounds like the
> potentials/forces will be multipled by that factor, thereby increasing the
> potentials/forces (if fudge > 1) or decreasing the potentials/forces (if
> fudge < 1).  In the OPLS-AA forcefield.itp file, fudgeLJ = 0.5 and fudgeQQ =
> 0.5 are used.  What does this mean?   Does this mean that my
> potentials/forces are being reduced?  Can you give me any advice about if
> this is reasonable, or do you know of other documentation that may describe
> the gen-pairs and fudge directives in more detail?  
> 

A 1-4 interaction occurs between atoms that are connected via three bonds.  For 
water, this is irrelevant.  Different force fields base their potential 
functions on different assumptions and practical concerns.  For some force 
fields, LJ parameters work out quite nicely for most intermolecular interactions 
except those at very short range (i.e., exerted between bonded atoms a few bonds 
away) so they need to be scaled down or else the molecule goes blasting apart. 
It's sort of a hack (hence "fudgeLJ/QQ" for a "fudge factor" to be applied in 
only certain circumstances).

> 3c.  One more thought I have is that since this forcefield.itp file is
> located in the oplsaa.ff folder, does this mean that the parameters set in
> this forcefield.itp file are those specific to, or prescribed by, the
> OPLS-AA force field?  This file is located in the "share" directory
> ("/usr/local/gromacs/share/gromacs/top/oplsaa.ff"), and I am not allowed to
> change these files; they are read only.  (Although I could change them by
> copying them, altering the contents, and editing the #include commands to
> refer to a local, modified version.)  So this makes me think that, yes, I
> should just trust that the values of gen-pairs and fudge used are prescribed
> by, and appropriate for, the OPLS-AA force field.  Do you know if this is
> correct?  
> 

Leave the file alone.  The forcefield.itp correctly describes the necessary 
intrinsic parameters for using the force field.  If you modify it, you're not 
really using OPLS-AA, you're using some bonded and nonbonded parameters, but 
then hacking the potential function, which invalidates the whole thing.

> 4.  In ffnonbonded.itp, why are both sigma and epsilon set to zero for HW
> (opls_117)?  This seems to imply that, as far as Lennard-Jones interactions
> are concerned, the hydrogens on the waters don't exist.  Or, in other words,
> in the absence of charges, the hydrogens don't "feel" the hydrogens, the
> hydrogens don't "feel" the oxygens, and the oxygens don't "feel" the
> hydrogens.  In other words, the hydrogens interact with the world only via
> electrostatic (Coulombic) interactions.  Is this a correct interpretation?
> 

Correct.  Many force fields do this.

-Justin

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================