[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 10:12:38 CEST 2015


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.





More information about the gromacs.org_gmx-users mailing list