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

[Berk Hess]

Hi,

I have an optimal leap-frog v-rescale thermostat implemented for CVS head.
This will be in Gromacs 4.1.
Gromacs 4.1 will probably also have a velocity verlet integrator,
for which we can also remove the last very small drift in the v-rescale implementation.

Berk

> From: servaas.michielssens at student.kuleuven.be
> To: gmx-users at gromacs.org
> Date: Tue, 12 May 2009 13:25:02 +0200
> Subject: [gmx-users] Re: v-rescale - harmonic oscillator
> 
> 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
> > 
> > 
> > ------------------------------
> 
> _______________________________________________
> gmx-users mailing list    gmx-users at 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-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php

_________________________________________________________________
See all the ways you can stay connected to friends and family
http://www.microsoft.com/windows/windowslive/default.aspx
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20090512/a89be139/attachment.html>