[gmx-developers] Pressure calculation in Gromacs

[David van der Spoel]

Dmitri Rozmanov wrote:
> 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?
> 

gromacs pressure is correct. but you should compare identical models.
check the description in the appendix, and note that gromacs uses an 
atomic virial rather than a molecular virial which your code might do.
note that you can also inspect the virial from the energy file.

> Thank you very much. Any suggestions will be appreciated.
> 
> With best wishes and Happy New Year!
> 
>     --- Dmitri.
> 
> 
> 
> 
> 
> 
> _______________________________________________
> gmx-developers mailing list
> gmx-developers at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-developers
> Please don't post (un)subscribe requests to the list. Use the www 
> interface or send it to gmx-developers-request at gromacs.org.


-- 
David van der Spoel, Ph.D., Professor of Biology
Molec. Biophys. group, Dept. of Cell & Molec. Biol., Uppsala University.
Box 596, 75124 Uppsala, Sweden. Phone:	+46184714205. Fax: +4618511755.
spoel at xray.bmc.uu.se	spoel at gromacs.org   http://folding.bmc.uu.se