[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