[gmx-developers] Pressure calculation in Gromacs
Dmitri Rozmanov
rozmanov at gmail.com
Fri Jan 2 11:03:12 CET 2009
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.
More information about the gromacs.org_gmx-developers
mailing list