[gmx-users] Performing energy minimisation for an alanine scan

Justin Lemkul jalemkul at vt.edu
Mon Jul 9 00:34:43 CEST 2018



On 7/8/18 6:00 PM, Anthony Nash wrote:
> Just looked through the log file of one of the energy minimisations and I have a few questions.
>
>
> Recap, my attempt is to immediately represent each structure as stateB i.e., the alanine substitution. I've used the .mdp values below, switching all lambdas into 1.0, by only having a single value to use and putting the initial lambda state as 0; the first column entry and the last given that there is only one lambda value to use.
>
>
> In the log file I am seeing energies:
>
>
>     Energies (kJ/mol)
>
>             Bond          Angle    Proper Dih.  Improper Dih.          LJ-14
>
>      1.29731e+03    6.85669e+03    2.14601e+04    7.28711e+02    7.88880e+03
>
>       Coulomb-14        LJ (SR)  Disper. corr.   Coulomb (SR)   Coul. recip.
>
>      9.15559e+04    6.07089e+05   -1.12079e+04   -4.54144e+06    7.23600e+03
>
>        Potential Pres. DC (bar) Pressure (bar)      dVcoul/dl       dVvdw/dl
>
>     -3.80853e+06    0.00000e+00   -2.86070e+04   -5.72155e+01   -1.65947e+03
>
>      dVbonded/dl dVrestraint/dl
>
>      0.00000e+00    0.00000e+00
>
>
> (not sure whether copy and paste worked well via this mailing list)
>
>
> I am guessing I am seeing coul and vdw energies with respect to their lambda value, my concern is the lack of bonded energy with respect to lambda. This was going from a big amino acid into alanine, so I thought it might just not be using certain bonds. So I then tried the other way round, glycine going into alanine. This dVbonded/dl are still zero.
>
>
> I might be barking up the wrong tree, is this a concern?

Are you still getting this from grompp:

"Some parameters for bonded interaction involving perturbed atoms are
specified explicitly in state A, but not B - copying A to"

If so, it's an indication that you don't have the right bonded parameters (i.e., they're not explicitly supplied, therefore you're getting A-state bonded parameters and B-state nonbonded).

-Justin



