[gmx-users] REMD and GBSA

[Ben Reynwar]

Hi gromacs list,

I'm about to start some REMD simulations using generalized Born
solvent on a protein of about 5000 atoms.  I have two questions, the
first of which is about gromacs, the second more about REMD in
general.

(1)
I'm getting some pretty ugly energy drift (300K->500K in 1 ns) for an
NVE MD test simulation using a 2 fs time step.  It goes away if I use
1 fs, however I was under the impression that 2 fs is normally OK.  I
was wondering whether that could be caused by the use of the cut-off
method which is required with the coloumb and VdW interactions when
using GBSA?  Or perhaps I'm doing something else wrong.  I'll include
the mdp file I'm using at the bottom of the email in case anyone feels
like pointing out my foolishness to me.

(2)
I've been analyzing some data from an REMD simulation by my
predecessor and see very slow replica flow rates.  They are about two
orders of magnitude smaller than the idealized rate of the exchange
attempt frequency multiplied by the acceptance fraction (exchanges are
attempted every 2 ps with a 0.4 acceptance fraction).  If I look at
the energy distribution for a given replica/temperature combination
over a time scale of around 1 ns, it is clearly shifted from the
average energy distribution for that temperature.  The timescale for
changes in this energy shift is around 10 ns.  My current theory for
the slow rate of replica flow is that the slow fluctuations in the
energy of the protein are limiting replica flow, since a replica with
lower than average energy will tend to remain at the bottom of the
temperature range, while those with higher than average energies will
tend to remain at the top.  Has anyone else observed this kind of
behavior?  Is my reasoning wrong in any obvious way?

Cheers,
Ben

; An energy drift simulation of TaHSP.

; 7.3.2 Preprocssing
; ------------------
; defines pass to the preprocessor
define =

; 7.3.3 Run Control
; -----------------
; group(s) for center of mass motion removal
comm_grps = System
; Use the MD integrator (as opposed to minimization).
integrator = md
; maximum number of steps to integrate
nsteps = 10000
; remove center of mass translation and rotation around centre of mass
comm_mode = Angular
; [ps] time step for integration
dt = 0.002
; [steps] frequency of mass motion removal
nstcomm = 10
; [ps] starting time for run
tinit = 0

; 7.3.8 Output Control
; --------------------
; [steps] freq to write velocities to trajectory
nstvout = 0
; [steps] freq to write energies to log file
nstlog = 100
; Write to energy file frequently.
nstenergy = 100
; group(s) to write to xtc trajectory
xtc_grps = System
; [real] precision to write xtc trajectory
xtc_precision = 1000
; [steps] freq to write coordinates to xtc trajectory
nstxtcout = 1000
; [steps] freq to write coordinates to trajectory
nstxout = 0
; group(s) to write to energy file
energygrps = System
; [steps] freq to write forces to trajectory
nstfout = 0

; 7.3.9 Neighbour Searching
; -------------------------
; [nm] cut-off distance for the short-range neighbor list
; For Generalized Born this must be equal to the cut-off length for
; the born radius calculation.
rlist = 1.6
; method of updating neighbor list
ns_type = grid
; [steps] freq to update neighbor list
nstlist = 1
; no periodic boundary conditions
pbc = no

; 7.3.10 Electrostatics
; ---------------------
; apply a cut-off to electostatic
coulombtype = cut-off
; coloumb cutoff radius
rcoulomb = 1.6

; 7.3.11 VdW
; ----------
; Dispersion correction makes no sense without box size.
DispCorr = no
; twin-range cut-off with rlist where rvdw >rlist
vdwtype = cut-off
; Increasing VdW cutoff to same as everything else.
rvdw = 1.6

; 7.3.13 Ewald
; ------------
; [nm] grid spacing for FFT grid when using PME
fourierspacing = 0.12
; relative strength of Ewald-shifted potential at rcoulomb
ewald_rtol = 1e-05
; interpolation order for PME, 4 cubic
pme_order = 4

; 7.3.14 Temperature Coupling
; ---------------------------
; Temperature.
ref_t = 300
; temperature coupling frequency
nsttcouple = 1
; Doing NVE simulation so no thermostat.
tcoupl = no
; couple everything to same bath
tc_grps = System
; [ps] time constant for coupling
tau_t = 0.1

; 7.3.15 Pressure Coupling
; ------------------------
; [bar] reference pressure for coupling
ref_p = 1.0
; pressure coupling in x-y-z directions
pcoupltype = isotropic
; [ps] time constant for coupling
tau_p = 2.0
; no pressure coupling if using generalized born
pcoupl = no
; [bar^-1] compressibility
compressibility = 4.5e-05

; 7.3.17 Velocity Generation
; --------------------------
; velocity generation turned off
gen_vel = no

; 7.3.18 Bonds
; ------------
; apply constraints to the start configuration
continuation = yes
; [degrees] maximum angle that a bond can rotate before LINCS will complain
lincs_warnangle = 30
; highest order in the expansion of the contraint coupling matrix
lincs_order = 4
; number of iterations to correct for rotational lengthening
lincs_iter = 1
; LINear Constraint Solver
constraint_algorithm = LINCS
; Bonds to H atoms are constraints
constraints = hbonds

; 7.3.22 NMR refinement (dihedral restraints not in manual)
; ---------------------------------------------------------
; scaling factor for force constant for dihedral constraints
dihre-fc = 10
; No dihedral constraints are set.
dihre = no

; 7.3.27 Implicit Solvent
; -----------------------
; The cutoff radius for calculating the Born radii (must be equal to rlist).
rgbradii = 1.6
; Use a Generalized Born Solvent
implicit_solvent = GBSA
; Use the Onufriev-Bashford-Case method to calculate the Born radii
gb_algorithm = OBC
; Dielectric constant for implicit solvent.
; Default is 80 but we use 78.5 for consistency with amber.
gb_epsilon_solvent = 78.5