[gmx-users] REMD - replicas sampling in temperatures beyond the assigned range
Justin Lemkul
jalemkul at vt.edu
Thu Jun 30 15:25:12 CEST 2016
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
==================================================
More information about the gromacs.org_gmx-users
mailing list