>
> Thanks
>
> Anthony
>
>
> Kind regards
> Anthony Nash PhD MRSC
> Department of Physiology, Anatomy, and Genetics
> University of Oxford
> ________________________________
> From: Anthony Nash
> Sent: 08 July 2018 22:45:39
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Performing energy minimisation for an alanine scan
>
>
> Thanks Mark, sounds like a good idea.
>
>
> Kind regards
> Anthony Nash PhD MRSC
> Department of Physiology, Anatomy, and Genetics
> University of Oxford
> ________________________________
> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Mark Abraham <mark.j.abraham at gmail.com>
> Sent: 08 July 2018 22:38:41
> To: gmx-users at gromacs.org
> Subject: Re: [gmx-users] Performing energy minimisation for an alanine scan
>
> Hi,
>
> The log and energy output reports the lambda value. And of course you can
> compare the value from a single point energy from a single-topology
> equivalent coordinate input.
>
> Mark
>
> On Sun, Jul 8, 2018 at 10:17 PM Anthony Nash <anthony.nash at dpag.ox.ac.uk>
> wrote:
>
>> Hi Justin and Mark,
>>
>>
>> Thank you for your replies. Very helpful indeed.
>>
>> In answer to Justin's question, yes, I simply used the Free Energy .mdp
>> options as a hack to get the structure into state B. This was a shortcut to
>> the horror which would have been opening 600+ structures in something like
>> Maestro and make the change by hand.
>>
>>
>> As a point of validation, besides blind faith in Gromacs, what is there
>> that I can do to ensure stateB was used? I have noticed a warning:
>>
>>
>>    "Some parameters for bonded interaction involving perturbed atoms are
>> specified explicitly in state A, but not B - copying A to"
>>
>> This was from aspartic acid(stateA) -> Alanine(stateB), after using .mdp
>> options which should have immediately put the structure into stateB.
>>
>> Mark, that's a good idea. I think it would be an excellent feature of pmx,
>> I'll drop them an email. I am not too concerned by the overall performance.
>> I am only using the lambda endpoint 1. The 600+ simulations that I will be
>> performing is to map a drug binding site volume (could also be overall
>> protein SASA if the entry site closes) as a function of a single point
>> mutation in the protein sequence. I won't be pushing each structure for
>> very long after equilibrium.
>>
>>
>> Thanks both
>>
>> Anthony
>>
>>
>> Kind regards
>> Anthony Nash PhD MRSC
>> Department of Physiology, Anatomy, and Genetics
>> University of Oxford
>> ________________________________
>> From: gromacs.org_gmx-users-bounces at maillist.sys.kth.se <
>> gromacs.org_gmx-users-bounces at maillist.sys.kth.se> on behalf of Mark
>> Abraham <mark.j.abraham at gmail.com>
>> Sent: 08 July 2018 19:21:17
>> To: gmx-users at gromacs.org
>> Subject: Re: [gmx-users] Performing energy minimisation for an alanine scan
>>
>> Hi,
>>
>> Justin is definitely right about performance for lambda not 0 or 1, but
>> pretty much all of the forms of slowness for FE calculations are trivially
>> avoidable for the endpoints. I just don't know that the code does that, and
>> I can't check the dozens of places to be sure, so I encourage Anthony to
>> compare performance, if he wants to proceed with this approach. As Justin
>> suggests, it would certainly be simpler if there was a way to generate the
>> coordinates of the mutated state as a single-topology structure. Do get in
>> touch with the authors of pmx, though, as if it's not currently possible it
>> sounds to me like a useful feature.
>>
>> Mark
>>
>> On Sun, Jul 8, 2018 at 8:12 PM Justin Lemkul <jalemkul at vt.edu> wrote:
>>
>>>
>>> On 7/7/18 6:09 PM, Anthony Nash wrote:
>>>> I'm trying to minimise 600+ structure. They differ by an alanine
>>> substitution, which has been put in place using stateA/stateB topologies.
>>> The technique used to generate the initial structure was taken from:
>>>
>> http://pmx.mpibpc.mpg.de/output.pl?jobid=model_scanA_amber99sb-star-ildn-mut.ff_3329_20180703151957
>>>>
>>>> This produces a topology of:
>>>>
>>>> stateA -> wildtype amino acid
>>>>
>>>> stateB -> alanine substitution
>>>>
>>>>
>>>> Given the number of structures, I am trying to automate everything.
>>> Centering, solvating, genion etc., all done. Now I come to an initial
>>> energy minimisation for each structure. Unfortunately, emtol is aren't
>>> falling < 1000, which I find to be a value sufficient to move on from
>>> energy minimisation to more dynamical equilibrium methodologies. One
>> could
>>> argue that I should look at the trajectory, but I want the parameters
>>> correct first, especially given how I am automating this.
>>>>
>>>> Since the alanine scan methodology generated stateA/stateB topologies
>>> per structure, I figured that I need to use free energy parameters in my
>>> .mdp file to ensure that I am immediately only ever using the stateB
>>> (alanine) topology.
>>>>
>>>> I basically need someone to sanity check the .mdp values to ensure that
>>> I only use stateB (the alanine topology) during the energy minimisation.
>>> The complete free energy .mdp values (from the mdout.mdp) file are:
>>>>
>>>> ; Free energy variables
>>>> free-energy              = yes
>>>> couple-moltype           =
>>>> couple-lambda0           = vdw-q
>>>> couple-lambda1           = vdw-q
>>>> couple-intramol          = no
>>>> init-lambda              = -1
>>>> init-lambda-state        = 0
>>>> delta-lambda             = 0
>>>> nstdhdl                  = 10
>>>> fep-lambdas              =
>>>> mass-lambdas             =
>>>> coul-lambdas             = 1.0
>>>> vdw-lambdas              = 1.0
>>>> bonded-lambdas           = 1.0
>>>> restraint-lambdas        = 1.0
>>>> temperature-lambdas      =
>>>> calc-lambda-neighbors    = 1
>>>> init-lambda-weights      =
>>>> dhdl-print-energy        = no
>>>> sc-alpha                 = 0
>>>> sc-power                 = 1
>>>> sc-r-power               = 6
>>>> sc-sigma                 = 0.3
>>>> sc-coul                  = no
>>>> separate-dhdl-file       = yes
>>>> dhdl-derivatives         = yes
>>>> dh_hist_size             = 0
>>>> dh_hist_spacing          = 0.1
>>>>
>>>> Please note, key-value pairs like coul-lambdas and vdw-lambdas are 1.0
>>> as I only want to model state B and I am not interested in slow growth.
>>>>
>>>> Are these parameters sufficient for modelling stateB upon immediate
>>> execution of an energy minimisation (and NPT/NVT dynamics)? I've had very
>>> little experience with the lambda implementation in gromacs and I turn to
>>> more experience.
>>>
>>> They should, yes, but are you actually computing free energy changes? If
>>> this is just a hack to use a single starting structure and model Ala in
>>> each position, that can certainly work, but the simulations will be
>>> slower (because there is overhead associated with the free energy code)
>>> and you may face convergence issues like you're mentioning. If you don't
>>> need to compute DDG, then why not just mutate each residue and have
>>> normal PDB files to work with? Presumably the automation would be just
>>> as easy, if not easier.
>>>
>>> -Justin
>>>
>>> --
>>> ==================================================
>>>
>>> Justin A. Lemkul, Ph.D.
>>> Assistant Professor
>>> Virginia Tech Department of Biochemistry
>>>
>>> 303 Engel Hall
>>> 340 West Campus Dr.
>>> Blacksburg, VA 24061
>>>
>>> jalemkul at vt.edu | (540) 231-3129
>>> http://www.thelemkullab.com
>>>
>>> ==================================================
>>>
>>> --
>>> 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.

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================



More information about the gromacs.org_gmx-users mailing list