[gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?
hannes.loeffler at stfc.ac.uk
hannes.loeffler at stfc.ac.uk
Thu Jul 30 09:57:51 CEST 2015
The argument regards dihedral that was made towards me was about convergence.
When you say "_can_ remain constant", what would be the alternative? I would still think "not zero" for bonds and angles. That is really the crucial point here.
________________________________________
From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of asaffarhi at post.tau.ac.il [asaffarhi at post.tau.ac.il]
Sent: 30 July 2015 08:45
To: gromacs.org_gmx-users at maillist.sys.kth.se
Subject: Re: [gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?
Hi,
It has been proven both analytically (http://arxiv.org/abs/1310.2112,
Section II) and in simulations (http://arxiv.org/abs/1310.2112,
Section IV and http://arxiv.org/abs/1501.01514) that the dihedral
terms can remain constant in the transformation in relative free
energy calculations. In the second of these articles the dihedral term
has been removed in a transformation both in vacuum and water and gave
the same free energy difference, which means that removing them is an
unnecessary operation (the contributions cancel out in the
Thermodynamic Cycle). If there is an energy barrier in the dihedral
potential significantly higher than k_B*T and there are several
valleys it may or may not need to be multiplied by a value between 0
and 1 to help with convergence (see http://arxiv.org/abs/1310.2112,
Section VI). In addition any bonded potential function which relates
between the decoupled (dummy) atoms and the rest of the system -
harmonic or non-harmonic (quadratic or non-quadratic) can remain
constant. Also, bond angle potentials of sub-molecules with coupled
degrees of freedom such as the methyl group can remain constant even
if they relate between the decoupled atoms and the rest of the
molecule (http://arxiv.org/abs/1310.2112, Section II). Please cite
these articles if you are using them.
Thanks,
Best regards,
Asaf
Quoting hannes.loeffler at stfc.ac.uk:
> Hi,
>
> you wouid _not_ set bonded terms to zero. That would mean that you
> have essentially a non-bonded end state allowing the atom to be
> "everywhere" (reading Stefan Boresch' papers from 1999 on this may
> be helpful). The dihedrals are debatable (a collaborator of mine
> thinks that having them zero may help in replicate methods or so).
> Our strategy with FESetup (http://www.hecbiosim.ac.uk/fesetup) is to
> copy those terms over from the state without the dummies.
>
> The mass-lambdas are best kept at one end state as discussed
> previously to avoid any bad interactions with constraints (you say
> you have X-C-H3 to X-H-DU3).
>
> Cheers,
> Hannes.
>
> ________________________________________
> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se
> [gromacs.org_gmx-users-bounces at maillist.sys.kth.se] on behalf of
> Julian Zachmann [FrankJulian.Zachmann at uab.cat]
> Sent: 29 July 2015 14:10
> To: gmx-users at gromacs.org
> Subject: [gmx-users] Alchemical free energy calculations - right
> choice of lambda points? What else might I do wrong?
>
> Dear Gromacs-Users,
>
> I calculate relative binding energies using free energy perturbation (FEP)
> and mbar. I have the binding data for 4 receptors and two ligands which are
> almost identical. The only difference is the change from a methyl group to
> a hydrogen. I alchemically convert a methyl group into one hydrogen atom
> with 3 dummy atoms - inside the receptor and in water. First I have several
> equilibration steps; the production run takes 4ns. Unfortunately my results
> don't match the experimental results and I would like to ask for input
> about what I might do wrong.
>
> First the results:
>
> Experimental | Theoretical (kcal) mbar calculated with pymbar
> Receptor 1: 0.85 | 0.57
> Receptor 2: 0.43 | -1.94
> Receptor 3: 1.13 | 0.78
> Receptor 4: 2.41 | 0.29
>
> I am using the following lambda points:
> ; 0 1 2 3
> 4 5 6 7 8 9 10
> 11 12 13 14 15 16 17
> 18 19
> fep_lambdas = 0.0000 0.2000 0.4000 0.5000 0.6000 0.7000 0.8000
> 0.9000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
> 1.0000 1.0000 1.0000
> vdw_lambdas = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
> 0.0000 0.0000 0.2000 0.4000 0.6000 0.7000 0.8000 0.9000 0.9300 0.9600
> 0.9900 0.9950 1.0000
>
> I think this might be a mistake as the graphic of free energy differences
> <https://drive.google.com/open?id=0B2M9aqeJrxnYMWRXbGZ4dkxrb00> shows there
> is a huge difference between lambda points 7-8 - when fep_lambdas turn from
> 0.9 to 1.0. I think of either introducing another lambda point such as
> fep_lambdas = 0.95 or to split up fep_lambdas and converge bonded_lambdas,
> coul_lambdas, mass_lambdas differently fast to 1.0. I think
> temperature_lambdas and restraint_lambdas are not important for my system
> so I my gradually convert them from 0.0 to 1.0 over all 20 lambda points. I
> would be very thankful if you could advise me how it would be best to
> converge the bonded-, coul-, mass_lambdas, and vdw_lambdas to 1.0.
>
> Other guesses about what I might do wrong are in the mol.itp
> <https://drive.google.com/file/d/0B2M9aqeJrxnYeGFBMl9XTmprZXM/view?usp=sharing>
> file:
>
> I am not 100% sure, but I think it is correct to put the bond, angle, and
> dihedral strengths for dummy atoms to zero.
>
> [ bonds ]
> ; ai aj funct r k
> (...)
> 16 18 1 1.0920e-01 2.8225e+05 1.0920e-01 2.8225e+05
> 13 14 1 1.0910e-01 2.8342e+05 1.0910e-01 2.8342e+05
> 13 15 1 1.0910e-01 2.8342e+05 1.0910e-01 2.8342e+05
> 8 9 1 1.0910e-01 2.8342e+05 1.0910e-01 0.0000
> 8 10 1 1.0910e-01 2.8342e+05 1.0910e-01 0.0000
> 8 11 1 1.0910e-01 2.8342e+05 1.0910e-01 0.0000
> 7 12 1 1.0330e-01 3.0878e+05 1.0330e-01 3.0878e+05
> 4 5 1 1.0910e-01 2.8342e+05 1.0910e-01 2.8342e+05
> 4 6 1 1.0910e-01 2.8342e+05 1.0910e-01 2.8342e+05
> (...)
>
> [ angles ]
> ; ai aj ak funct theta cth
> (...)
> 13 16 18 1 1.1005e+02 3.8802e+02 1.1005e+02 3.8802e+02
> 12 7 13 1 1.1011e+02 3.8652e+02 1.1011e+02 3.8652e+02
> 10 8 11 1 1.1074e+02 3.2669e+02 1.1074e+02 0.0000
> 9 8 10 1 1.1074e+02 3.2669e+02 1.1074e+02 0.0000
> 9 8 11 1 1.1074e+02 3.2669e+02 1.1074e+02 0.0000
> 8 7 12 1 1.1011e+02 3.8652e+02 1.0811e+02 3.3907e+02
> 7 8 9 1 1.0791e+02 4.1020e+02 1.0791e+02 0.0000
> 7 8 10 1 1.0791e+02 4.1020e+02 1.0791e+02 0.0000
> 7 8 11 1 1.0791e+02 4.1020e+02 1.0791e+02 0.0000
> 7 13 14 1 1.0791e+02 4.1020e+02 1.0791e+02 4.1020e+02
> 7 13 15 1 1.0791e+02 4.1020e+02 1.0791e+02 4.1020e+02
> (...)
>
> [ dihedrals ]
> ;i j k l func C0 ... C5
> (...)
> 12 7 13 15 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000
> 12 7 13 16 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000
> 11 8 7 12 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 11 8 7 13 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 10 8 7 12 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 10 8 7 13 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 9 8 7 12 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 9 8 7 13 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000
> 0.00000 0.00000
> 8 7 13 14 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000
> 8 7 13 15 3 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000 0.65270 1.95811 0.00000 -2.61082
> 0.00000 0.00000
> (...)
>
>
> My feeling tells me, that mdp files
> <https://drive.google.com/open?id=0B2M9aqeJrxnYdDdYYUJoTmk3MDA> are
> correct. I sticked close to the tutorial
> <http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/01_theory.html>
> on
> the Gromacs page.
>
> I would be very thankful for any input you could give me!
>
> Best wishes,
> Julian
> --
> 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.
> --
> 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.
--
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.
More information about the gromacs.org_gmx-users
mailing list