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.

Reply via email to