[gmx-users] REMD and GBSA

[Ben Reynwar]

On Fri, Oct 14, 2011 at 9:11 AM, Mark Abraham <Mark.Abraham at anu.edu.au> wrote:
> 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.
>
Great. Thank you.

>>
>> (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
>

Do you know if anyone has done any studies looking at replica flow in
well-defined, comparatively low-dimensional landscapes to get a
qualitative feel for these kinds of effects?  It would be interesting
to see what effect replica exchange settings can have on replica flow
beyond the simple random walk models, when you take into account the
fact that different regions of configuration space could have
different potential energy distributions for the same temperature.

Cheers,
Ben

>>
>> 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 at 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>