[gmx-users] Dubious results with NPT

[Justin Lemkul]

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

==================================================