GROMACS version:2021.6
GROMACS modification: No
Here post your question
Hi all,
I’m attempting to run non-eq switching transitions from long equilibrated end states, in order to calculate the free energy difference between lig A and lig B in both complex and solvent. For this Non-EQ switching I’m following the protocol described by Gapsys et al (DOI: 10.1039/C9SC03754C)
In order to increase sampling, I want to run longer equilibration end states. To achieve this with finite computational resources I have decided to increase my timestep from 2fs to 4fs. I therefore repartitioned all hydrogen atoms (excluding solvent) for both protein and ligand, using the default 4x scaling factor. I constrained h-bonds using Lincs algorithm, and following recommendations in this blog post, I also increased lincs-order from 4 to 6 and removed lincs constraints during energy minimisation.
However, when running equilibration after a successful minimisation I get the following warning:
WARNING 2 [file top_complex.top, line 44289]:
The bond in molecule-type MOL between atoms 29 C and 35 DC has an
estimated oscillational period of 1.7e-02 ps, which is less than 5 times
the time step of 4.0e-03 ps.
Maybe you forgot to change the constraints mdp option.
So I then changed constraints from h-bonds to all-bonds, and the simulation appears to be running without breaking. It’s my understanding that this will increase computation cost and I was wondering if there are alternative methods to deal with large timesteps? I have seen a lot of discussion over the choice of integrator used during large time step MD, however it’s not clear to me as to whether there is consensus.
This lengthy OpenMM github thread suggests using a BAOABLangevinIntegrator with longer timesteps and HMR, which is supposedly equivalent to the gromacs SD integrator according to this recent paper - DOI 10.1021/acs.jctc.2c00585. Does anyone have suggestions over choice of integrator? And would switching from my current MD integrator perhaps allow me to remove the constraints=all-bonds flag?
My current MDP settings are listed as follows:
; RUN CONTROL PARAMETERS
integrator = md
; Start time and timestep in ps
tinit = 0
dt = 0.004
nsteps = 5000000
;20 ns
ld_seed = -1
; For exact run continuation or redoing part of a run
init_step = 0
; mode for center of mass motion removal
comm-mode = Linear
; number of steps for center of mass motion removal
nstcomm = 100
nstcalcenergy = 100
nstdhdl = 10000
; group(s) for center of mass motion removal
comm-grps =
; ENERGY MINIMIZATION OPTIONS
; Force tolerance and initial step-size
emtol = 100
emstep = 0.01
; Max number of iterations in relax_shells
niter = 0
; Step size (1/ps^2) for minimization of flexible constraints
fcstep = 0
; Frequency of steepest descents steps when doing CG
nstcgsteep = 1000
nbfgscorr = 10
; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout = 0
nstvout = 0
nstfout = 0
; Checkpointing helps you continue after crashes
nstcheckpoint = 1000
; Output frequency for energies to log file and energy file
nstlog = 10000
nstenergy = 23500
; Output frequency and precision for xtc file
nstxtcout = 23500
xtc-precision = 1000
; This selects the subset of atoms for the xtc file. You can
; select multiple groups. By default all atoms will be written.
xtc-grps =
; Selection of energy groups
energygrps =
; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
cutoff-scheme = verlet
nstlist = 10
; ns algorithm (simple or grid)
ns-type = Grid
; Periodic boundary conditions: xyz (default), no (vacuum)
; or full (infinite systems only)
pbc = xyz
; nblist cut-off
rlist = 1.2
domain-decomposition = no
; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype = PME
rcoulomb-switch = 0
rcoulomb = 1.1
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-r = 1
; Method for doing Van der Waals
vdw-type = switch
; cut-off lengths
rvdw-switch = 1.0
rvdw = 1.1
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Extension of the potential lookup tables beyond the cut-off
table-extension = 1
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
; EWALD/PME/PPPM parameters
pme_order = 4
ewald_rtol = 1e-05
ewald_geometry = 3d
epsilon_surface = 0
optimize_fft = no
; GENERALIZED BORN ELECTROSTATICS
; Algorithm for calculating Born radii
gb_algorithm = Still
; Frequency of calculating the Born radii inside rlist
nstgbradii = 1
; Cutoff for Born radii calculation; the contribution from atoms
; between rlist and rgbradii is updated every nstlist steps
rgbradii = 2
; Salt concentration in M for Generalized Born models
gb_saltconc = 0
; IMPLICIT SOLVENT (for use with Generalized Born electrostatics)
implicit_solvent = No
; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl = v-rescale
; Groups to couple separately
tc-grps = System
; Time constant (ps) and reference temperature (K)
tau-t = 0.1
ref-t = 298
; Pressure coupling
Pcoupl = Parrinello-Rahman
Pcoupltype = Isotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 5
compressibility = 4.6E-5
ref-p = 1
; Random seed for Andersen thermostat
andersen_seed = 815131
; SIMULATED ANNEALING
; Type of annealing for each temperature group (no/single/periodic)
annealing = no
; Number of time points to use for specifying annealing in each group
annealing_npoints = 2
; List of times at the annealing points for each group
annealing_time = 0 250
; Temp. at each annealing point, for each group.
annealing_temp = 0 298
; GENERATE VELOCITIES FOR STARTUP RUN
gen-vel = no
gen-temp = 298
gen-seed = 173529
; OPTIONS FOR BONDS
constraints = all-bonds
; Type of constraint algorithm
constraint-algorithm = Lincs
; Do not constrain the start configuration
unconstrained-start = yes
; Use successive overrelaxation to reduce the number of shake iterations
Shake-SOR = no
; Relative tolerance of shake
shake-tol = 1e-04
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 6
; Number of iterations in the final step of LINCS. 1 is fine for
; normal simulations, but use 2 to conserve energy in NVE runs.
; For energy minimization with constraints it should be 4 to 8.
lincs-iter = 2
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30
; Convert harmonic bonds to morse potentials
morse = no
; NMR refinement stuff
; Distance restraints type: No, Simple or Ensemble
disre = No
; Force weighting of pairs in one distance restraint: Conservative or Equal
disre-weighting = Equal
; Use sqrt of the time averaged times the instantaneous violation
disre-mixed = no
disre-fc = 1000
disre-tau = 0
; Output frequency for pair distances to energy file
nstdisreout = 100
; Orientation restraints: No or Yes
orire = no
; Orientation restraints force constant and tau for time averaging
orire-fc = 0
orire-tau = 0
orire-fitgrp =
; Output frequency for trace(SD) to energy file
nstorireout = 100
; Dihedral angle restraints: No, Simple or Ensemble
dihre = No
dihre-fc = 1000
dihre-tau = 0
; Output frequency for dihedral values to energy file
nstdihreout = 100
; Free energy control stuff
free-energy = yes
init-lambda = 0
delta-lambda = 0
sc-alpha = 0.3
sc-sigma = 0.25
sc-power = 1
sc-coul = yes
;alpha_Q = 0.3
;sigma_Q = 1.0
refcoord-scaling = com
; Non-equilibrium MD stuff
cos-acceleration = 0
Many thanks,
Noah