[gmx-users] How the pressure is computed

[Justin A. Lemkul]

Andrew DeYoung wrote:
> 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?  
> 

I suspect that within very few simulation steps, there will be no difference. 
Under NPT, vacuum regions compress very rapidly.

You may need a different approach, for instance "pbc = xy" for an infinite 
z-dimension.  There are limitations to this approach, as stated in the manual, 
but it may be more sound.

-Justin

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

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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