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

[servaas]

Dear Dr. Bussi,

Thank you very much for the explanation of this problem. 

So I understand that a problem remains with the integrator used in
GROMACS for small systems and that Berk is fixing this (and also a small
error is comming from using leap frog instead of velocity verlet). 

I don't use this thermostat in my hybrid monte carlo algoritm I just
compared the results from a hybrid monte carlo run the results with
v-rescale. So for me there is no hurry the implement the different
integrator, since Berk is already working on it I think it is not very
useful if I also start doing this. I am prepared to help with it or test
it if this is required.


Kind regards,

Servaas

> -----------------------------
> 
> Message: 1
> Date: Tue, 12 May 2009 09:06:52 +0200
> From: Giovanni Bussi <giovanni.bussi at gmail.com>
> Subject: [gmx-users] Re: v-rescale - harmonic oscillator
> To: gmx-users at gromacs.org
> Message-ID:
> 	<3297c3ca0905120006l2f8e630dv34c2ed4b6cd02801 at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
> 
> 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
> 
> 
> ------------------------------