[gmx-users] Pressure Calculation in Gromacs

[Justin Lemkul]

On 7/16/18 9:55 PM, Brandon Wiebe wrote:
> Hi there,
>
> I am wondering if anyone can break down exactly how the instantaneous pressure of any given frame is calculated in Gromacs?
>
> According to section 3.4.3 of the manual, the pressure tensor is given by: P = (2/V) * (Ek - W), where Ek is kinetic energy and W is the virial.
> I have taken a particular frame and calculated the volume, kinetic energy and virial.
> My calculated quantities all match quite closely with the values reported by Gromacs, however, the pressure that I get is significantly lower than the value reported by Gromacs.
> The elements of my calculated pressure tensor are all approximately 16.5 times smaller than those reported by Gromacs, even though the component quantities that go into the pressure calculation are nearly identical.

There's a unit conversion factor (see src/gromacs/math/units.h):

/* Pressure in MD units is:
  * 1 bar = 1e5 Pa = 1e5 kg m^-1 s^-2 = 1e-28 kg nm^-1 ps^-2 = 1e-28 / 
AMU amu nm^1 ps ^2
  */
#define BAR_MDUNITS      (1e5*NANO*PICO*PICO/AMU)
#define PRESFAC          (1.0/BAR_MDUNITS)

PRESFAC works out to be 16.6.

-Justin

> I have repeated this calculation for a number of frames, and have tried using Gromacs 4.6.7 and 5.0.4 to no avail.
> Clearly there is some systematic error in my methodology, where could this discrepancy come from?
>
> Regards,
> Brandon

-- 
==================================================

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

==================================================