[gmx-users] Using the md integrator for calculating free energy of solvation
Justin Lemkul
jalemkul at vt.edu
Mon May 15 14:12:27 CEST 2017
On 5/12/17 2:52 PM, Dan Gil wrote:
> I ran equilibrium simulations, and then a longer simulations with the
> recommended sd integrator. Still, the mass_lambda dH/dL is nonzero and
> sometimes it is the dominant contribution at certain lambdas. I don't think
> it is a convergence issue because there aren't huge fluctuations in dH/dL
> values seen in dhdl.xvg.
>
> I also ran simulations where only the mass changed as a function of lambda.
> Others (i.e. Electrostatic, Van der Waals, bonded, and restraints) were
> kept constant. The mass still has an effect on the Helmholtz free energy
> change. Is there possibly a mistake with how I've prepared the topology?
>
Can you provide the section of your .mdp with free energy options? I don't see
it in this chain of emails. It has been suggested previously to leave
mass-lambdas in one of the end states (A or B, doesn't matter) while doing a
relative free energy calculation on the assumption that the contributions
cancel. Perhaps this is not the case, but I have very little else to contribute
on the topic because I am not intimately familiar with those aspects of the
code. Maybe it's a constraint issue? You're going from H to F and presumably
the C-H bonds are to be constrained with the C-F should not be. I doubt this is
it (I have successfully done H->Cl transformations and the results match between
GROMACS and CHARMM), but it's something to check.
-Justin
> [ moleculetype ]
> HEPT 3
>
> [ atoms ]
> 1 opls_135 1 HEPT CH3 1 -0.180 12.011
> opls_961 0.360 12.011
> 2 opls_136 1 HEPT CH2 2 -0.120 12.011
> opls_962 0.240 12.011
> 3 opls_136 1 HEPT CH2 3 -0.120 12.011
> opls_962 0.240 12.011
> ...
> 21 opls_140 1 HEPT H 7 0.060 1.008
> opls_965 -0.120 18.998
> 22 opls_140 1 HEPT H 7 0.060 1.008
> opls_965 -0.120 18.998
> 23 opls_140 1 HEPT H 7 0.060 1.008
> opls_965 -0.120 18.998
>
> ... [ bonds ] ... [ angles ] ... etc...
>
> On Thu, May 11, 2017 at 11:33 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>>
>>
>> On 5/11/17 9:49 AM, Dan Gil wrote:
>>
>>> Thank you Dr. Lemkul,
>>>
>>> The simulations ran, but I have some a question about the results.
>>>
>>> Is it natural to have a nonzero mass_lambda dH/dL? My intuition was that
>>> the potential energy does not depend on mass, and my thinking was that the
>>> mass contribution should be very small, if not zero.
>>>
>>>
>> It should. Maybe it's a convergence issue?
>>
>>
>> -Justin
>>
>> Best Regards,
>>>
>>> Dan
>>>
>>> On Wed, May 3, 2017 at 10:46 AM, Justin Lemkul <jalemkul at vt.edu> wrote:
>>>
>>>
>>>>
>>>> On 5/3/17 10:43 AM, Dan Gil wrote:
>>>>
>>>> Hi Dr. Kausar,
>>>>>
>>>>> If I want the molecule to have all of its nonbonded interactions with
>>>>> the
>>>>> rest of the system (solvent) at each lambda, don't I want to have vdwq
>>>>> option ON for both couple-lambda0 and couple-lambda1?
>>>>>
>>>>> What I am attempting to do is mutate molecule A to molecule B through
>>>>> the
>>>>> simulation.
>>>>>
>>>>>
>>>>> Relative free energy calculations are a little bit of a strange beast;
>>>> you
>>>> do need to set the B-state (couple-lambda1) to "none" and the B-state
>>>> types
>>>> and charges are read from [atoms].
>>>>
>>>>
>>>> -Justin
>>>>
>>>> Best Regards,
>>>>
>>>>>
>>>>> Dan
>>>>>
>>>>> On Wed, May 3, 2017 at 1:04 AM, Tasneem Kausar <
>>>>> tasneemkausar12 at gmail.com
>>>>>
>>>>>>
>>>>>> wrote:
>>>>>
>>>>> In the last line of mdp file couple-lambda0 and couple-lambda1 have the
>>>>>
>>>>>> same vdwq variable. It indicates A and B state of system. So grompp is
>>>>>> complaining about the identical states.
>>>>>>
>>>>>> You can also see the alchemistry.org for detail information about the
>>>>>> free
>>>>>> energy calculations.
>>>>>>
>>>>>> On Wed, May 3, 2017 at 2:14 AM, Dan Gil <dan.gil9973 at gmail.com> wrote:
>>>>>>
>>>>>> Thank you Dr. Lemkul,
>>>>>>
>>>>>>>
>>>>>>> Following your advice, I was able to calculate the solvation free
>>>>>>> energy
>>>>>>>
>>>>>>> of
>>>>>>
>>>>>> two molecules that I am interested in. The results, unfortunately, were
>>>>>>>
>>>>>>> not
>>>>>>
>>>>>> what I expected. I want to try an alternative method, that is, the
>>>>>>> alchemical thermodynamic integration method.
>>>>>>>
>>>>>>> If I understand this correctly, I should be able to calculate the free
>>>>>>> energy difference between two systems A and B as I change the topology
>>>>>>> of
>>>>>>> the molecule from A to B in small steps.
>>>>>>>
>>>>>>> I wasn't able to find a tutorial online, but I attempted to try the
>>>>>>>
>>>>>>> method
>>>>>>
>>>>>> and I obtained this error from grompp:
>>>>>>>
>>>>>>> WARNING 1 [file grompp.mdp, line 43]:
>>>>>>> The lambda=0 and lambda=1 states for coupling are identical
>>>>>>>
>>>>>>> I was hoping if you could help me understand what I am missing in my
>>>>>>> topology or mdp file.
>>>>>>>
>>>>>>> Topology:
>>>>>>> [ moleculetype ]
>>>>>>> HEPT 3
>>>>>>>
>>>>>>> [ atoms ]
>>>>>>> 1 opls_135 1 HEPT CH3 1 -0.180
>>>>>>> 12.011
>>>>>>> opls_961 0.360 12.011
>>>>>>> 2 opls_136 1 HEPT CH2 2 -0.120
>>>>>>> 12.011
>>>>>>> opls_962 0.240 12.011
>>>>>>> 3 opls_136 1 HEPT CH2 3 -0.120
>>>>>>> 12.011
>>>>>>> opls_962 0.240 12.011
>>>>>>> 4 opls_136 1 HEPT CH2 4 -0.120
>>>>>>> 12.011
>>>>>>> opls_962 0.240 12.011
>>>>>>> 5 opls_136 1 HEPT CH2 5 -0.120
>>>>>>> 12.011
>>>>>>> opls_962 0.240 12.011
>>>>>>> 6 opls_136 1 HEPT CH2 6 -0.120
>>>>>>> 12.011
>>>>>>> opls_962 0.240 12.011
>>>>>>> 7 opls_135 1 HEPT CH3 7 -0.180
>>>>>>> 12.011
>>>>>>> opls_961 0.360 12.011
>>>>>>> 8 opls_140 1 HEPT H 1 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 9 opls_140 1 HEPT H 1 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 10 opls_140 1 HEPT H 1 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 11 opls_140 1 HEPT H 2 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 12 opls_140 1 HEPT H 2 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 13 opls_140 1 HEPT H 3 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 14 opls_140 1 HEPT H 3 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 15 opls_140 1 HEPT H 4 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 16 opls_140 1 HEPT H 4 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 17 opls_140 1 HEPT H 5 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 18 opls_140 1 HEPT H 5 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 19 opls_140 1 HEPT H 6 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 20 opls_140 1 HEPT H 6 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 21 opls_140 1 HEPT H 7 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 22 opls_140 1 HEPT H 7 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> 23 opls_140 1 HEPT H 7 0.060
>>>>>>> 1.008
>>>>>>> opls_965 -0.120 18.998
>>>>>>> ... And so on with the bonds, pairs, angles, and dihedrals directive.
>>>>>>>
>>>>>>> MDP File:
>>>>>>>
>>>>>>> ;Integration Method and Parameters
>>>>>>> integrator = md
>>>>>>> nsteps = 10000
>>>>>>> dt = 0.002
>>>>>>> nstenergy = 1000
>>>>>>> nstlog = 5000
>>>>>>>
>>>>>>> ;Output Control
>>>>>>> nstxout = 0
>>>>>>> nstvout = 0
>>>>>>>
>>>>>>> ;Cutoff Schemes
>>>>>>> cutoff-scheme = group
>>>>>>> rlist = 1.0
>>>>>>> vdw-type = cut-off
>>>>>>> rvdw = 2.0
>>>>>>>
>>>>>>> ;Coulomb interactions
>>>>>>> coulombtype = pme
>>>>>>> rcoulomb = 1.0
>>>>>>> fourierspacing = 0.4
>>>>>>>
>>>>>>> ;Constraints
>>>>>>> constraints = all-bonds
>>>>>>>
>>>>>>> ;Temperature coupling
>>>>>>> tcoupl = v-rescale
>>>>>>> tc-grps = system
>>>>>>> tau-t = 0.1
>>>>>>> ref-t = 300
>>>>>>>
>>>>>>> ;Pressure coupling
>>>>>>> pcoupl = parrinello-rahman
>>>>>>> ref-p = 1.01325
>>>>>>> compressibility = 4.5e-5
>>>>>>> tau-p = 5
>>>>>>>
>>>>>>> ;Free energy calculation
>>>>>>> free-energy = yes
>>>>>>> init-lambda = 0
>>>>>>> delta-lambda = 0.00001
>>>>>>> couple-moltype = HEPT
>>>>>>> couple-lambda0 = vdwq
>>>>>>> couple-lambda1 = vdwq
>>>>>>>
>>>>>>> On Mon, Apr 3, 2017 at 2:05 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>>> wrote:
>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>> On 4/3/17 2:01 PM, Dan Gil wrote:
>>>>>>>>
>>>>>>>> Thank you Dr. Lemkul,
>>>>>>>>
>>>>>>>>>
>>>>>>>>> I am trying to run grompp with the md integrator, but I am getting
>>>>>>>>>
>>>>>>>>> this
>>>>>>>>
>>>>>>>
>>>>>> error:
>>>>>>>
>>>>>>>> "For proper sampling of the (nearly) decoupled state, stochastic
>>>>>>>>>
>>>>>>>>> dynamics
>>>>>>>>
>>>>>>>
>>>>>>> should be used"
>>>>>>>>
>>>>>>>>>
>>>>>>>>> Should I ignore this warning with the -maxwarn option and try
>>>>>>>>> running
>>>>>>>>>
>>>>>>>>> it?
>>>>>>>>
>>>>>>>
>>>>>>> I
>>>>>>>>
>>>>>>>>> will see if I obtain comparable values for ethanol.
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> People have already done such a comparison and that's why grompp is
>>>>>>>>>
>>>>>>>> telling you this - it's better to use the Langevin integrator.
>>>>>>>> You'll
>>>>>>>>
>>>>>>>> get
>>>>>>>
>>>>>>> better sampling, particularly towards the end states. Using -maxwarn
>>>>>>>>
>>>>>>>> tells
>>>>>>>
>>>>>>> grompp "you're trying to prevent me from making a mistake, but I want
>>>>>>>>
>>>>>>>> to
>>>>>>>
>>>>>>
>>>>>> do
>>>>>>>
>>>>>>> it anyway" :) You'd better have a really, really good reason to try
>>>>>>>> to
>>>>>>>> override it.
>>>>>>>>
>>>>>>>> -Justin
>>>>>>>>
>>>>>>>>
>>>>>>>> Best Regards,
>>>>>>>>
>>>>>>>>
>>>>>>>>> Dan
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> On Mon, Mar 27, 2017 at 12:51 PM, Justin Lemkul <jalemkul at vt.edu>
>>>>>>>>>
>>>>>>>>> wrote:
>>>>>>>>
>>>>>>>
>>>>>>>
>>>>>>>>
>>>>>>>>>
>>>>>>>>> On 3/26/17 9:40 PM, Dan Gil wrote:
>>>>>>>>>>
>>>>>>>>>> Hi,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>> I am following Dr. Sander Pronk's and Dr. Justin Lemkul's tutorial
>>>>>>>>>>>
>>>>>>>>>>> on
>>>>>>>>>>
>>>>>>>>>
>>>>>> calculating free energy of solvation. Is it possible and
>>>>>>>
>>>>>>>>
>>>>>>>>>>> theoretically
>>>>>>>>>>
>>>>>>>>>
>>>>>> sound to use the md integrator instead of the sd integrator for
>>>>>>>
>>>>>>>>
>>>>>>>>>>> these
>>>>>>>>>>
>>>>>>>>>
>>>>>> calculations?
>>>>>>>
>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> Langevin dynamics gives better sampling so it is frequently used
>>>>>>>>>>> for
>>>>>>>>>>>
>>>>>>>>>>> free
>>>>>>>>>> energy calculations. You may get comparable results with the
>>>>>>>>>>
>>>>>>>>>> leap-frog
>>>>>>>>>
>>>>>>>>
>>>>>> integrator, but I haven never done a side-by-side comparison.
>>>>>>>
>>>>>>>>
>>>>>>>>>> -Justin
>>>>>>>>>>
>>>>>>>>>> I have already done a considerable amount of work using md
>>>>>>>>>>
>>>>>>>>>> integration,
>>>>>>>>>
>>>>>>>>
>>>>>> and
>>>>>>>
>>>>>>>>
>>>>>>>>>> I want to make sure that the free energy values I calculate are
>>>>>>>>>>
>>>>>>>>>>> consistent
>>>>>>>>>>> with my previous work.
>>>>>>>>>>>
>>>>>>>>>>> If using the md integrator is not sound, is there an alternative
>>>>>>>>>>> way
>>>>>>>>>>>
>>>>>>>>>>> of
>>>>>>>>>>
>>>>>>>>>
>>>>>>> calculating solvation energy that will be consistent?
>>>>>>>>
>>>>>>>>>
>>>>>>>>>>> Best Regards,
>>>>>>>>>>>
>>>>>>>>>>> Dan
>>>>>>>>>>>
>>>>>>>>>>>
>>>>>>>>>>> --
>>>>>>>>>>>
>>>>>>>>>>> ==================================================
>>>>>>>>>>
>>>>>>>>>> Justin A. Lemkul, Ph.D.
>>>>>>>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>>>>>>>
>>>>>>>>>> Department of Pharmaceutical Sciences
>>>>>>>>>> School of Pharmacy
>>>>>>>>>> Health Sciences Facility II, Room 629
>>>>>>>>>> University of Maryland, Baltimore
>>>>>>>>>> 20 Penn St.
>>>>>>>>>> Baltimore, MD 21201
>>>>>>>>>>
>>>>>>>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>>>>>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>>>>>>>
>>>>>>>>>> ==================================================
>>>>>>>>>> --
>>>>>>>>>> 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.
>>>>>>>
>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --
>>>>>>>>>>
>>>>>>>>> ==================================================
>>>>>>>>
>>>>>>>> Justin A. Lemkul, Ph.D.
>>>>>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>>>>>
>>>>>>>> Department of Pharmaceutical Sciences
>>>>>>>> School of Pharmacy
>>>>>>>> Health Sciences Facility II, Room 629
>>>>>>>> University of Maryland, Baltimore
>>>>>>>> 20 Penn St.
>>>>>>>> Baltimore, MD 21201
>>>>>>>>
>>>>>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>>>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>>>>>
>>>>>>>> ==================================================
>>>>>>>> --
>>>>>>>> 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.
>>>>>>
>>>>>>
>>>>>> --
>>>> ==================================================
>>>>
>>>> Justin A. Lemkul, Ph.D.
>>>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>>>
>>>> Department of Pharmaceutical Sciences
>>>> School of Pharmacy
>>>> Health Sciences Facility II, Room 629
>>>> University of Maryland, Baltimore
>>>> 20 Penn St.
>>>> Baltimore, MD 21201
>>>>
>>>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>>>> http://mackerell.umaryland.edu/~jalemkul
>>>>
>>>> ==================================================
>>>> --
>>>> 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.
>>>>
>>>>
>> --
>> ==================================================
>>
>> Justin A. Lemkul, Ph.D.
>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
>>
>> Department of Pharmaceutical Sciences
>> School of Pharmacy
>> Health Sciences Facility II, Room 629
>> University of Maryland, Baltimore
>> 20 Penn St.
>> Baltimore, MD 21201
>>
>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
>> http://mackerell.umaryland.edu/~jalemkul
>>
>> ==================================================
>> --
>> 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.
>>
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list