On 4/18/14, 4:57 AM, Ankita Naithani wrote:
Hi,
I am trying to run a simulation with the following optins in my production
mdp file. Previously I have run the simulations with the production run
parameters using V-rescale thermostat and berendsen barostat. But now
recently I read in discussion that Justin had advised using
Parinello-Rahman barostat as it would be much more suitable and provides
adequate sampling. Also, I had to use the Verlet cut-off scheme because the
cluster where I run my production simulation has updated gomacs to the
recent version 4.6.5 which requires one to use the verlet scheme.
A few things to clarify:
1. Parrinello-Rahman is superior to Berendsen for production runs, not
necessarily equilibration. In fact, P-R is often unstable during equilibration
(see below).
2. The Verlet scheme is *not* required in version 4.6.5. It is the new default,
as the group scheme will eventually be phased out, but that is not the case at
the moment.
My mdp file is below
title = production MD
; Run parameters
integrator = md ; leap-frog algorithm
nsteps = 50000000 ; 0.002 * 50000000 = 100000 ps or 100 ns
dt = 0.002 ; 2 fs
; Output control
nstxout = 0 ; save coordinates every 2 ps
nstvout = 0 ; save velocities every 2 ps
nstxtcout = 1000 ; xtc compressed trajectory output every 5 ps
nstenergy = 1000 ; save energies every 5 ps
nstlog = 1000 ; update log file every 5 ps
; Bond parameters
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds)
constraine
d
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cells
nstlist = 5 ; 25 fs
rlist = 1.0 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.4 ; short-range van der Waals cutoff (in nm)
rlistlong = 1.0 ; long-range neighborlist cutoff (in nm)
cutoff-scheme = Verlet
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostat
ics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
nstcomm = 10 ; remove com every 10 steps
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more
accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 318 318 ; reference temperature, one for each group,
in
K
; Pressure coupling is off
pcoupl = Parrinello-Rahman ; pressure coupling is on
for NP
T
pcoupltype = isotropic ; uniform scaling of box vectors
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is on
gen_temp = 318 ; reference temperature, for protein in K
ERROR 1 [file md.mdp]:
With Verlet lists rcoulomb!=rvdw is not supported
NOTE 1 [file md.mdp]:
With Verlet lists the optimal nstlist is >= 10, with GPUs >= 20. Note
that with the Verlet scheme, nstlist has no effect on the accuracy of
your simulation.
NOTE 2 [file md.mdp]:
nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting
nstcomm to nstcalcenergy
WARNING 1 [file md.mdp]:
You are generating velocities so I am assuming you are equilibrating a
system. You are using Parrinello-Rahman pressure coupling, but this can
be unstable for equilibration. If your system crashes, try equilibrating
first with Berendsen pressure coupling. If you are not equilibrating the
system, you can probably ignore this warning.
Velocities were taken from a Maxwell distribution at 318 K
Removing all charge groups because cutoff-scheme=Verlet
There were 2 notes
There was 1 warning
-------------------------------------------------------
Program grompp_d, VERSION 4.6.5
Source code file:
/home/y07/y07/gmx/4.6.5-phase1/source/src/kernel/grompp.c, line: 1610
Fatal error:
There was 1 error in input file(s)
For more information and tips for troubleshooting, please check the GROMACS
website at http://www.gromacs.org/Documentation/Errors
-------------------------------------------------------
So, I guess, NOTES 1 and 2 are harmless and can be ignored, right?
Correct. It's merely advisory information for rather low-level settings.
Regarding error 1, I think I would have to change my rvdw to 1.0 and that
should sort the problem out.
You should be setting your cutoffs based on what the force field requires.
I wanted to get advice on the warning too, should I ignore and proceed? Or
should I go back to berendsen as I have been using for my previous
simulations. For NVT and NPT, I have used berendsen. As this is a very
crucial simulation run, I just wanted to confirm since I have not had any
previous experience with Parinello-Rahman so, any advice would surely help
me.
If, as the title implies, this is a production simulation, you shouldn't be
re-generating velocities. Doing so negates the purpose of prior equilibration
completely. Using P-R with velocity generation is unstable, because the
velocities can push things in bad directions. To preserve the previous
equilibration, you should be passing a .cpt file to grompp -t and setting
"gen_vel = no."
-Justin
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 601
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalem...@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-requ...@gromacs.org.
