[gmx-users] Md simulation mdp options error
Ankita Naithani
ankitanaithani at gmail.com
Fri Apr 18 14:55:51 CEST 2014
Thank you Justin for your reply. Yes indeed this is a production run. So
would the following mdp file be okay now?
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.0 ; 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 = no ; Velocity generation is off
gen_temp = 318 ; reference temperature, for protein in K
grompp -f md.mdp -c npt.gro -p complex.top -o md.tpr?
Could you please advice?
Kind regards,
Ankita
On Fri, Apr 18, 2014 at 12:58 PM, Justin Lemkul <jalemkul at vt.edu> wrote:
>
>
> 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
>
> ==================================================
> --
> 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.
>
--
Ankita Naithani
More information about the gromacs.org_gmx-users
mailing list