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

[Anthony Nash]

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.


Many thanks

Anthony


Kind regards
Anthony Nash PhD MRSC
Department of Physiology, Anatomy, and Genetics
University of Oxford