Re: [ESPResSo-devel] Thoughts on implementing an improved thermostat for ESPResSo

[Ulf Schiller]

Evyatar,

What your are suggesting in 1)-5) is actually not so much related to the
thermostat but in fact a modified integration scheme. The ESPResSo
integrator is a Verlet algorithm with two explicit updates of the
velocity. Hence there is no need to keep the forces from the previous
configuration. The important thing to note is that, in case of velocity
dependent forces, this integrator is only first order, see
http://www2.mpip-mainz.mpg.de/~andrienk/journal_club/integrators.pdf.
For Langevin dynamics, this can be easily turned into a second order
scheme by noting that the stochastic momentum update is equivalent to an
Ornstein-Uhlenbeck process. In my opinion, the order of the integrator
is much more important for stability than the thermostat (which
typically has a stabilising effect anyway).

I think the improvement reported in the paper you refer to is based on
an atomistic model and the comparison is for specific algorithms
implemented in GROMACS, LAMMPS, and AMBER.

If you want to try and implement an updated scheme in ESPResSo, keep in
mind that there are a number of other algorithms involved, such as an
NPT integrator and fluid-particle coupling, which all have to be
addressed consistently. I think the best way to go about it is to get a
thorough understanding of the ESPResSo integrator first. The PDF linked
above may help you to sort the various steps out.

Hope this helps,
Ulf

On 02/06/15 11:23, Evyatar Arad wrote:
> Hello,
> 
> My name is Evyatar Arad, and as part of my studies Im working on an
> implementation of a newly released thermostat, G-JF, in the ESPResSo
> package.
> The new thermostat exhibits stable results with use of much longer time
> steps (compared to existing options, as shown in attached article) and
> thus may save the ESPResSo users community a-lot of computational time
> if integrated properly (see attached article, p3).
> I’m rather new to the ESPResSo c++ code and package and might need some
> help regarding the best way to properly integrate it.
> I’ve taken a look around and saw the 4 phases of Velocity-Verlet
> integration in the ESPResSo code.
> 
> The way I see the G-JF integration loop should go something like:
> 
> 1) Randomise gaussian numbers and from them derive Beta^(n+1) for every
> particle.
> 2) propagate pos by calculating eq 4 in attached article
> 3) recalc forces to find f^(n+1) (while keeping f^n)
> 4) propagate velocities by calculating eq. 5 in attached article
> 5) f^(n+1) => f^n
> 
> The problem I’ve been having is of an orderly nature. I’m just not sure
> where to place each step in the tangled code and was wondering how do
> you think is best to go about it?
> 
> I’ve enclosed a 4 pages article for reference.
> Application of the G-JF discrete-time thermostat for fast and accurate
> molecular simulations:
> http://www.bgu.ac.il/~ofarago/cpc-gjf.pdf
> 
> Thank you.
> 
> Evyatar Arad,
> Department Of Biomedical Engineering,
> Ben Gurion University.
> 
> 


-- 
Dr Ulf D Schiller
Centre for Computational Science
University College London
20 Gordon Street
London WC1H 0AJ
United Kingdom

Phone: +44 (0)20 7679 5300