How to integrate numerically Nosé Hoover equation?

Question

In NVT molecular dynamics, Nosé Hoover thermostat is a method defining an extended system. I can understand perfectly how the Nosé Hoover differential equations are derived (Frenkel&Smit's book). However, when it comes to practice, I have no idea how to solve it numerically. I tried a method similar to velocity Verlet to integrate the full equations, but it does not really work. What is mysterious is that all the textbooks I have seen so far all stop at the differential equations and do not explain how to numerically integrate it, and this is the hard part. Could anyone suggest any references that clearly explain how to discretize Nosé Hoover equations, or any written codes to look at?

Thank you!

Answer 1

You can find explicit pseudo-code for Nosé-Hoover Chains in Section 4 of:

Martyna, G. J., Tuckerman, M. E., Tobias, D. J., & Klein, M. L. (1996). Explicit reversible integrators for extended systems dynamics. Molecular Physics, 87(5), 1117–1157. Retrieved from http://www.tandfonline.com/doi/abs/10.1080/00268979600100761

Here is a short example of the Nosé-Hoover-Langevin thermostat on a 1D harmonic oscillator implemented in MATLAB ($q$ is position, $p$ is momentum, $\xi$ is the degree of freedom associated to the thermostat, $\gamma$ is the friction, $kT$ is the temperature times the Boltzmann constant, and $dt$ is the time-step length):

function [q p xi] = nhl(q, p, xi, gamma, kT, dt)

mass = 1.0;

kinetic = @(p) 1/(2*mass) * p^2;

force = @(q) -q;

N = 1;                                  % Number of dimensions

xi = exp(-gamma*dt/2) * xi + sqrt(kT/mass*(1 - exp(-gamma*dt))) * randn;

p = p + dt/2 * force(q);

p = p * exp(-xi*dt/4);
xi = xi + dt/2 * (2*kinetic(p) - N*kT)/mass;
p = p * exp(-xi*dt/4);

q = q + dt/mass * p;

p = p * exp(-xi*dt/4);
xi = xi + dt/2 * (2*kinetic(p) - N*kT)/mass;
p = p * exp(-xi*dt/4);

p = p + dt/2 * force(q);

xi = exp(-gamma*dt/2) * xi + sqrt(kT/mass*(1 - exp(-gamma*dt))) * randn;

You can find the general Nosé-Hoover-Langevin method described explicitly in:

Leimkuhler, B., Noorizadeh, E., & Theil, F. (2009). A Gentle Stochastic Thermostat for Molecular Dynamics. Journal of Statistical Physics, 135(2), 261–277. http://doi.org/10.1007/s10955-009-9734-0

Since these two methods (NHC and NHL) are closely related to the original NH method, understanding how to get from the ODEs to the actual integrator in any of these methods will lead you to understand how it's done for NH.

Answer 2

tl;dr

The Nose-Hoover equations are normally defined by their Hamiltonian. The two partial derivatives of the Hamiltonian define a partitioned set of ODEs. From there the partitioned ODEs are either solved by an integrator for 1st order ODEs or some symplectic method.

Details shown using DifferentialEquations.jl

As a rather thorough example, let's look at the architecture for how DifferentialEquations.jl handles it. It allows you to define a dynamical ODE problem directly using the Hamiltonian. What this is really doing is defining the partitioned ODE version using the partial derivatives of the Hamiltonian, i.e. the equations of motion for the position and momentum given by

$$\dot{p} = -\frac{\partial H}{\partial q}$$ $$\dot{q} = \frac{\partial H}{\partial p}$$

from Hamilton's relations. For systems of equations this is the gradient, so you can see from the source code how this is done using autodifferentiation and define it similarly in any language.

Once you have the partitioned set of ODEs, then there are many things you can do. If you group them together, you can call this one giant system of first order ODEs. In that case, you can use any standard ODE solver on this system (which for DifferentialEquations.jl it's this list, or in MATLAB that would be something like ode45, or in SciPy you could use LSODA, etc.).

But you "know more structure" about the equation, and so theoretically methods could be made to use this structure and be more efficient, right? Yes, that is the case. There's a big list of them here on the DifferentialEquations.jl page, but we can split them into two classes: Runge-Kutta Nystrom methods and Symplectic integrators. Runge-Kutta Nystrom methods are extensions of Runge-Kutta methods specifically developed for second order ODEs, and this structure you have here (a position equation and a momentum equation) is this second order ODE structure. These are able to be more efficient than the first order methods on this specific class of problems by cutting out some function evaluations.

Symplectic integrators are more interesting. These integrators satisfy the property that their integration conserves the symplectic two-form. This commonly makes people believe that they conserve energy, that is not true. However, for long term integration they will stay on the symplectic manifold and thus not drift away, leading to better estimates of things like energy and angular momentum in the long run. However, these methods have a few extra costs associated with them (of course, because they have to satisfy more properties!). One big thing is that they are fixed timestep (there are adaptive timestep variants, but the recent research in the topic doesn't show that they are that efficient, mostly because there's a lot of extra stuff that needs to go on for adaptive timestep to stay symplectic).

Recap

So let's recap a little bit. You either analytically, symbolically, numerically, or via autodifferentiation calculate the partitioned ODE for $(\dot{p},\dot{q})$: which methods do you choose? If your problem isn't too hard, just throw that into your standard ODE solver and you're good. If you're doing long time integration and need to preserve some extra properties like energy and angular momentum "better", use some symplectic integrator (here's a whole list, a popular one for this is Velocity Verlets because it's easy to implement and reduces force calculations, which isn't actually super efficient but I'll get to that some other time...). If need it efficient and aren't doing long time integration, one of the Runge-Kutta Nystrom methods is what you're looking for.

Available software

So DifferentialEquations.jl has it all of course as shown above. But if you're not using Julia, the papers for each of these methods can be found on the citations page. If you need a symplectic integrator and need to roll it yourself, Velocity Verlet is a good thing to try first. Other than that, you can find a few of these methods in Fortran, the Dormand-Prince Nystrom method in MATLAB, a few symplectic integrators in Python which wrap REBOUND, and a symplectic integrator in Mathematica, but I don't know of another large set of methods which has it all together to benchmark with so Julia is the choice if you want/need to do that.

Concluding benchmarks and examples

I'll end by leading to a few benchmarks and examples which could help you choose the right method. These are taken from DiffEqTutorials.jl and DiffEqBenchmarks.jl, but let me point you to ones which are relevant:

1. This is a quick demonstration of the difference in energy error between a symplectic integrator and a standard 1st order integrator on a few classic problems.
2. This problem shows the "drift" that occurs in non-symplectic integrators over long time integrations in the Kepler problem
3. This notebook benchmarks high-precision solutions using various methods when high energy accuracy is needed
4. This is another high accuracy benchmark.

While none of those are the Nose-Hoover equation itself, those are all benchmarks of other problems defined directly from the Hamiltonian and display the differences between the types of methods, so that should help you to extrapolate from there to see what you need.

Ending quick recommendation

If I had to say one recommendation without knowing exactly what you're doing, and assuming you have to implement it yourself (otherwise... try them all and benchmark!), I'd say go with the Velocity Verlet (on the partitioned set of ODEs defined from the Hamiltonian). It's pretty standard here. Given the benchmarks I have done, symplectic is usually unnecessary and overused (people overplay what symplectic actually means, and in this class Verlets aren't even the most efficient), but they are easy to implement, standard, and give you safety for long-time integrations.

Edit

Oops, I mixed it up with the Nose system which is a Hamiltonian system. The Nose-Hoover system is non-Hamiltonian, so the conversion to 1st order ODEs for standard integrators works and so does a form of the Verlet method, but standard symplectic integrators do not apply. We'll get some implementations to benchmark this field "soon enough".

Answer 3

The simplest way is to use fourth-order Runge-Kutta integration.