fix_nh possible missing factor in nve_x

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

f1 = exp(2*arg),
f2 = exp(arg) * sinh(arg)/arg
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.


Aidan can answer this specifically. But nve_x() is
only part of the update within FixNH.


The factor f2 is applied directly to the velocities in function
nh_v_press(). We truncated the Maclaurin series at the first term:
sinh(x)/x = 1 + O(x^2)..