[gmx-users] Dubious results with NPT
sujithkakkat .
sujithks58 at gmail.com
Mon Feb 24 05:37:32 CET 2014
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 <vvchaban at 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 <jalemkul at 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
> >
> > 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.
> --
> 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.
>
More information about the gromacs.org_gmx-users
mailing list