[gmx-users] Gromacs Relative Free Energy Calculation Issue

Hannes Loeffler Hannes.Loeffler at stfc.ac.uk
Tue Oct 4 10:51:48 CEST 2016


On Mon, 3 Oct 2016 19:57:07 +0000
Guanglin Kuang <guanglin at kth.se> wrote:

> Dear Gromacs users,
> 
> Has any of you managed to use Gromacs to do relative free energy
> calculations?  I have some technical questions that would need your
> suggestions.
> 
> I am trying to reproduce the Amber free energy calculation using
> PMEMD (http://ambermd.org/tutorials/advanced/tutorial9/#overview).

Keep in mind that this is just a tutorial for demonstration purposes
only.  It is not intended to produce absolutely correct results.


> First, I generated the topology and coordinate files of the ligands
> (benzene and phenol) using alchemistry_setup.py, developed by Dr.
> David Mobley (https://github.com/MobleyLab/alchemical-setup), as well
> as antechamber with AM1/BCC charge.

I believe David's tool requires you to provide a mapping (in fact, the
idea is to use Lomap) but at first glance your topology seems fine.
That's why I have suggested to have a look into FESetup which can do
that all for you.  Current downside is that it supports AMBER force
fields only.


> Then I set up the mdp file and topology file for relative free energy
> calculations. I used the single topology together with the one-step
> and three-step transformations, respectively. Please have a look at
> the relevant files  attached for the details.
> 
> Finally, I used g_bar to calculate the free energy
> (alchemistry_analysis.py can give similar results). From the results
> I found that one-step transformation can give good results, only if
> some restraints are applied to the protein-ligand complex, because I
> found that benzene/phenol are not strong binders to lyszome, and
> would drift around in the binding pocket, thus give poor results
> deviating from experimental data.

You should not need to use restraints for this type of relative AFE
simulation.  If you still do I would think you would have to take
the energy cost of the restraint into account as well.  The ligands may
be weak binders but you should see the same effect in other codes as
well.  The tutorial is run only for a very short time.  Making it
converge may have a large effect on the final ddG too.


> However, in the three-step  strategy, the results are much worse,
> this is probably due to the fact that a dual topology is used in
> Amber but a single topology is used in Gromacs (my case), we may not
> be able to do the changes directly (as shown in the attached topology
> files).

The implementation details here should not matter.  With AMBER you can
and you obviously did map atoms which at the end means that you
effectively use a "single topology" (AMBER does energy scaling while
Gromacs may do parameter scaling but I would have to check).  The real
difference is that AMBER does not require you to use explicit "dummy"
atoms and that's how the A9 tutorial has been set up.

I think the real issue with A9 is that for a reason unknown to me the
original author set up the transformation such that a hydrogen (in
benzene) and the -OH group (in phenol) are duplicated (so you _are_
actually using partial dual topology with both codes).  I do not see why
that would be beneficial in this case and would map the hydrogen to the
oxygen to only have a single disappearing atom.  This would also avoid
the awkward creation of partial charges and simplify the protocol.
When there are only disappearing atoms, setup with Gromacs is very
easy: you only need a single topology file and can vary electrostatics
and vdW lambda separately in a single mdp file.  Results should be the
same.


> Even though one-step transformation seems to give a better result in
> this case, it may be due to luck, because it is suggested that
> one-step transformation is usually not as good as three-step
> transformation, and can even give misleading results in some
> occasions (what occasions?).

The one-step protocol may be a problem with AMBER and relative AFEs.
I have tested this with a set of small organic molecules to compute the
free energy of hydration.  It may be that this depends on the number of
dummy atoms or where the dummies are located or both or...  I have not
found out yet why that is.  The other problem is that AMBER has
difficulties reproducing the end-point geometries e.g. in some cases
bond-length involving dummies are too long.  This is particularly
strange given that AMBER does actually use correct end-points.

In a quick test with Gromacs this did not seem a problem i.e. the
one-step protocol appeared to work satisfactorily.  Maybe also that with
AMBER the situation would be different with binding free energies
(maybe any errors would just happen to cancel out in the thermodynamic
cycle).  But that could be easily checked with the A9 tutorial...


Cheers,
Hannes.


More information about the gromacs.org_gmx-users mailing list