[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
Sun Aug 2 11:34:51 CEST 2015
p.s
The discussion about convergence is relevant only if Replica Exchange
(H-REMD/H-PT) is used and there is a high energy barrier in the
dihedral potential. If the energy barrier is not higher than
~4.5*k_B*T you are welcome to try a calculation with and without
"zeroing" the dihedral term and see that you get the same result.
Thanks,
Best,
Asaf
Quoting asaffarhi at post.tau.ac.il:
> 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.
>
>
>
> --
> 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