Barostat issues with two-phase system

Hello all,

I am having an issue with a coexistence study involving a gas phase and a solution phase. My box is a rectangular prism with the gas phase as a “slab”. The key issue is that the pressure tensor elements are not equal on the diagonal. With only fluid phases, one would expect that Pset =

= = = , where Pset is what pressure we set our barostat to and the bracketed Ps correspond to the average pressures (overall and in the x, y, z dimensions, respectively). However, the results that I am seeing indicate that

is sometimes close to Pset, but that is often greater than Pset (though not always) and that = < Pset. The x dimension is longer than the other two dimensions (by ~ a factor of 3), as well as being the direction which is normal to the gas/solution interface. I have tried several different barostatting scenarios; in general all give the aforementioned pressure relationships. These include:

fix 1 all npt temp $T T {Tdamp} x $p p {Pdamp}
(only the x box dimension is changed)
fix 1 all npt temp $T T {Tdamp} iso $p p {Pdamp}
(x, y, and z box dimensions are changed, each by same percentage)

Finally, the only simulation that I have gotten to successfully maintain the pressure elements as approximately equal is one in which I used anisotropic pressure control:

fix 1 all npt temp $T T {Tdamp} aniso $p p {Pdamp} couple yz
(x, y, and z dimensions are all changed, and can change independently)

However, the anisotropic case has the unfavorable characteristic that my box aspect ratios get severely distorted (x dimension becomes very long at the expense of y and z, probably to reduce the interfacial area between the two phases).

Probably the strangest part is that this inequality of pressure elements occurs even for NVT conditions

fix 1 all nvt temp $T T {Tdamp}

I am trying out some other scenarios right now, but these results seem odd. Do you have any idea why they might be occurring and how they might be fixed?

Thank you very much, I greatly appreciate any insight,

p.s. I’ve attached a sample “in” file with forcefield file and restart file

in.coex (891 Bytes)

MeWa.sw (800 Bytes)

prod2.restart (688 KB)

Maybe Aidan wants to comment on this. It's not
clear to me that the 3 diagonal components should
be equal for this system, so there may be nothing
you can (or should) do about it.



For a system with multiple interfaces, you should get exactly the behavior that you describe below.

The assumption of isotropic pressures is correct only when the system consists of a homogeneous phase. A two-phase solution, even if both phases are fluids, will have an active interface, and therefore a surface tension. This of course requires that the pressure tensor components not be equal.

You can see Israelachvili’s book or any number of papers in the literature on this point.


One additional point: you do realize that you
can use fix npt to barostat just a single dimension
of the box, against a single component of
the pressure tensor. E.g. just barostat the dimension
normal to the interface if you like.



We call the condition of Pxx=Pyy=Pzz hydrostaticity. It only occurs under
certain conditions. Those include homogeneous fluids. You fluid is
inhomogeneous, with two interfaces oriented normal to the x direction.

"The x dimension is longer than the other two dimensions (by ~ a factor of
3), as
well as being the direction which is normal to the gas/solution interface."

Hence it is natural to expect that Pxx != Pyy = Pzz. Look up any text book
on MD simulations of liquid interfaces for an explanation of why.


Thank you Aidan, Steve and Ahmed - your responses are all very helpful. Obviously, I was caught up on fluids having isotropic pressure tensors and neglecting the surface tension between two fluids. I’m going to go read my introductory textbooks now.

Thanks again,