[gmx-users] pressure fluctuations

[Justin A. Lemkul]

Hassan Shallal wrote:
> Dear Gromacs users,
>  
> I have some concerns about the both the pressure fluctuations and 
> averages I obtained during the equilibration phase. I have already read 
> through several similar posts as well as the following link 
> http://www.gromacs.org/Documentation/Terminology/Pressure. I understand 
> the pressure is a macroscopic rather than instantaneous property and the 
> average is what really matters. I also found out through similar posts 
> that negative average pressure indicates the system tendency to contract.
>  
> In the above link, it mentioned that pressure fluctuations should 
> decrease significantly with increasing the system's size. In my cases, I 
> have a fairly big systems (case_1 with *17393* water molecules 
> and case_2 with *11946 *water molecules). However, the pressure still 
> has huge fluctuations (around 500 bars) from the reference value (1 
> bar). Here are the average pressure and density values resulting from 
> the equilibration phases of two cases, please notice the negative 
> average pressure values in both cases...
>  
> Case_1_pressure:
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Pressure                   *-2.48342*       0.92    369.709   -4.89668  
> (bar)
> Case_1_density:
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Density                     1022.89       0.38     3.8253    2.36724  
> (kg/m^3)
> Case_2_pressure:
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Pressure                   *-8.25259*        2.6    423.681   -12.1722  
> (bar)
> Case_2_density:
> Energy                      Average   Err.Est.       RMSD  Tot-Drift
> -------------------------------------------------------------------------------
> Density                     1034.11       0.37    2.49964    1.35551  
> (kg/m^3)
>  
> So I have some questions to address my concerns:
> 1- each of the above systems has a protein molecule, NaCl to give 0.15 M 
> system and solvent (water) molecules... Could that tendency to contract 
> be an artifact of buffering the system with sodium and chloride ions?
>  

I suppose anything is possible, but given that these are fairly standard 
conditions for most simulations, I tend to doubt it.  My own (similar) systems 
do not show this problem.

> 2- how to deal with the tendency of my system to contract?  Should 
> I change the number of water molecules in the system?  
> or
> Is it possible to improve the average pressure of the above systems by 
> increasing the time of equilibration from 100 ps to may be 500 ps or 
> even 1 ns?
>  
> 3- Is there a widely used range of average pressure (for ref_p = 1 bar) 
> that indicates acceptable equilibration of the system prior to the 
> production?
>  

To answer #2 and #3 simultaneously - equilibration is considered "finished" when 
your system stabilizes at the appropriate conditions (usually temperature and 
pressure).  Your results indicate that your equilibrium is insufficient.

> 4- I can't understand how the system has a tendency to contract whereas 
> the average density of the solvent is already slightly higher than it 
> should be (1000 kg/m^3).

The contraction causes the density to rise.  Pressure and density are not 
independent; density is a result of pressure.

> I would like to ignore the pressure based judgement of the above 
> equilibration given that the average density values are very close to 
> the natural value (1000 kg/m^3) (by the way I am using tip3p water model 
> with CHARMM27 ff) Any comment!!
>  

It is not guaranteed that simulation water models will reproduce the real 
(experimental) density of water.  If memory serves, the expected density of 
TIP3P should be ~0.98 g mL^{-1}, but I could be wrong.

> 5- Is the huge fluctuation of the pressure values of the above system 
> despite thier large sizes still acceptable? or large fluctuation is only 
> acceptable for small size systems and is unacceptable for large size 
> systems?
> If it is unacceptable, any idea of how could it be alleviated or minimized?
>  

The fluctuations seem reasonable.  You might try equilibrating with the 
Berendsen barostat first, then switching to Parrinello-Rahman.  The P-R barostat 
allows for wider fluctuations, so if the system is poorly equilibrated, your 
system will be slower to converge.

-Justin

> I am including the .mdp used in the above equilibration in case it is 
> needed.
>  
> Any feedback or response to the above questions is so much appreciated..
>  
> Great regards
> Hassan
>  
> .mdp used for the above equilibration
> define  = -DPOSRES ; position restrain the protein
> ; Run parameters
> integrator = md  ; leap-frog integrator
> nsteps  = 25000  ; 4 * 25000 = 100 ps
> dt  = 0.004  ; 4 fs, virtual sites along with heavy hydrogens are used
> ; Output control
> nstxout  = 100  ; save coordinates every 0.2 ps
> nstvout  = 100  ; save velocities every 0.2 ps
> nstenergy = 100  ; save energies every 0.2 ps
> nstlog  = 100  ; update log file every 0.2 ps
> ; Bond parameters
> continuation = yes  ; Restarting after NVT
> constraint_algorithm = lincs ; holonomic constraints
> constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
> lincs_iter = 1  ; accuracy of LINCS
> lincs_order = 6  ; also related to accuracy, changed to 6 because of 
> using virtual sites along with a larger time step
> ; Neighborsearching
> ns_type  = grid  ; search neighboring grid cells
> nstlist  = 5  ; 20 fs
> rlist  = 1.2  ; short-range neighborlist cutoff, equal to rcoulomb to 
> allow for PME electrostatics  (in nm)
> ; Lennard-Jones
> vdwtype         = switch        ; VDW interactions are switched of 
> between 1 and 1.2
> rvdw_switch     = 1             ;
> rvdw            = 1.2           ; short-range vdw cutoff, optimal for 
> CHARMM27 ff  (in nm)
> ; Electrostatics
> coulombtype = PME    ; Particle Mesh Ewald for long-range electrostatics
> pme_order = 4  ; cubic interpolation
> rcoulomb = 1.2  ; short-range electrostatic cutoff, optimal for CHARMM27 
> ff  (in nm)
> ; Temperature coupling is on
> tcoupl  = V-rescale ; modified Berendsen thermostat
> tc-grps  = Protein Non-Protein ; two coupling groups - more accurate
> tau_t  = 0.1 0.1 ; time constant, in ps
> ref_t  = 300  300 ; reference temperature, one for each group, in K
> ; Pressure coupling is on
> pcoupl  = Parrinello-Rahman ; Pressure coupling on in NPT
> pcoupltype = isotropic ; uniform scaling of box vectors
> tau_p  = 1  ; 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 switch vdW scheme
> ; Velocity generation
> gen_vel  = no  ; Velocity generation is off
>  
> 

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

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

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