[gmx-users] Performing energy minimisation for an alanine scan
Anthony Nash
anthony.nash at dpag.ox.ac.uk
Mon Jul 9 00:00:21 CEST 2018
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?
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.
More information about the gromacs.org_gmx-users
mailing list