# How to integrate numerically Nosé Hoover equation?

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!

• In the current state, I don't think that your question is going to have a lot of attention. Some tips: provide the full bibliographic information of the book, describe the method for those that don't know it, try to describe better what you have done and the errors that you get (consider adding some snippets of code). I found this code when looking in GitHub. It also seems to be implemented in LAMMPS. Oct 14 '16 at 15:08

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.

## 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:

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".

• The feature set of DifferentialEquations.jl is very impressive. Despite the fact that there is such a thing as the Nose Hamiltonian, the Nose-Hoover thermostat is not a Hamiltonian system. Oct 30 '17 at 13:54
• Yes it is a Hamiltonian system, if the one we are discussing is the one that is mentioned here, here, here, etc. If most of those sources are correct, then it's different from the Langivan-type systems and this one is defined by the Hamiltonian given in those sources. Oct 30 '17 at 14:11
• However, I could be wrong: my research is in the methods and I know about this problem as "something these methods are applied to", so please let me know if there's a difference between the Nose-Hoover mentioned in the OP and the Nose-Hoover Hamiltonian systems mentioned in those systems. If it is a Langivan equation with additive noise, then I'll edit the post to be about solvers for additive noise SDEs, but I always heard that under the name Brownian dynamics. Oct 30 '17 at 14:13
• I use DifferentialEquations.jl as an example because I am not sure where else you can find all of those methods together to actually disucss the difference between them, other than a few scattered around Hairer's books (specifically, there are a few Nystrom methods in Hairer I and the symplectic integrators are Hairer III (Geometric Methods), but most of the important ones like the adaptive Dormand-Prince Nystrom methods are left out). I hope to collect this info soon, but for now the docs+source is the most complete discussion that I know about. Oct 30 '17 at 14:16
• Hi Chris. The Nose-Hoover thermostat is deterministic (although there are variants such as Nose-Hoover-Langevin that are not), as you say. However, the corresponding dynamical system formulation of N-H does not follow from the Nose Hamiltonian that you point to. As a matter of fact, in the second reference that you provide in the comments there are two different sections: one for the Nose thermostat (which is a Hamiltonian system) and another section for the Nose-Hoover thermostat (which is not a Hamiltonian system). Oct 30 '17 at 14:30

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