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

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


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.

More information about the gromacs.org_gmx-users mailing list