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