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