Dear Steve,

I think I found a bug in velocity Verlet integrator in fix_nvt.cpp file.

In that file, the velocity is updated as

void FixNVT::final_integrate()

{

:

for (int i = 0; i < nlocal; i++) {

if (mask[i] & groupbit) {

dtfm = dtf / mass[type[i]] * factor;

v[i][0] = v[i][0]*factor + dtfm*f[i][0];

v[i][1] = v[i][1]*factor + dtfm*f[i][1];

v[i][2] = v[i][2]*factor + dtfm*f[i][2];

:

}

Instead, it should be updated as

v[i][0] = (v[i][0] + dtfm*f[i][0])*factor;

v[i][1] = (v[i][1] + dtfm*f[i][1])*factor;

v[i][2] = (v[i][2] + dtfm*f[i][2])*factor;

I checked the reversibility and also found the same pattern in

fix_nph.cpp and fix_npt.cpp files.

Happy New Year,

Keonwook