Large pressure fluctuations in molecular dynamics

Question

I'm performing MD simulations, using NPT ensemble. The problem is that the pressure fluctuates severely and doesn't average at required level.

For example, I set the pressure at 1 bar and it averages around -60 bar jumping from -500 to 500 bar during the run.

I'm using Parrinello-Rahman barostat with $\tau$ parameter equal to 20 ps. I've tried to adjust $\tau$, but nothing changed.

Here is the part of mdp file:

Pcoupl                   = Parrinello-Rahman
Pcoupltype               = isotropic
tau_p                    = 20
compressibility          = 3e-5
ref_p                    = 1.0

Is this an important thing to consider if density and temperature looks OK?

Answer 1

Large pressure fluctuations are normal for nanoscale simulations of aqueous solutions (and other condensed systems without a gas phase).

Pressure fluctuations are usually enormous in simulations of aqueous solutions because water is nearly incompressible.

Let's say your system is a cube of water $(2~\mathrm{nm})^3$. Smaller systems will have larger fluctuations. At equilibrium, the variance of the pressure (its fluctuation) is

$$\tag{1}\sigma_P^2 = -kT \left( \frac{\partial P}{\partial V} \right)_S = \frac{kT}{\kappa_S V},$$

where $\kappa_S$ is the isentropic compressibility. The isentropic compressibility is related to the density $\rho$ and speed of sound ($c$) by $\kappa_S = 1/(\rho c^2)$, so

$$\tag{2}\sigma_P^2 = \frac{\rho c^2 k T}{V}.$$

The standard deviation of the pressure should then be $$\tag{3}\sigma_P = \sqrt{\frac{\rho c^2 k T}{V}}.$$

Taking the values for water near standard conditions $\rho=997~\mathrm{kg/m^3}, c=1500~\mathrm{m/s}, T=298.15~\mathrm{K}$ and using GNU units,

units 'sqrt(997*kg/m^3*(1500*m/s)^2*k*298*K/((2*nm)^3))' 'bar'

I get $\sigma_P \approx 340~\mathrm{bar}$. So the fluctuation of $\pm500~\mathrm{bar}$ is pretty reasonable.

You say you get a mean pressure of $-60~\mathrm{bar}$. What is the uncertainty in that calculation? The uncertainty of the mean pressure can be estimated by $\delta P = \sigma_P/\sqrt{N}$, where $N$ is the number of independent data points for the pressure. So the number of independent (uncorrelated) data points you need to get to get the uncertainty $\delta P < 60~\mathrm{bar}$ is given by $N=\sigma_P^2/\delta P^2 \approx 30$. How long was your simulation? What is the autocorrelation time of the pressure (which determines the number of independent pressure values you have).

My guess is that your results are statistically consistent with a mean pressure of $1~\mathrm{bar}$.

Finally, for most molecular interactions in water, pressure–volume work ($p~\mathrm{d}V$) is negligible, so you don't have to worry much about the exact pressure as long as you don't have any weird vacuum bubbles that sometimes appear in simulations.

Aside: You might say, I'm using a barostat, shouldn't it force the pressure be constant? First, you have to realize that the instantaneous pressure is not well-defined and there isn't a unique way to calculate it. Some barostats might make some measures of instantaneous pressure constant, but this isn't true for the barostat I usually use. The likewise is true for thermostats. Instantaneous temperature is usually calculated from the average kinetic energy of a single simulation frame, which doesn't necessarily correspond to the thermodynamic temperature. Consider that a system with no momentum degrees of freedom (and no kinetic energy) can still have a thermodynamic temperature (e.g. the Ising model).

Answer 2

Thermodynamics is being bossy?

After discussions in the chatroom I don't see anything wrong with your simulation settings. I think the culprit is phase equilibria.

Either your system wants to be in two phases, or, it wants to be a single phase but you are starting with a box size large enough (density low enough) that you start with two phases and even after 10-50 ns runs, do not truly equilibrate.

If the problem is the latter, that you are starting with two phases and not converging to 1, the solution is simple: start with a higher density so that you start as a liquid system, and then hopefully, never leave - and this was noted by WaterMolecule.

If the problem is that at P=1 bar your system wants to be in two phases, that may explain why the pressure is so wild. An NPT should be able to accommodate multiple phases, but, it may take a long time. I am not a big fan of two phases in 1 simulation box, but your molecules are likely too big for Gibbs Ensemble, and the best method for MD phase equilibria won't be published by some colleagues of mine for a couple more months.

An interesting question is whether there should be 2 phases or 1 at these conditions? The force-field appears to think that state point is at the least, close to a phase boundary. Force-fields are often wrong though, so that in itself would be good for you to know.