[gmx-users] Re: How the pressure is computed

[Andrew DeYoung]

Justin, thank you very much for your time and help.  I am using pbc = xy
(with two walls, one on each side), but my box actually expands wildly,
which is very strange, rather than contracting and compressing the vacuum
regions.  So I am not sure what to think.  Maybe I should switch to pbc =
xyz and see what happens with the box size under NPT.  

If you have time, I have another, related question.  I assume that the
density stored in the .edr file (and printed using g_energy) takes into
account the entire box volume, not just the part of the box which is not
vacuum.  Is there any way to compute the density locally, i.e. in the part
of the box actually occupied by molecules?

Thank you so very much for your time! 

Andrew
Carnegie Mellon University 

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

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