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.bu...@gmail.com> > Subject: [gmx-users] Re: v-rescale - harmonic oscillator > To: gmx-users@gromacs.org > Message-ID: > <3297c3ca0905120006l2f8e630dv34c2ed4b6cd02...@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 > > > ------------------------------ _______________________________________________ gmx-users mailing list gmx-users@gromacs.org http://www.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php