[gmx-users] Dubious results with NPT
Dr. Vitaly Chaban
vvchaban at gmail.com
Sun Feb 23 18:07:51 CET 2014
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.
More information about the gromacs.org_gmx-users
mailing list