[gmx-users] Md simulation mdp options error

[Justin Lemkul]

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

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

==================================================