# Hamiltonian for a restricted isoenergetic thermostat

## Abstract

Nonequilibrium molecular dynamics simulations often use mechanisms called thermostats to regulate the temperature. A Hamiltonian is presented for the case of the isoenergetic (constant internal energy) thermostat corresponding to a tunable isokinetic (constant kinetic energy) thermostat, for which a Hamiltonian has recently been given.

Thermostats are modifications to the equations of motion of a classical system to simulate thermal interaction of a system with the environment. The Nosé-Hoover thermostat is used to simulate fluctuations in energy of an equilibrium system corresponding to the canonical ensemble of statistical mechanics, and the Nosé-Hoover and Gaussian thermostats, among others, are used to remove heat from a system driven by external forces into a nonequilibrium stationary state . There has been a recent interest in thermostatted equations of motion, focussed on the symplectic structure of the equations of motion, and the related pairing of the Lyapunov exponents. Both a Hamiltonian and pairing of Lyapunov exponents are known for Nosé-Hoover and Gaussian isokinetic (GIK: constant kinetic energy) thermostats . Numerical evidence against pairing (and hence the existence of a Hamiltonian) are discussed in  for the GIK thermostat applied to shearing systems and in  for the Gaussian isoenergetic (GIE: constant internal energy) thermostat. The latter paper does, however show that in a special case of the GIE thermostat, involving one rather than two arbitrary potentials, the Lyapunov exponents are paired. The purpose of this short note is present a Hamiltonian for this case.

The GIE thermostat has equations of motion of the form

$`{\displaystyle \frac{d𝐱_i}{dt}}={\displaystyle \frac{𝐩_i}{m_i}},`$ $`{\displaystyle \frac{d𝐩_i}{dt}}={\displaystyle \frac{\mathrm{\Phi }^{(ext)}}{𝐱_i}}{\displaystyle \frac{\mathrm{\Phi }^{(int)}}{𝐱_i}}\alpha 𝐩_i,`$ (1)   
$`\alpha `$ $`=`$ $`{\displaystyle \frac{_i\frac{𝐩_i}{m_i}\frac{\mathrm{\Phi }^{(ext)}}{𝐱_i}}{_i𝐩_i𝐩_i/m_i}},`$ (2)   
where $`\mathrm{\Phi }^{(ext)}`$ is the external driving potential, $`\mathrm{\Phi }^{(int)}`$ the interparticle potentials, and $`\alpha `$ is the thermostat term which ensures that the equations conserve internal energy $`E=_i𝐩_i^2/(2m_i)+\mathrm{\Phi }^{(int)}`$. The equations reduce to no thermostat when $`\mathrm{\Phi }^{(ext)}=0`$ and to GIK when $`\mathrm{\Phi }^{(int)}=0`$. A more general example of a limit involving only one arbitrary potential is the case $`\mathrm{\Phi }^{(ext)}=\gamma \mathrm{\Phi }`$, $`\mathrm{\Phi }^{(int)}=(1\gamma )\mathrm{\Phi }`$, leading to the equations

$$\frac{d𝐱_i}{dt}=\frac{𝐩_i}{m_i},\frac{d𝐩_i}{dt}=\frac{\mathrm{\Phi }}{𝐱_i}+\gamma \frac{_i\frac{𝐩_i}{m_i}\frac{\mathrm{\Phi }}{𝐱_i}}{_i𝐩_i𝐩_i/m_i}𝐩_i,$$
(3)   
which conserve energy $`E=_i𝐩_i^2/(2m_i)+(1\gamma )\mathrm{\Phi }`$. Here, $`\gamma `$ effectively controls the strength of the thermostat from no thermostat ($`\gamma =0`$), to the GIK thermostat ($`\gamma =1`$). For any $`\gamma `$ the Lyapunov exponents are paired , suggesting the existence of a Hamiltonian.

Following the GIK case , the conservation law is enforced by setting the numerical value of the Hamiltonian equal to the conserved energy, assigned the value zero by a shift in the potential energy. This allows the kinetic energy term in the denominator of (2) to be replaced by minus the potential energy (note $`\mathrm{\Phi }<0`$),

$$\alpha =\frac{\gamma }{2(1\gamma )}\underset{i}{}\frac{𝐩_i}{m_i}\frac{}{𝐱_i}\mathrm{ln}|\mathrm{\Phi }|.$$
(4)   
Another aspect of a Hamiltonian description of thermostatted systems is that in the physical variables $`(𝐱,𝐩)`$ there is a phase space volume contraction rate proportional to $`\alpha `$, while in the canonical variables $`(𝐱,𝝅)`$ phase space volume is conserved. This means that $`𝝅`$ must be greater than $`𝐩`$ by a factor equal to $`\mathrm{exp}\left(\alpha 𝑑t\right)=|\mathrm{\Phi }|^{\gamma /(2(1\gamma ))}`$. Multiplying the zero energy by an arbitrary power of $`|\mathrm{\Phi }|`$ we have

$$H_\beta (𝐱,𝝅,\lambda )=|\mathrm{\Phi }|^{\frac{\gamma }{1\gamma }+\beta }\underset{i}{}\frac{𝝅_i^2}{2m_i}+(1\gamma )\mathrm{\Phi }|\mathrm{\Phi }|^\beta ,$$
(5)   
which, combined with the constraint $`H_\beta =0`$ and the identifications $`dt=|\mathrm{\Phi }|^{\gamma /(2(1\gamma ))+\beta }d\lambda `$ and $`𝐩_i=|\mathrm{\Phi }|^{\gamma /(2(1\gamma ))}𝝅_i`$ leads to the equations of motion, (3). Interesting cases are $`\beta =\gamma /(2(1\gamma ))`$ for which there is no time scaling, $`\beta =0`$ has a certain simplicity, $`\beta =\gamma /(1\gamma )`$ yields the familiar form of kinetic plus potential energy, and $`\beta =1`$ for which the Hamiltonian is that of a geodesic in a conformally flat space, see .

The author is grateful for discussions with W. G. Hoover.