[gmx-users] Questions about itp files

[Andrew DeYoung]

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

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?

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?

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?  

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?  

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?  

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?

--

Thank you very much for your time!

Sincerely,

Andrew DeYoung
Carnegie Mellon University