[gmx-developers] Pressure calculation in Gromacs

[Dmitri Rozmanov]

Hello All,

I am comparing the results of Gromacs with our in-house MD simulator. It seems that I have 
general agreement with energies, and forces, but the pressure is somehow different.

I am doing one point calculation of a system of 2400 TIP4P water molecules in a rectangular 
elongated box.

I was actually able to separate most of the pressure contributions Gromacs calculates and
have some troubles with understanding how the shake/links/settles virials are obtained and used.

All force virials are the same, the trouble is with the constraints virial.

Our code is for rigid dynamics and it explicitly uses translations and rotations of the 
molecules, waters in this case. To debug the pressure I removed the charges from TIP4P sites 
and zeroed LJ params as well. So, there are no nonbonded interactions in the system at all, and 
all the molecules are rigid bodies.

I can assume, that due to the equi-partitioning the kinetic energy is split between the 
rotations and translations equally and only the translational part of it will contribute to the 
pressure. The rotational kinetic virial should be canceled in the pressure due to constraints 
in the system. Right?

So, we should have someting like this: P = 16.6 * 2 * Ekin_trn / 3 / V

In my topology I have
	EKinetic Trn (kJ/mol) = 8709.4539759232
	EKinetic Rot (kJ/mol) = 8688.6908125891
	EKinetic Tot (kJ/mol) = 17398.1447885123

Wich gives me P = 1273.3 bar

This is the number I get in Gromacs if I remove the rotations from the starting topology:
	
KE = 8709.45 kJ/mol --> P = 1273.3 bar.

However, I both the translations and rotations are retained I get this:

KE = 17398.1 kJ/mol --> P = 1883.9 bar = 2543.6 - 659.7

Where, the first is the pressure due to the total KE, and the second is the pressure due to the 
constraints. The kinetic part is correctly about twice the translational only contribution, but
the constraint virial does not cancel the rotational virial!

The difference is about 600 bar... Well...

If all the interactions are there, it becomes even worse:

My calculation P = -112.3 bar, Gromacs gives -2238.98 bar

The pressure is negative due to the sightly different water model used to obtain the topology. 
This is normal.

Now I have 2100 bar difference.

I changed the tip4p.itp in a way that Gromacs uses shake/lincs instead of settles, and the 
pressures are almost the same for all three of them.

So, am I doing something wrong? Which pressure is more correct? If Gromacs pressure is not 
correct, is it a big deal?

Thank you very much. Any suggestions will be appreciated.

With best wishes and Happy New Year!

	--- Dmitri.