[gmx-users] Using the md integrator for calculating free energy of solvation
Dan Gil
dan.gil9973 at gmail.com
Fri May 12 20:52:17 CEST 2017
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?
[ 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.
>
More information about the gromacs.org_gmx-users
mailing list