[gmx-users] REMD - replicas sampling in temperatures beyond the assigned range
Mark Abraham
mark.j.abraham at gmail.com
Thu Jun 30 16:42:37 CEST 2016
Hi,
Best practice is to read and learn others practice from publications that
are similar to what you want to do, rather than making ad-hoc changes. In
this case, the GROMACS defaults are pretty close to the de facto standard,
and supported by analysis work done by other members of the community.
Mark
On Thu, Jun 30, 2016 at 4:16 PM NISHA Prakash <nishnith20591 at gmail.com>
wrote:
> Dear Justin,
>
> Thanks a lot for pointing out the issues. I now understand why there were
> such high oscillations.
>
> Could you please also tell me if there are any ideal values for pme_order
> and fourier spacing with respect to the cut offs' value of 1.4.
>
> Does the following Note imply I can raise the fourier grid spacing to 0.25?
>
> NOTE 2 [file sim-new.mdp]:
> The optimal PME mesh load for parallel simulations is below 0.5
> and for highly parallel simulations between 0.25 and 0.33,
> for higher performance, increase the cut-off and the PME grid spacing
>
> Thank you again,
>
> Nisha
>
>
> On Thu, Jun 30, 2016 at 6:55 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
> >
> >
> > On 6/30/16 9:16 AM, NISHA Prakash wrote:
> >
> >> Dear Justin,
> >>
> >> Thank you for your reply.
> >> It is a protein carbohydrate system. Including the solvent, the number
> of
> >> atoms is 43499.
> >> I have minimized the system for 200 ps followed by NPT and NVT
> simulations
> >> for 200 ps respectively
> >>
> >>
> > Given that your temperature output started from 0 K, then you did not
> > continue from the equilibration properly by supplying the checkpoint file
> > to grompp -t. This is important to get right, otherwise you're basically
> > starting over from some random point (an equilibrated structure without
> any
> > velocities likely isn't a physically realistic state).
> >
> > Below is the .mdp file.
> >>
> >>
> >> ; VARIOUS PREPROCESSING OPTIONS
> >> title = REMD Simulation
> >> define = -DPOSRES
> >>
> >>
> >> ; RUN CONTROL PARAMETERS
> >> integrator = md-vv ; velocity verlet algorithm -
> >> tinit = 0 ;
> >> dt = 0.002 ; timestep in ps
> >> nsteps = 5000000 ;
> >> simulation-part = 1 ; Part index is updated automatically on
> >> checkpointing
> >> comm-mode = Linear ; mode for center of mass motion
> removal
> >> nstcomm = 100 ; number of steps for center of mass
> motion
> >> removal
> >> comm-grps = Protein_Carb Water_and_Ions ; group(s)
> for
> >> center of mass motion removal
> >>
> >>
> > In a solvated system, you should not be separating these groups. This
> > could explain the sudden jump in temperature - you could have things
> > clashing badly over the course of the simulation.
> >
> >
> >
> >> ; ENERGY MINIMIZATION OPTIONS
> >> emtol = 10 ; Force tolerance
> >> emstep = 0.01 ; initial step-size
> >> niter = 20 ; Max number of iterations in relax-shells
> >> fcstep = 0 ; Step size (ps^2) for minimization of
> >> flexible constraints
> >> nstcgsteep = 1000 ; Frequency of steepest descents steps
> >> when
> >> doing CG
> >> nbfgscorr = 10
> >>
> >>
> >> ; OUTPUT CONTROL OPTIONS
> >> nstxout = 50000 ; Writing full precision coordinates
> >> every
> >> ns
> >> nstvout = 50000 ; Writing velocities every nanosecond
> >> nstfout = 0 ; Not writing forces
> >> nstlog = 5000 ; Writing to the log file every step
> 10ps
> >> nstcalcenergy = 100
> >> nstenergy = 5000 ; Writing out energy information every
> >> step 10ps
> >> nstxtcout = 2500 ; Writing coordinates every 5 ps
> >> xtc-precision = 1000
> >> xtc-grps = Protein_Carb Water_and_Ions ; subset of
> >> atoms for the .xtc file.
> >> energygrps = Protein_Carb Water_and_Ions ; Selection
> of
> >> energy groups
> >>
> >>
> >> ; NEIGHBORSEARCHING PARAMETERS
> >> nstlist = 10 ; nblist update frequency-
> >> ns-type = Grid ; ns algorithm (simple or grid)
> >> pbc = xyz ; Periodic boundary conditions: xyz,
> >> no,
> >> xy
> >> periodic-molecules = no
> >> rlist = 1.4 ; nblist cut-off
> >> rlistlong = -1 ; long-range cut-off for switched
> >> potentials
> >>
> >>
> >> ; OPTIONS FOR ELECTROSTATICS
> >> coulombtype = PME ; Method for doing electrostatics
> >> rcoulomb = 1.4 ;
> >> epsilon-r = 1 ; Relative dielectric constant for the
> >> medium
> >> pme_order = 10;
> >>
> >>
> >> ; OPTIONS FOR VDW
> >> vdw-type = Cut-off ; Method for doing Van der Waals
> >> rvdw-switch = 0 ; cut-off lengths
> >> rvdw = 1.4 ;
> >> DispCorr = EnerPres; Apply long range dispersion
> >> corrections for Energy and Pressure
> >> table-extension = 1 ; Extension of the potential lookup
> tables
> >> beyond the cut-off
> >> fourierspacing = 0.08 ; Spacing for the PME/PPPM FFT grid
> >>
> >>
> > This small Fourier spacing, coupled with the very high PME order above,
> is
> > going to unnecessarily slow your system down. Is there some reason you
> > have set these this way?
> >
> >
> >> ; GENERALIZED BORN ELECTROSTATICS
> >> gb-algorithm = Still ; Algorithm for calculating Born
> radii
> >> nstgbradii = 1 ; Frequency of calculating the Born
> radii
> >> inside rlist
> >> rgbradii = 1 ; Cutoff for Born radii calculation
> >> gb-epsilon-solvent = 80 ; Dielectric coefficient of the
> implicit
> >> solvent
> >> gb-saltconc = 0 ; Salt concentration in M for
> Generalized
> >> Born models
> >>
> >>
> >> ; Scaling factors used in the OBC GB model. Default values are OBC(II)
> >> gb-obc-alpha = 1
> >> gb-obc-beta = 0.8
> >> gb-obc-gamma = 4.85
> >> gb-dielectric-offset = 0.009
> >> sa-algorithm = Ace-approximation
> >> sa-surface-tension = -1 ; Surface tension (kJ/mol/nm^2) for the
> >> SA
> >> (nonpolar surface) part of GBSA - default -1
> >>
> >>
> > Implicit solvent should not be used if you have explicit solvent, though
> > it looks like these options are probably off since the default for the
> > implicit-solvent keyword is "no," but be aware that these are extraneous.
> >
> >
> >>
> >> ; Temperature coupling
> >> tcoupl = nose-hoover
> >> nsttcouple = 10 ;
> >> nh-chain-length = 10
> >> tc-grps = Protein_Carb Water_and_Ions ; Groups to
> >> couple separately
> >> tau-t = 10 10; Time constant (ps)-
> >>
> >
> > tau-t = 10 is far too permissive and will allow the temperature to
> > oscillate very widely. A value of 1 or so should generally be used with
> > N-H.
> >
> > -Justin
> >
> >
> > ref-t = 270.0 270.0; reference temperature (K)
> >>
> >>
> >> ; pressure coupling
> >> pcoupl = no ;-
> >>
> >>
> >> ; GENERATE VELOCITIES FOR STARTUP RUN
> >> gen-vel = no
> >> gen-temp = 270.0
> >> gen-seed = 173529
> >>
> >>
> >> ; OPTIONS FOR BONDS
> >> continuation = yes ; constrain the start configuration
> >>
> >> constraints = all-bonds
> >> constraint-algorithm = lincs ; Type of constraint algorithm-
> >> lincs-order = 4
> >> lincs-iter = 1
> >> lincs-warnangle = 30
> >>
> >>
> >> Thank you for your help.
> >>
> >> Nisha
> >>
> >>
> >>
> >> On Thu, Jun 30, 2016 at 6:21 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
> >>
> >>
> >>>
> >>> On 6/30/16 8:46 AM, NISHA Prakash wrote:
> >>>
> >>> Dear all,
> >>>>
> >>>> I have conducted a 10ns REMD simulation for a protein ligand complex
> >>>> with
> >>>> the temperature range - 270 to 350 K, however the temperature
> >>>> distribution
> >>>> plot of the replicas show that the sampling has occurred at higher
> >>>> temperatures as well that is beyond 350K -
> >>>> Below is an excerpt from the temperature xvg file
> >>>>
> >>>>
> >>>> @ title "Gromacs Energies"
> >>>> @ xaxis label "Time (ps)"
> >>>> @ yaxis label "(K)"
> >>>> @TYPE xy
> >>>> @ view 0.15, 0.15, 0.75, 0.85
> >>>> @ legend on
> >>>> @ legend box on
> >>>> @ legend loctype view
> >>>> @ legend 0.78, 0.8
> >>>> @ legend length 2
> >>>> @ s0 legend "Temperature"
> >>>> 0.000000 0.000000
> >>>> 10.000000 350.997864
> >>>> 20.000000 353.618927
> >>>> 30.000000 350.068481
> >>>> 40.000000 353.921753
> >>>> 50.000000 359.485565
> >>>> 60.000000 353.463654
> >>>> 70.000000 352.015778
> >>>> 80.000000 350.657898
> >>>> 90.000000 351.927155
> >>>> 100.000000 354.539429
> >>>> 110.000000 354.287720
> >>>> 120.000000 349.436096
> >>>> 130.000000 352.960541
> >>>> 140.000000 351.631317
> >>>> 150.000000 354.217407
> >>>> 160.000000 350.185852
> >>>> 170.000000 350.294434
> >>>> 180.000000 350.980194
> >>>> 190.000000 350.914429
> >>>> ....
> >>>> ....
> >>>> 470.000000 349.224060
> >>>> 480.000000 350.819458
> >>>> 490.000000 348.541748
> >>>> 500.000000 350.393127
> >>>> 510.000000 398.775208
> >>>> 520.000000 444.802856
> >>>> 530.000000 470.899323
> >>>> 540.000000 466.652740
> >>>> 550.000000 465.600677
> >>>> 560.000000 469.225555
> >>>> 570.000000 470.548370
> >>>> 580.000000 470.011566
> >>>> 590.000000 470.643951
> >>>> 600.000000 472.433197
> >>>> 610.000000 470.451172
> >>>> 620.000000 469.991699
> >>>> 630.000000 469.073090
> >>>> 640.000000 467.259521
> >>>> 650.000000 464.561798
> >>>> 660.000000 468.416901
> >>>> 670.000000 468.754913
> >>>> 680.000000 469.259613
> >>>> 690.000000 467.641144
> >>>> 700.000000 468.542328
> >>>>
> >>>>
> >>>> Temperature coupling was done using Nose hoover algorithm.
> >>>>
> >>>> Does this imply the sampling is wrong or insufficent?
> >>>> Any help / suggestion is appreciated.
> >>>>
> >>>>
> >>>> How large is your system, and what is it? What were your (full) .mdp
> >>> settings? The fact that your temperature started at 0 K and ramped up
> >>> suggests that you did not equilibrate prior to the run, did not
> generate
> >>> appropriate velocities, or did not continue properly. The sudden jump
> in
> >>> temperature later suggests instability, and could be due to incorrect
> >>> settings. N-H allows for large oscillations, but I wouldn't expect a
> >>> stable system to that degree.
> >>>
> >>> -Justin
> >>>
> >>> --
> >>> ==================================================
> >>>
> >>> Justin A. Lemkul, Ph.D.
> >>> Ruth L. Kirschstein NRSA Postdoctoral Fellow
> >>>
> >>> Department of Pharmaceutical Sciences
> >>> School of Pharmacy
> >>> Health Sciences Facility II, Room 629
> >>> University of Maryland, Baltimore
> >>> 20 Penn St.
> >>> Baltimore, MD 21201
> >>>
> >>> jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> >>> http://mackerell.umaryland.edu/~jalemkul
> >>>
> >>> ==================================================
> >>> --
> >>> Gromacs Users mailing list
> >>>
> >>> * Please search the archive at
> >>> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> >>> posting!
> >>>
> >>> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >>>
> >>> * For (un)subscribe requests visit
> >>> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> >>> send a mail to gmx-users-request at gromacs.org.
> >>>
> >>>
> > --
> > ==================================================
> >
> > Justin A. Lemkul, Ph.D.
> > Ruth L. Kirschstein NRSA Postdoctoral Fellow
> >
> > Department of Pharmaceutical Sciences
> > School of Pharmacy
> > Health Sciences Facility II, Room 629
> > University of Maryland, Baltimore
> > 20 Penn St.
> > Baltimore, MD 21201
> >
> > jalemkul at outerbanks.umaryland.edu | (410) 706-7441
> > http://mackerell.umaryland.edu/~jalemkul
> >
> > ==================================================
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
More information about the gromacs.org_gmx-users
mailing list