Dear colleagues,<br><br>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&#39;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).<br>

<br>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&#39;s own &quot;kinetic energy&quot;. And in the PR  case, the temperature drop/total energy loss is indeed slower, but it&#39;s still there. Anyway, the temperature is
supposed to remain constant in this process, but it doesn&#39;t.<br><br>At this point I should probably mention that I&#39;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&#39;t really believe that cutoff issues could explain an effect that strong.<br>

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

<br>Thanks in advance,<br>Vasilii<br>