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 <jalem...@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 > > 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. > -- Ankita Naithani -- 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.