[gmx-users] potential energy of one molecule using GAFF: confused regarding Ryckaert-Bellemans and 1-4 potentials

[Vedat Durmaz]

hey guys,

i'm a really confused at the moment. reading the gromacs manual 4.5.4 
(GM) and the paper introducing GAFF did not provide insight into this 
issue which drives me crazy! maybe some advanced gromacs user can do so.

up to now, i have simulated very few different systems using amber force 
fields and gromacs (4.0.7). assume you have simulated the global minimum 
conformation of some handsome molecule "mol" mainly consisting of a 
simple cyclic aliphatic structure (about 13 atoms long) containing a 
double bond resulting in an E-mol and Z-mol isomer.

i know from radicalic interconversion experiments, that Z-mol is about 
100.000 times likelier than E-mol (due to the cyclic structure).

from the simulation, which i've performed in water as well as in vacuum, 
i find that my intramolecular LJ14 and Coul14 energies indeed reflect 
the experimental steady state distribution. LJ-SR and Coul-SR just show 
slight differences only for the isomers. now, in order to see, whether 
the high 14-difference is caused by solvent effects or not, i extracted 
several energies from the vacuum simulation (E-mol left and Z-mol right) 
allowing me additionally to extract the potential energy (as well as 
bonded terms) of the molecule under consideration only:

1  Bond              69.2    67.3

2  Angle             96.1    104.9

3  Proper Dih.        3.1    4.6

4  Ryckaert-Bell.    63.9    106.3

5  LJ-14             48.4    46.1

6  Coulomb-14     -1036.6    -1065.3

7  LJ (SR)           20.6    19.1

8  Coulomb (SR)     555.9    556.5

9  Potential       -179.4    -160.5


the potential energy (line 9) obviously equals the sum of all other 
terms (line 1-8) and reveals, in contrast to the Coulomb-14 (C14) term, 
a completely different steady state distribution of cis- and trans-mol. 
this is mainly due to the Ryckaert-Bellemans (RB) term.

however, as i've understood from the GM, the RB potential for dihedrals 
is not incorporated by every force field because often, the "periodic 
type" is used instead. and indeed, when comparing the equations from the 
gromacs manual with the Wang-paper in jCompChem-2004 concerning GAFF, it 
seems neither GAFF uses RB.


and here are my lovely questions:

Q1 am i right with my assumption whether GAFF uses RB or not? if not, 
will someone please explain it to me in a few words?

Q2 why does the md.edr file provide two different terms concerning the 
same thing (proper dihedrals): "proper dihedrals" and "RB" how can i 
interpret these values?


GM says: "With the periodic GROMOS potential a special 1-4 LJ-interaction
must be included; with the Ryckaert-Bellemans potential for alkanes the 
1-4 interactions must be
excluded from the non-bonded list. Note: Ryckaert-Bellemans potentials 
are also used in e.g. the
OPLS force field in combination with 1-4 interactions"

Q3 does anyone know, how this is handled in GAFF?

Q4 does this mean i have to exclude 1-4 manually, after having 
parameterized my molecules using acpype (antechamber)?


GM says: "For alkanes, the following proper dihedral potential [RB] is 
often used"

Q5 what does it depend on, whether "RB" is used for alkanes or not?


GM says: "For third neighbors, the normal Lennard-Jones repulsion is 
sometimes still too strong, ... especially the case for carbon-carbon 
interactions in a cis-conformation ... for some of these interactions, 
the Lennard-Jones repulsion has been reduced in the GROMOS force field, 
which is implemented by keeping a separate list of 1-4 and normal 
Lennard-Jones parameters. In other force fields, such as OPLS [88], the 
standard Lennard-Jones parameters are reduced by a factor of two, ..."

Q6 how this is handled with GAFF in gromacs?

Q7 in general, can i "still trust" the given energies as listed above 
and compare them quantitatively? the simulation should actually result 
in a favorable Z-mol (due to strain) and less favorable E-mol. are the 
force fields and simulation software able to capture these kind of 
differences concerning the steady-state distribution of cis/trans, when 
the respective potentials are modified in such various ways among the 
force fields?

Q8 does anyone have an idea, how to perform the simulation and on which 
energy terms to concentrate in order to get reliable results?

thanks in advance and sorry for the long text ... as i've already told 
you, i'm really confused.

kind regards

vedat durmaz