On Mon, 3 Oct 2016 19:57:07 +0000 Guanglin Kuang <guang...@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-requ...@gromacs.org.