[gmx-users] NPE simulations: where does all the energy go?

[Vasilii Artyukhov]

Dear colleagues,

There is probably some stupidity involved in my question, but I want to use
GROMACS to explore water phase equilibria using the phase coexistence
technique. This means that my simulation cell contains a liquid slab and an
ice slab in contact, and what I'm trying to do is to see which way the
interface moves depending on the external conditions. The system is
preequilibrated (100-ps NPT run at T slightly above or below the phase
transition temperature).

Now, normally one would want to use NVE simulations to just let the system
reach its equilibrium temperature, however, the density difference between
the liquid and the solid would result in all kinds of issues with pressure.
So it seems that a barostat is required. But what I see in NPE runs is that
my system loses energy. The total energy decreases perfectly linearly, with
both the potential energy and the temperature going down with time. I do
realize that my choice of subject to study might be rather unlucky in this
respect: since liquid water is more dense than ice, my system expands as its
energy decreases (freezing), and it also loses energy as it expands due to
mechanical work. As I understand it, some of this work might become lost
when using the first-order Berendsen barostat; however, what I would expect
from the Parrinello-Rahman barostat is to store this work in the barostat's
own "kinetic energy". And in the PR  case, the temperature drop/total energy
loss is indeed slower, but it's still there. Anyway, the temperature is
supposed to remain constant in this process, but it doesn't.

At this point I should probably mention that I'm using cutoffs for
intermolecular interactions with PME and the long-range dispersion
correction, and not PME-Switch/shifted which are advised for constant-energy
runs; however, I don't really believe that cutoff issues could explain an
effect that strong.

Sorry for the probably cumbersome explanations, but it seems I really need
some help from someone more experienced with NPE simulations to figure
out *where
does all my energy go*. Could there be some algorithmic caveat I'm unaware
of, or perhaps this behavior is completely natural and I'm just missing some
simple basic physics?

Thanks in advance,
Vasilii
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20091018/73a72953/attachment.html>