On 14/10/2011 10:12 AM, Ben Reynwar wrote:
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.
Drift is normal if you use 2fs with hbond constraints. all-bond
constraints are necessary for 2fs.
(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?
A replica with "lower than average energy" *for that temperature* will
tend to drift down the temperature ladder in favour of another.
One can observe blockages in replica flow. If all the replicas below a
given temperature are in regions of configuration space that cannot
access PE high enough to have a significant chance of exchanging above
that temperature, then flow does not occur (and vice-versa, of course).
If one were to sample a FES that had two minima that should be sampled
in a 2:1 ratio, but started from a 1:1 ratio and did not have a high
enough temperature range to cross the barrier, then the exchange
acceptance rate can look good when nothing useful is occurring - the
observation will necessarily be that these minima are equally likely.
The two groups are actually engaging in disjoint flow, and one needs to
look at metrics other than the acceptance rate to observe it. The only
way to deal with such a bottleneck is to have replicas at a high enough
temperature that both groups can exchange to those temperatures - only
now can barrier crossing occur. These kinds of phenomena can certainly
occur over short time scales in localized regions of configuration and
temperature space.
Mark
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
--
gmx-users mailing list gmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/Search before posting!
Please don't post (un)subscribe requests to the list. Use the
www interface or send it to gmx-users-requ...@gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
