[gmx-users] Re: v-rescale - harmonic oscillator

[Giovanni Bussi]

Dear Servaas,

the problem that you found is not related to the stochastic rescaling
itself, but to the way the effective energy is calculated in gromacs.
In particular, the increment in the integral part is not consistent
with the applied scaling. You can correct for this by changing
src/mdlib/coupling.c, line 462:

// NEW:
     if (Ek + dEk <= 0) {
        ekind->tcstat[i].lambda = 0.0;
        therm_integral[i] -= -Ek;
      } else {
        ekind->tcstat[i].lambda = sqrt((Ek + dEk)/Ek);
        ekind->tcstat[i].lambda = max(min(1.25,ekind->tcstat[i].lambda),0.8);
        therm_integral[i] -=
Ek*(ekind->tcstat[i].lambda*ekind->tcstat[i].lambda-1.0);
      }
// OLD:
//      if (Ek + dEk <= 0) {
//      ekind->tcstat[i].lambda = 0.0;
//      } else {
//      ekind->tcstat[i].lambda = sqrt((Ek + dEk)/Ek);
//      ekind->tcstat[i].lambda = max(min(1.25,ekind->tcstat[i].lambda),0.8);
//      }
//      therm_integral[i] -= dEk;

Note that this is a problem only when lambda is very different from
one, and this happens in very small systems. The effect for large
systems should be negligible. Furthermore, you should combine this
with the fix that Berk posted a few days ago, which is in file
src/mdlid/update.c, line 165:
// OLD:
//          vv             = lg*(vn + f[n][d]*w_dt);
// NEW:
          vv             = lg*vn + f[n][d]*w_dt;

These fixes will give you a much better energy conservation.

Consider also the following comments:

* For very small systems, the special cases for lambda (max(min...))
introduce small artifacts in the ensemble which are out of control.
Thus, I don't think that you can really use this scheme to perform
rigorous hybrid Monte Carlo. A possible solution is to use a better
integrator for the thermostat (the one discussed in the appendix of
Bussi, Donadio and Parrinello, JCP 2007), which works for every choice
of taut and any number of degrees of freedom, thus it does not need
any special case for lambda. I can send you an implementation of it
for gromacs. Berk is also working on that, and will probabily
introduce it soon in the official gromacs (4.1 or the cvs).

* Even with all these fixes (including the better integrator), there
still a (small) drift in the effective energy for the harmonic
oscillator. I did not test, but I guess that it originates from the
fact that gromacs uses leapfrog and not velocity-Verlet, and leapfrog
is not easily combined with our stochastic scheme which is ultimately
based on Trotter decomposition. The drift should *completely*
disappear for a harmonic system integrated with velocity-Verlet. If
you are interested, we have proven this for a different thermostat
(Bussi Parrinello, PRE 2007, for Langevin dynamics). The demonstration
can be straightforwardly ported to the velocity-rescaling thermostat.

* A part from the two previous implementation issues, I do not expect
any problem with the stochastic rescaling applied to the harmonic
oscillator. In particular, lack of ergodicity in harmonic systems
(which is the big pain in Nose-Hoover) should not be a problem at all.
Furthermore, the algorithm should work and provide the proper ensemble
also on very small systems.

Best regards,

Giovanni Bussi