NPT, damping and pressure

Hi all,

after extensive testing, I suspect there may be one more off-by-one error somewhere in the fix_nh implementation. But I'm not 100% sure, which is why I want to discuss here before opening a bug report.

Under long enough ensemble averages, fix npt should reach a state corresponding to the prescribed T and p conditions. Large systems will require less time averaging (relatively smaller fluctuations), small systems will require longer time averages (more fluctuations). The intensive results should be equivalent if converged and no other errors were made.

As others (such as Steven Vandenbrande in [1]) have discovered, the LAMMPS implementation of the Nose-Hoover thermostat seems to be very sensitive to the choice of pressure relaxation time \tau_p. On the side of large \tau_p, it may take "forever" (longer than any reasonable simulation time) to reach pressure equilibrium, on the other side, there will be catastrophic oscillations if \tau_p < 5\tau_T for many metals because of their high bulk modulus. This can be somewhat alleviated by adding damping using the "drag" option. I am aware that this changes the velocity *distributions*, but it shouldn't affect the long-term averages of quantities.

Now, the issue: if drag>~0.5 and mtk=yes, the average of the pressure reported by thermo press will always be too low by *exactly* kBT/vol, as in:
     <press> = p_set - k_B <T> / <vol>

As far as I can see, all other quantities are correct, especially <temp>, <etotal>, <enthalpy> and <vol> are unaffected by the setting of drag. This leads me to believe that only the reported value of press is wrong, but the internal setpoint and dynamics are in fact correct?

Setting mtk=no removes this effect, which leads me to believe something in the MTK correction terms is missing in a form that only shows up with damped cell dynamics?

This is consistently reproducible using this [2] sample input (not quite minimal, sorry -- directly derieved from a real problem), for various system sizes, temperatures and pressures. Note it is a binary if: there is no scaling relationship to the value of drag, it's a hard switch from "permanently oscillating and 'correct'" to "converging and constant offset". ploop and pchain do not affect the outcome, I used *chain=1 because plots are easier to read.

LAMMPS version is current unstable, but I also tested versions before last fall's fixes to fix_nh.cpp and they show the same behaviour (so, not a regression).

May I ask if someone has encountered this before, this is expected or... something?

Regards,

Sebastian Hütter

[1] <https://lammps.sandia.gov/threads/msg74862.html>
[2] <https://puu.sh/CDi8V/2c60ec0101.zip/testnpt.zip>

Hi Sebastian,

It is not clear that this is a problem with the code:

1. The two issues identified by Steven Vandenbrande that you cite were confirmed by him (on the same thread) to be resolved by this pull request:

https://github.com/lammps/lammps/pull/1259

2. I would not trust any results generated using "drag" option. It is a useful but completely ad hoc feature.

Aidan

Hi Aidan,

1. The two issues identified by Steven Vandenbrande that you cite were confirmed by him (on the same thread) to be resolved by this pull request:

https://github.com/lammps/lammps/pull/1259

And I even commented on it, that's embarrassing.

2. I would not trust any results generated using "drag" option. It is a useful but completely ad hoc feature.

I know. And for being that, it is better than one might suspect: <http://puu.sh/CE85a/291848b5d0.png> Tracks almost exactly the same pressures near the beginning, and then does not start to oscillate.

And I must correct myself: mtk=no does not really fix the output - the difference between these two just happens to be the exact same value as well, so it cancels out. For example:

mtk drag <p> <H> <vol>
1 0 1 -17276.585 44474.059
1 5 -0.5 -17276.672 44473.822
0 0 2.5 -17276.653 44473.988
0 5 1 -17276.571 44473.944

This seems odd. Either I'm missing something fundamental, or the computed pressure value does not always make sense? If the resulting volume is the same, pressure can't be different unless the pair forces are different, right?

Regards,

Sebastian Hütter