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
========================================
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
