<html>
<head>
<style>
.hmmessage P
{
margin:0px;
padding:0px
}
body.hmmessage
{
font-size: 10pt;
font-family:Verdana
}
</style>
</head>
<body class='hmmessage'>
Hi,<br><br>I have an optimal leap-frog v-rescale thermostat implemented for CVS head.<br>This will be in Gromacs 4.1.<br>Gromacs 4.1 will probably also have a velocity verlet integrator,<br>for which we can also remove the last very small drift in the v-rescale implementation.<br><br>Berk<br><br>> From: servaas.michielssens@student.kuleuven.be<br>> To: gmx-users@gromacs.org<br>> Date: Tue, 12 May 2009 13:25:02 +0200<br>> Subject: [gmx-users] Re: v-rescale - harmonic oscillator<br>> <br>> Dear Dr. Bussi,<br>> <br>> Thank you very much for the explanation of this problem. <br>> <br>> So I understand that a problem remains with the integrator used in<br>> GROMACS for small systems and that Berk is fixing this (and also a small<br>> error is comming from using leap frog instead of velocity verlet). <br>> <br>> I don't use this thermostat in my hybrid monte carlo algoritm I just<br>> compared the results from a hybrid monte carlo run the results with<br>> v-rescale. So for me there is no hurry the implement the different<br>> integrator, since Berk is already working on it I think it is not very<br>> useful if I also start doing this. I am prepared to help with it or test<br>> it if this is required.<br>> <br>> <br>> Kind regards,<br>> <br>> Servaas<br>> <br>> > -----------------------------<br>> > <br>> > Message: 1<br>> > Date: Tue, 12 May 2009 09:06:52 +0200<br>> > From: Giovanni Bussi <giovanni.bussi@gmail.com><br>> > Subject: [gmx-users] Re: v-rescale - harmonic oscillator<br>> > To: gmx-users@gromacs.org<br>> > Message-ID:<br>> > <3297c3ca0905120006l2f8e630dv34c2ed4b6cd02801@mail.gmail.com><br>> > Content-Type: text/plain; charset=ISO-8859-1<br>> > <br>> > Dear Servaas,<br>> > <br>> > the problem that you found is not related to the stochastic rescaling<br>> > itself, but to the way the effective energy is calculated in gromacs.<br>> > In particular, the increment in the integral part is not consistent<br>> > with the applied scaling. You can correct for this by changing<br>> > src/mdlib/coupling.c, line 462:<br>> > <br>> > // NEW:<br>> > if (Ek + dEk <= 0) {<br>> > ekind->tcstat[i].lambda = 0.0;<br>> > therm_integral[i] -= -Ek;<br>> > } else {<br>> > ekind->tcstat[i].lambda = sqrt((Ek + dEk)/Ek);<br>> > ekind->tcstat[i].lambda = max(min(1.25,ekind->tcstat[i].lambda),0.8);<br>> > therm_integral[i] -=<br>> > Ek*(ekind->tcstat[i].lambda*ekind->tcstat[i].lambda-1.0);<br>> > }<br>> > // OLD:<br>> > // if (Ek + dEk <= 0) {<br>> > // ekind->tcstat[i].lambda = 0.0;<br>> > // } else {<br>> > // ekind->tcstat[i].lambda = sqrt((Ek + dEk)/Ek);<br>> > // ekind->tcstat[i].lambda = max(min(1.25,ekind->tcstat[i].lambda),0.8);<br>> > // }<br>> > // therm_integral[i] -= dEk;<br>> > <br>> > Note that this is a problem only when lambda is very different from<br>> > one, and this happens in very small systems. The effect for large<br>> > systems should be negligible. Furthermore, you should combine this<br>> > with the fix that Berk posted a few days ago, which is in file<br>> > src/mdlid/update.c, line 165:<br>> > // OLD:<br>> > // vv = lg*(vn + f[n][d]*w_dt);<br>> > // NEW:<br>> > vv = lg*vn + f[n][d]*w_dt;<br>> > <br>> > These fixes will give you a much better energy conservation.<br>> > <br>> > Consider also the following comments:<br>> > <br>> > * For very small systems, the special cases for lambda (max(min...))<br>> > introduce small artifacts in the ensemble which are out of control.<br>> > Thus, I don't think that you can really use this scheme to perform<br>> > rigorous hybrid Monte Carlo. A possible solution is to use a better<br>> > integrator for the thermostat (the one discussed in the appendix of<br>> > Bussi, Donadio and Parrinello, JCP 2007), which works for every choice<br>> > of taut and any number of degrees of freedom, thus it does not need<br>> > any special case for lambda. I can send you an implementation of it<br>> > for gromacs. Berk is also working on that, and will probabily<br>> > introduce it soon in the official gromacs (4.1 or the cvs).<br>> > <br>> > * Even with all these fixes (including the better integrator), there<br>> > still a (small) drift in the effective energy for the harmonic<br>> > oscillator. I did not test, but I guess that it originates from the<br>> > fact that gromacs uses leapfrog and not velocity-Verlet, and leapfrog<br>> > is not easily combined with our stochastic scheme which is ultimately<br>> > based on Trotter decomposition. The drift should *completely*<br>> > disappear for a harmonic system integrated with velocity-Verlet. If<br>> > you are interested, we have proven this for a different thermostat<br>> > (Bussi Parrinello, PRE 2007, for Langevin dynamics). The demonstration<br>> > can be straightforwardly ported to the velocity-rescaling thermostat.<br>> > <br>> > * A part from the two previous implementation issues, I do not expect<br>> > any problem with the stochastic rescaling applied to the harmonic<br>> > oscillator. In particular, lack of ergodicity in harmonic systems<br>> > (which is the big pain in Nose-Hoover) should not be a problem at all.<br>> > Furthermore, the algorithm should work and provide the proper ensemble<br>> > also on very small systems.<br>> > <br>> > Best regards,<br>> > <br>> > Giovanni Bussi<br>> > <br>> > <br>> > ------------------------------<br>> <br>> _______________________________________________<br>> gmx-users mailing list gmx-users@gromacs.org<br>> http://www.gromacs.org/mailman/listinfo/gmx-users<br>> Please search the archive at http://www.gromacs.org/search before posting!<br>> Please don't post (un)subscribe requests to the list. Use the <br>> www interface or send it to gmx-users-request@gromacs.org.<br>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php<br><br /><hr />See all the ways you can stay connected <a href='http://www.microsoft.com/windows/windowslive/default.aspx' target='_new'>to friends and family</a></body>
</html>