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

[Mark Abraham]

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