[gmx-users] REMD and GBSA

[Mark Abraham]

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