[gmx-users] Performing energy minimisation for an alanine scan
jalemkul at vt.edu
Sun Jul 8 20:12:08 CEST 2018
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 A. Lemkul, Ph.D.
Virginia Tech Department of Biochemistry
303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061
jalemkul at vt.edu | (540) 231-3129
More information about the gromacs.org_gmx-users