[gmx-users] REMD - replicas sampling in temperatures beyond the assigned range
NISHA Prakash
nishnith20591 at gmail.com
Thu Jun 30 16:15:42 CEST 2016
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.
>
More information about the gromacs.org_gmx-users
mailing list