[gmx-users] Performing energy minimisation for an alanine scan
anthony.nash at dpag.ox.ac.uk
Sun Jul 8 22:17:23 CEST 2018
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.
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
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.
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:
> > 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 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
> Gromacs Users mailing list
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> * 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