Hello, Thank you both for the comments. I am using gromos96 forcefield . I read a little bit and as you said the nonbonded cutoff has to be higher. The tau_p=5ps was chosen , since the manual mentions that the value has to be raised by 4-5 times on going from berendsen to parrinello-rahman barostat, though I did not completely follow the reasons behind it. I will try with lower values.
I had run an earlier simulation with the same parameters for a pure water system for 5ns and reference pressure 5bar, and things worked fine there with the average pressure at 5.15bar. I guess the sigma for the case of water was low and therefore the small cut-off of 0.8nm did not matter. However the case of cyclohexane alone remains to be tried. I guess Dr Vitaly was saying about using a switch/shift function. I will try the simulation with the new settings and see. Regards, Sujith. On Sun, Feb 23, 2014 at 10:37 PM, Dr. Vitaly Chaban <vvcha...@gmail.com>wrote: > There is such thing in simulations as energy conservation... > > If you use "vdwtype= cut-off" and this cut-off happens at 0.8nm, while > sigma for the largest atom is ~0.34nm, the problems are inevitable. > > Your cut-off should not be smaller than 0.90nm, and you need to apply > a more polite method to bring pairwise energy term down to zero at > this cut-off. > > > Dr. Vitaly V. Chaban > > > On Sun, Feb 23, 2014 at 3:07 PM, Justin Lemkul <jalem...@vt.edu> wrote: > > > > > > On 2/23/14, 8:30 AM, sujith wrote: > >> > >> Dear GROMACS users, > >> > >> I am new to GROMACS, and recently started using the version 4.6.5. > >> > >> I have seen a lot of NPT related issues raised earlier in this forum, > but > >> in > >> my case the error looks much more severe. > >> > >> I am following Justin Lemkul's tutorial > >> > >> ( > http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/biphasic/01_genconf.html > ) > >> on Biphasic system, containing cyclohexane and water.Everything went > well > >> till the NPT simulation to bring system to reference pressure of 1 bar. > >> > >> Here are the details of the system , .mdp file , what I did and where I > >> find the problem. > >> > >> SYSTEM : 512 cyclohexane + 2720 water molecules. > >> > >> CURRENT STATUS: > >> > >> (1) Energy minimization : energy converged. > >> > >> (2) NVT (tcoupl=berendsen , tau_t=0.1ps / ref_t=298K) > >> completed and target pressure of 298 K attained. > >> > >> (3) NPT (tcoupl=nose-hoover, tau_t=0.5ps / > >> pcoupl=berendsen, > >> tau_p=1ps / ref_p = 1bar , total time=15 ns) completed with the > following > >> result: > >> > >> Energy Average Err.Est. RMSD > >> Tot-Drift > >> > >> > >> > ------------------------------------------------------------------------------- > >> Pressure 1.08924 0.67 114.66 > >> -1.74033 (bar) > >> > >> Since berendsen barostat doesn't generate the true ensemble, an > >> equilibration with pcoupl=parrinello-rahman is performed in the next > >> step. > >> > >> (4) NPT (pcoupl=parrinello-rahman, tau_p=5ps, total > >> time=5ns) > >> ----> PROBLEM FACED: The average pressure too high as shown below > which I > >> feel is not going to improve. > >> > >> Energy Average Err.Est. RMSD Tot-Drift > >> > >> > ------------------------------------------------------------------------------- > >> Pressure 23.6651 0.53 144.464 -1.00634 > >> (bar) > >> > >> > >> > >> I am aware of the fact that pressure is subject to large fluctuations > >> in > >> small sized systems and that this may affect the average value of the > >> pressure. But , here the average pressure looks too large to be ignored. > >> The > >> pressure-vs-time graph doesn't show any upward trend, and the pressure > >> looks > >> like fluctuating about a central value. > >> > >> > >> Here is the .mdp file. Only changes made from the > >> previous .mdp for berendsen pressure coupling , are with pcoupl and > >> tau_p. > >> > >> ; 7.3.3 Run Control > >> integrator = md > >> tinit = 0 > >> dt = 0.002 > >> nsteps = 2500000 > >> comm_mode = Linear > >> nstcomm = 1 > >> comm_grps = CHX SOL > >> > >> ; 7.3.8 Output Control > >> nstxout = 2500 > >> nstvout = 2500 > >> nstfout = 2500 > >> nstlog = 2500 > >> nstenergy = 100 > >> nstxtcout = 1000 > >> xtc_precision = 1000 > >> xtc_grps = System > >> energygrps = System > >> > >> ; 7.3.9 Neighbor Searching > >> nstlist = 1 > >> ns_type = grid > >> pbc = xyz > >> rlist = 0.8 > >> > >> ; 7.3.10 Electrostatics > >> coulombtype = PME > >> rcoulomb = 0.8 > >> ; 7.3.11 VdW > >> vdwtype = cut-off > >> rvdw = 0.8 > > > > > > What force field are you using? If it's Gromos96 like my tutorial, the > > value of rvdw is much too short and can lead to artifacts. Using 0.8 for > > rlist/rcoulomb is fine, though 0.9 is more common. rvdw should be 1.4. > > > >> DispCorr = EnerPres > >> ; 7.3.13 Ewald > >> fourierspacing = 0.12 > >> pme_order = 4 > >> ewald_rtol = 1e-5 > >> > >> ; 7.3.14 Temperature Coupling > >> tcoupl = nose-hoover > >> tc_grps = CHX SOL > >> tau_t = 0.5 0.5 > >> ref_t = 298 298 > >> > >> ; 7.3.15 Pressure Coupling > >> pcoupl = parrinello-rahman > >> pcoupltype = isotropic > >> tau_p = 5.0 > >> compressibility = 4.5e-5 > >> ref_p = 1.0 > >> > > > > The value of tau_p seems a bit long to me; does changing it to 1.0 > improve > > the results? > > > >> ; 7.3.17 Velocity Generation > >> gen_vel = no > >> > >> ; 7.3.18 Bonds > >> constraints = all-bonds > >> constraint_algorithm = LINCS > >> continuation = yes > >> lincs_order = 4 > >> lincs_iter = 1 > >> lincs_warnangle = 30 > >> "npt.mdp" 63L, 4080C > >> > >> I guess there is something seriously wrong in the choice of > >> methods/parameters in the .mdp file, which I cant figure out. Kindly go > >> through and let me know your comments. > >> I would be happy to give any further details. Any help would be > >> appreciated. > > > > > > An even better test is to simplify the system. Run pure water or pure > > cyclohexane with the existing settings, then try the modifications above. > > That should help root out whether you're having a problem with the .mdp > > settings, force field, or maybe even the combination of both. > > > > -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. > -- > 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. > -- 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.