Dear Users and Developers,

in all likelihood I have just overseen something, but as I'm not sure I post this question.

In the member function FixNH::nve_x() atom positions are updated like:

x[i][0] += dtv * v[i][0];

But according to Molecular Physics 87.5 (1996): 1117-1157 the operator acting on the positions should be

exp( dt * \sum_{k=1}^{Natom} [ v_k + v_{\epsilon}*r_k] \nabla_{r_k}

which translates to (omitting _k)

r -> r*f1 + v*f2*dt

where

f1 = exp(2*arg),

f2 = exp(arg) * sinh(arg)/arg

and

arg = dt/2*v_{\epsilon},

and v is the velocity after the first velocity half step.

v_{\epsilon} matches omega_dot in the code if I'm correct.

f1 seems to be applied in member function FixNH::remap(), but I couldn't see where f2 is applied unless it's somehow hidden in dtv, which I could'nt retrace.

I'm not very familar with lammps code base, so I did a grep for one of the maclaurin series prefactors (1/362880) of sinh(x)/x, the only file showed up was fix_rigid_nh.h and .cpp where f2 is applied.

Perhaps someone can point out what I'm missing.

For the record, I'm referring to code of the latest stable tarball.

Best,

Richard