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

[Justin Lemkul]

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

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