[gmx-users] Gromacs Relative Free Energy Calculation Issue

Wed Oct 5 11:28:22 CEST 2016

Dear Hannes,

Part I:
――――
I have tried to use FESetup to generate the topology/coordinate files for the mutated ligands, but I met some problems. The relevant files can be found in this attachment:
https://drive.google.com/open?id=0B8f0-zVoaBXNWWtvYTIwU2ZBRGs

Please go to separate folders with different lambda values, the error information can be found in the .err files. The topology file is lig_FE.top and the coordinate file is ligand_wi.gro, the mdp files for minimization, heating, density equilibration and production run are also included.

I have carefully set up the system and they can be run in the energy minimization steps but all fail in the heating process with similar error messages, complaining that the position of some atoms cannot be set.

I have also tried alchemistry_setup.py, the resultant files would not have this problem and everything runs properly. I think the problem may be due to the fact a direct mutation is used with FESetup, namely mutate H to O. However, with alchemistry_setup.py, it is indirect mutation, namely H to dummy, and dummy to O, separately. This is the only difference I found with FESetup and alchemistry_setup.py. The mdp files are the same.

I have tried a time step at 1 fs and 2 fs, the results are similar for both.

Besides, can the current version of FESetup fit Gromacs 5.X, because I found the generated mdp file use some deprecated parameters, I didn't use them though, I use the latest suggested parameters for Gromacs 5.X.

Part II:
―――――
Even everything is fine with alchemistry_setup.py, I have just reproduced the relative (hydration) free energy difference in solvent. However, I failed to reproduce the relative free energy for the complex, because it seems that the ligand is more dynamic in the protein binding pocket in Gromacs simulation, the mutated ligand would drift around in the simulations. I have repeated and checked many times. I have also checked the simulation results of Amber and found that the ligand is indeed more stable in the binding pocket with Amber and give better results.

Maybe the situation with the complex is due to the fact that benzene and phenol is too weak a binder, or maybe there is something wrong with Gromacs (making the mutated ligand too dynamic), or maybe both reasons. I am trying another system and try to find out the problem.

I have also attach the files for simulation with the complex structure, please have a look at the relevant files. The organization of the files is same as that in the solvent.
https://drive.google.com/open?id=0B8f0-zVoaBXNWWtvYTIwU2ZBRGs

Thank you so much for your help!
Best regards!
Guanglin

-----Original Message-----
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [mailto:gromacs.org_gmx-users-bounces at maillist.sys.kth.se] On Behalf Of Hannes Loeffler
Sent: Tuesday, October 04, 2016 10:51 AM
To: gromacs.org_gmx-users at maillist.sys.kth.se
Cc: gmx-users at gromacs.org
Subject: Re: [gmx-users] Gromacs Relative Free Energy Calculation Issue

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.
--
Gromacs Users mailing list

* Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!

* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists

* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-request at gromacs.org.