On 19/10/2011 12:34 AM, Vedat Durmaz wrote:
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).
Don't exclusively read documentation that doesn't match the version of
the software. Things change.
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).
Part of that will be energy-based (E-mol has higher intra-atomic
repulsion), and part will be entropy-based (Z-mol has fewer available
configurations), but apparently a 13-atom ring is large enough that we
should see the former dominate.
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.
Maybe that is just fortuitous. Force fields are not parameterized to
necessarily be able to be decomposed in ways that correlate numbers with
chemical intuition. In particular, there is a balance between 1-4
charge-charge interaction and 1-4 dihedral interaction which depends on
how everything else is parameterized.
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,
Note that you can construct subsets from your water .tpr and .trr with
tpbconv and trjconv, and then use mdrun -rerun on these existing
configuration, rather than run a vacuum simulation that samples
different space.
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:
I'm assuming these are post-equilibration average values over at least
several nanoseconds. Otherwise they don't mean anything at all.
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.
One could equally observe that the contribution to the "wrong" relative
energy is about half due to the angle term and about half due to the sum
of dihedral and 1-4 terms.
Apparently your experimental results suggest the same relative stability
trend should be seen with the vacuum and condensed-phase systems,
however the parameters you are using for the simulation are optimized
for the condensed phase, so it is far from clear that they will produce
the correct ensemble.
If you'd done the mdrun -rerun I suggested and produced numbers like the
above, you'd have the problem that you'd have assumed that the solvent
potential energy is identical in each case.
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?
Read the GAFF paper for what mathematical function should be used.
GROMACS would then have the job of implementing that function, but it
need not do so with one that is named the same way. In particular,
before GROMACS 4.5, there was no way of implementing multiple periodic
dihedral types on the same four atoms (as used by CHARMM and AMBER force
fields), so these had to be combined into one equivalent RB form.
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?
Your .top file has both types, and it seems the latter result from the
above-mentioned combining.
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?
I think these observations are force-field-specific. I believe AMBER
force fields for 4.0.x used R-B as a hack. They should handle 1-4
interactions the way the AMBER force field documentation specifies.
Consult your .top file to see what it is doing.
Q4 does this mean i have to exclude 1-4 manually, after having
parameterized my molecules using acpype (antechamber)?
I doubt it, but now you have a third layer of documentation to get straight.
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?
Probably the force field.
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?
Not relevant, as above.
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?
The GROMACS manual has to attempt to describe the kinds of things it
does for a dozen different force fields. It cannot do so exhaustively
for every force field. It does compute them faithfully, but nomenclature
and exceptions must change from case to case. Ultimately, you have to
have confidence that your knowledge of chapter 5 of the manual means
that you know that the contents of your .top match how the force field
literature specifies the maths should work. General statements in (say)
chapter 4 cannot be specifically accurate for all cases. So you should
not pay strong attention to statements that don't pertain to your case,
and you should read documentation that matches your case.
Q8 does anyone have an idea, how to perform the simulation and on
which energy terms to concentrate in order to get reliable results?
You seem to be performing it OK, given that you've said very little
about any details...
I think the problem is poorly constructed. You have some experimental
data that gives a general understanding of the size of the free energy
difference between the isomers. You can't necessarily expect to
reproduce that from (average) potential energy differences between
conformations of those isomers. Measuring the free energy difference
with simulations is hard because you cannot model the intermediate
stages of bond breaking and forming. There are "alchemical" free energy
methods that could in principle treat this problem effectively, but
there will be some significant issues and you are best doing your own
homework there.
Mark
--
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
