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