[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