[gmx-users] Alchemical free energy calculations - right choice of lambda points? What else might I do wrong?
asaffarhi at post.tau.ac.il
asaffarhi at post.tau.ac.il
Thu Jul 30 17:09:39 CEST 2015
You're welcome. Please note that improper dihedral term which relates
between the dummy atoms and the rest of the molecule needs to be
removed ("zeroed") as is explained in Section II. Thanks.
Best regards,
Asaf
Quoting Julian Zachmann <FrankJulian.Zachmann at uab.cat>:
> Thanks a lot for your answers!!! I was really sure that the bonded terms to
> the dummy atoms should be zero so it would have taken me months to find
> this mistake!! I'm now running the simulations with the lambda points like
> this:
>
> ; 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.0500 0.1000 0.1500 0.2000 0.2500 0.3000
> 0.3500 0.4000 0.4667 0.5333 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500
> 0.9000 0.9500 1.0000
> coul_lambdas = 0.0000 0.1000 0.2000 0.3000 0.5000 0.7000 1.0000
> 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
> 1.0000 1.0000 1.0000
> restraint_lambdas = 0.0000 0.1000 0.2000 0.3000 0.5000 0.7000 1.0000
> 1.0000 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.1000 0.2500 0.4000 0.5000 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500
> 0.9000 0.9500 1.0000
> bonded_lambdas = 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
> 0.1000 0.2500 0.4000 0.5000 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500
> 0.9000 0.9500 1.0000
>
> Let's hope it works! Once again, thank you a lot. And of course I will
> quote the articles!
>
> 2015-07-30 10:12 GMT+02:00 <asaffarhi at post.tau.ac.il>:
>
>> Dear Hannes,
>>
>> It means that in many cases dihedral terms do not need to be transformed =
>> not be multiplied by lambda= remain constant. Removing them also reduces
>> the computational efficiency.
>>
>> The alternative is to remove them as is often done ("zero" them).
>>
>> We have shown "not zero" for dihedrals and all the other cases explained
>> in the previous email.
>>
>>
>> Best regards,
>> Asaf
>>
>> Quoting hannes.loeffler at stfc.ac.uk:
>>
>> 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.
>>> --
>>> 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.
>>
>
>
>
> --
> Julian Zachmann
> Laboratori de Medicina Computacional
> Unitat de Bioestadistica. Facultat de Medicina
> Universitat Autonoma de Barcelona
> 08193 Bellaterra (Barcelona). Spain
> Phone: (3493) 581 2797
> Fax: (3493) 581 2344
> E-mail: FrankJulian.Zachmann at uab.cat
> http://lmc.uab.es
> --
> 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