[gmx-users] potential energy of one molecule using GAFF: confused regarding Ryckaert-Bellemans and 1-4 potentials
Vedat Durmaz
durmaz at zib.de
Tue Oct 18 15:34:37 CEST 2011
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
More information about the gromacs.org_gmx-users
mailing list