[gmx-users] How the pressure is computed

[Andrew DeYoung]

Hi,

I am running Gromacs 4.5.5.  I am using an orthorhombic simulation box (all
box angles of 90 degrees: i.e., a rectangular prism) for which two of the
box vectors are approximately equal in length, but the third box vector is
much, much longer than the first two vectors:

Lz >> Lx ~ Ly

For example, my box vectors are Lz = 30 nm and Lx ~ 3 nm and Ly ~ 3 nm, so
my simulation box has approximate volume 270 nm^3.  However, my system
occupies only a small fraction of the volume of my simulation box: my system
occupies only about (6 nm)*(3 nm)*(3 nm) = 54 nm^3.  My system is terminated
at z=-3 nm and at z=+3 nm by large, periodic molecules (comprising a
membrane).  Between these barriers of large, periodic molecules are many
small molecules (comprising a liquid).  

The upshot of this is that a very large percentage of the volume is vacuum,
having no molecules.  The reason for this is a trick, by which my system is
effectively two-dimensional; Lz is so long so that the system will not "see"
its own periodic image along the z direction.  

But, my question is, what if I run an NPT simulation and then use g_energy
to print the system pressure as a function of time?  Is the calculated
pressure: (1) the pressure calculated by taking into account the entire box
volume (270 nm^3), or is it (2) the pressure calculated by taking into
account only the volume of the box which is not vacuum, i.e., actually
occupied by molecules (54 nm^3)?

My guess is that it is pressure (1) that is being calculated, and printed to
the .edr file.  Is this correct?

What if I would like to calculate the pressure (2), the pressure due to only
the non-vacuum part of the system?  Is there a way to specify that only part
of the simulation box be taken into account when computing the pressure?  Or
is there a way to calculate the "local pressure," just for that 54 nm^3 out
of 270 nm^3?  

Thank you very much for your time!

Andrew DeYoung
Carnegie Mellon University