Shouldn't econserve stay constant in NpT, and NpH?

I’m looking for a stable timestep for a MEAM NiTi potential. Usually, I’m looking at stability of temperature in long (~15ns, with dt=0.001ps it is 15 millions of steps) NVE, and NpH runs. And at stability of etotal = total energy (pe + ke) in long NpT, and NpH runs.

Then I stumbled upon econserve = pe + ke + ecouple = etotal + ecouple. I think, it should stay constant during NpT, and NpH runs, shouldn’t it?

In reality, I can have relatively stable etotal for small timestep, but econserve still drifts to large positive or negative values. It looks like for a stable timestep it drifts to large negative values and for unstable large timestep it drifts to positive values. Can someone explain this? Is econserve not a good quantity to look at, when testing the stability of the integration?

1 Like

This is a question for @athomps

Formally, because of the symplectic integrators that are coded in LAMMPS, Econserve for Nose-Hoover-style dynamics should exhibit the same characteristics as Etotal in NVE dynamics i.e. for sufficiently small timestep, fluctuations about an unknown but fixed shadow value should be of bounded magnitude, independent of simulation duration, and the fluctuation magnitude should scale as O(dt^2). In practice, for a variety of reasons, this behavior may not be observed over long simulations. For example, if the MEAM potential infrequently visits regions of high stiffness, this may violate the “sufficiently small timestep” requirement.

From a practical point of view, you should convince yourself that whatever science result you are getting from your simulation is insensitive to the choice of timestep size. If not, then you will either have to decrease your timestep size until that condition is satisfied, or you may have to completely rethink your simulation protocol.

1 Like