[gmx-users] Pressure Calculation in Gromacs
jalemkul at vt.edu
Tue Jul 17 04:23:31 CEST 2018
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.
> 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?
Justin A. Lemkul, Ph.D.
Virginia Tech Department of Biochemistry
303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061
jalemkul at vt.edu | (540) 231-3129
More information about the gromacs.org_gmx-users