Integration of the barostat equations of motion in fix_nh

Hello everyone,

I'm investigating some differences between trajectories generated by LAMMPS
using fix/npt and our own MD code. I have already identified two sources of
differences and proposed a fix in PR#942. However, I'm having a hard time
understanding the time integration scheme used for barostatting in LAMMPS and
I'd be very glad if anyone could shed some light into this.

The problem is, there are relatively large deviations in box shapes/
trajectories between our code and LAMMPS when fix/npt is used in the fully
triclinic (shape-changing) mode. The box parameters already exhibit small but
visible differences after a single integration step and are completely
different after ~100 steps. If I decouple the tilt factors from the barostat
(forcing them to stay at their original values), the trajectories stay
indistinguishable for much longer (~500 steps), so the original differences
are not caused just by numerical noise.

Looks like the evolution of omega_dot is almost exactly the same as in our
code, but the resulting coordinate/velocity scaling differs. I have noticed
that the way the coordinates and velocities are scaled/integrated by fix_nh is
completely different from the algorithm in the cited Tuckerman2006 paper or in
Yu et al. 2010 (doi:10.1016/j.chemphys.2010.02.014). Most importantly, fix_nh
doesn't do the eigenvector transformation to decouple the individual
components prior to integration. The hyperbolic sine factors are also not
present at all.

Is this a deliberate approximation? If so, is the integration scheme used in
LAMMPS published and analyzed in a paper somewhere?

Best regards,

Tomáš

Aidan (cc’d) would be the most competent LAMMPS developer to comment on this subject.

axel.

yes. I’m sure Aidan will look into this and respond.

However he is on extended travel at the moment.

Thanks for the careful explanation of your Qs.

Steve