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

Mark Abraham mark.j.abraham at gmail.com
Sun Jul 8 23:38:57 CEST 2018 

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.
>

More information about the gromacs.org_gmx-users mailing list