[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