# Large pressure fluctuations in molecular dynamics

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?

• I've moved the comments on this question to the GROMACS chat room
– Tyberius
Apr 1, 2021 at 15:10

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).

• gromacs usually gets the average pressure equal to the set pressure though. The Berendsen barostat for instance basically forces the simulation to the set pressure. You should join the chat, linked above below OP's question. More trouble shooters is better :) Apr 1, 2021 at 18:27
• Thank you for your explanation by the means of statistics. I have a polymer system, not aqueous solution, but the compressibility is near the same as for water. Also I have weird "vacuum bubbles" in the center of my box. Apr 2, 2021 at 9:44
• @romaichenko The weird vacuum bubbles usually result from starting the system with too large of an initial volume. Try starting with a smaller volume. You may want to use an initial structure with a higher density as well. Apr 2, 2021 at 17:22
• @B.Kelly I'm used to seeing the instantaneous measures of temperature and pressure fluctuate a lot and only converge after long simulations (hundreds of nanoseconds?), but maybe this is due to the methods I'm using (Langevin thermostat/Langevin piston). Instantaneous measures of temperature and pressure usually (and by necessity) consider only some degrees of freedom (momentum) and therefore should fluctuate even in an ideal NPT or NVT ensemble. Questions similar to the OQ appear often on MD forums, so I don't think this discussion necessarily belongs in Gromacs chat. Apr 2, 2021 at 17:36
• @WaterMolecule sure its a general problem, but the chat is where its at, I didn't make it. I also use Langevin, since I am after free energies. Gromacs uses the atomic virial to calculate pressure, since bond forces are included, and they fluctuate alot, the virial fluctuates alot. So yes, fluctations are normal, but, being 60 bar away from set pressure is not normal, even if it is within uncertainty, that is too large. Apr 2, 2021 at 17:50

## 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.