[gmx-users] Cannot get correct pressure value with MTTK pressure coupling

[Justin A. Lemkul]

On 6/6/12 11:07 AM, Bao Kai wrote:
> HI, all,
>
> I want to use MTTK for the pressure coupling for NPT simulation, while
> I can not get the correct pressure value.
>
> The following result is a simulation with 1000 SPC-E water molecules
> with temperature 303.15K and pressure 1 bar.
>
> Statistics over 1000001 steps [ 0.0000 through 1000.0000 ps ], 3 data sets
> All statistics are over 200001 points
>
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Temperature                 303.155      0.016    5.54033  0.0532334  (K)
> Pressure                    138.295        1.3    429.699   -5.68756  (bar)
> Density                     1002.82       0.24    8.04912   0.668867  (kg/m^3)
> .
>
> The following is the mdp file
>
> title = fizzy water
> ; Run parameters
> integrator = md-vv ; velocity-verlet integrator
> nsteps = 1000000; 0.001 * 1000000 = 1000 ps, 1 ns
> dt = 0.001 ; 2 fs
>
> ; Output control
> nstxout = 1000 ; save coordinates every 2 ps. Incr by *10
> nstvout = 1000 ; save velocities every 2 ps. Incr by *10
> nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps. Incr by *10
> nstenergy = 1000 ; save energies every 2 ps. Incr by *10
> nstlog = 1000 ; update log file every 2 ps. Incr by *10
>
> ; Bond parameters
> continuation = yes ;<-- Restarting after NPT
> constraint_algorithm = lincs ; holonomic constraints
> constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
> lincs_iter = 1 ; accuracy of LINCS
> lincs_order = 4 ; also related to accuracy
>
> ; Neighborsearching
> ns_type = grid ; search neighboring grid cells
> nstlist = 5 ; 10 fs
> rlist = 0.9 ; short-range neighborlist cutoff (in nm)
> rcoulomb = 0.9 ; short-range electrostatic cutoff (in nm)
> rvdw = 0.9 ; short-range van der Waals cutoff (in nm)
>
> ; Electrostatics
> coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
> pme_order = 4 ; cubic interpolation
> fourierspacing = 0.16 ; grid spacing for FFT
>
> ; Temperature coupling is on
> tcoupl = nose-hoover ; modified Berendsen thermostat
> tc-grps = System ; two coupling groups - more accurate
> tau_t = 0.1 ; time constant, in ps
> ref_t = 303.15 ; reference temperature, one for each group, in K
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; turning on pressure coupling
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;<-- Pressure coupling is on
> pcoupl = mttk; Pressure coupling on in NPT
> nstpcouple  = 1;
> pcoupltype = isotropic ; uniform scaling of box vectors
> tau_p = 2.0 ; time constant, in ps
> ref_p = 1.0 ; reference pressure, in bar
> compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
>
> ; Periodic boundary conditions
> pbc = xyz ; 3-D PBC
>
> ; Dispersion correction
> DispCorr = EnerPres ; account for cut-off vdW scheme
>
> ; Velocity generation
> gen_vel = no ;<-- Velocity generation is off
> ;
>
> Could anybody tell me anything wrong with my configuration?
>

Probably nothing.  Pressure is a very hard term to converge, as noted here:

http://www.gromacs.org/Documentation/Terminology/Pressure

Note that the manual warns that you will get very large oscillations in the 
pressure value if you don't start from a reasonable (i.e., previously 
equilibrated) value.  Using Berendsen for initial NPT equilibration is more 
common and more stable.

Your simulations are also very short.  Pressure (and other terms) may require 
more than just 1 ns to converge, but you should also bear in mind that the value 
you obtained is 138.295 ± 429.699.  Is that significantly different from the 
target value?  Given the magnitude of the standard deviation, I would say no.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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