Separating Transformations for Controlled-Variable Molecular Dynamics

Ben Leimkuhler Department of Mathematics 405 Snow Hall University of Kansas Lawrence, KS 66045 USA leimkuhl@math.ukans.edu

Abstract

We show that Nosé’s thermostat can be given a mechanical formulation (with Hamiltonian having a Euclidean kinetic energy metric). We also give a mechanical formulation for simultaneous control of temperature and pressure. This approach simplifies the construction of symplectic methods while providing a more intuitive perspective on the nature of the controlled variable molecular dynamics.

PACS Codes: 02.70.Ns, 05.20.-y, 82.20.Fd

Keywords: Andersen constant pressure method, Nosé thermostat, molecular dynamics, Hamiltonian systems, symplectic integrators

I. INTRODUCTION

Extended Hamiltonian formulations for controlling variation of some physical quantity during simulation originated with Andersen’s 1979 article on constant pressure simulation.

[1] Nosé later formulated a scheme for temperature control, [2] and these techniques are now widely used to perform classical molecular simulations in the canonical ensemble. [3–12] The mixing of positions and momenta in the kinetic energy complicates dynamical analysis of the Nosé and Andersen formulations and impedes the construction of efficient symplectic integration methods. [14] Since strong structural stability results under discretization exist for symplectic methods, [13] it seems desirable to stay within this class, whenever possible.

We point out that for the controlled variable models, there is an elementary separation of variables. Advantages of the resulting mechanical Hamiltonian compared to the original are numerous. It is relatively easy to compare the structure of the original and modified potential energies, thus simplifying the dynamical analysis. Explicit Störmer-Verlet discretization and general higher-order symplectic splitting methods [14] can only be applied to Hamiltonians of this class, holonomic constraints are easily handled (e.g. using SHAKE [15]), and advanced methods such as Rowland’s methods [16] and others based on the potential Hessian [17] are directly applicable. A singularity is eliminated by the separating transformation. The presence of such singularities can affect the growth of terms in the asymptotic Hamiltonian series expansion of a symplectic method, [13] and may therefore contribute to numerical instability.

In Nosé dynamics, the parameterization of time is dependent on the thermostatting variable. In most recent applications of thermostatting, a noncanonical change of variables is employed which eliminates the artificial time scaling (Nosé-Hoover dynamics). The symplectic structure is neglected in numerical simulation but time-reversal symmetry is typically preserved. [4,8] The alternative Poincaré transformation has recently been shown to be as effective, while allowing all calculations to be performed in the canonical Hamiltonian setting. [18] While the use of a separating transformation typically requires a modification of the natural timescale, this does not present much of an obstacle for computation. In Section V, we examine the numerical properties of the mechanical formulation.

II. SEPARATING VARIABLES

Consider an N-degree of freedom system with canonical Hamiltonian of the form

ℋ = g(q)\frac{p^{T} M^{−1}p}{2} + f(q),  (1)

where g is a positive scalar-valued function of position q, f is a potential energy function, and M is a constant mass matrix. A Poincaré time transformation [19] can be used to reduce (1) to mechanical form:

𝒢 = \frac{ℋ − E}{g(q)} = \frac{p^{T} M^{−1}p}{2} + \frac{f(q) − E}{g(q)} ,

where E is a constant representing the value of the energy integral. The problem of finding a separating the variables in the Hamiltonian thus reduces to finding a change variables to the form (1).

Obviously we cannot expect to be able to “flatten out” an arbitrary metric with coordinate and time transformations. However, this is possible for Andersen and Nosé dynamics. The starting point for controlled variable simulations is a mechanical system with N degrees of freedom and the Hamiltonian

H = \frac{p^{T} M^{−1}p}{2} + f(q),  (2)

where M is a positive definite and symmetric mass matrix, and f is a smooth potential energy function. In each of the controlled variable methods, the Hamiltonian is altered by the incorporation of a new control variable, say σ, in the potential, a scaling of the kinetic energy by a σ-dependent function, and the addition of a kinetic term for σ, resulting in

ℋ = g(σ)\frac{p^{T} M^{−1}p}{2} + \frac{π_{σ}^{2}}{2μ} + f̂(q, σ),  (3)

where g(σ) > 0, and μ > 0 is an artificial parameter. Such systems admit separating transformations of the form

σ = ϕ(θ), π_{σ} = ϕ′(θ)^{−1}π_{θ},

hence

ℋ = g(ϕ(θ))\frac{p^{T} M^{−1}p}{2} + \frac{π_{θ}^{2}}{2μ(ϕ′(θ))^{2}} + f̂(q, ϕ(θ)).

To obtain a position-dependent scaling of a constant metric, we must solve the scalar differential equation g(ϕ) = k(ϕ′)^{−2}, for some positive constant k. Since g is positive, the variables can always be separated in this case (at least locally).

III. NOSÉ DYNAMICS

The Nosé extended Hamiltonian based on (2) is

ℋ _{c.t.} = \frac{p^{T} M^{−1}p}{2σ^{2}} + \frac{π_{σ}^{2}}{2μ} + f(q) + gk_{B}T ln σ.  (4)

where k_{B} is the Boltzmann constant, g = N + 1, μ represents a thermal mass assigned to the thermostatting variable σ, and π_{σ} represents its momentum.

We perform a separating change of variables in the Nosé Hamiltonian, in this case replacing σ by its logarithm:

θ = ln σ, π_{σ} = exp(−θ)π_{θ}.

After a rescaling of time by e^{2θ}, the Hamiltonian becomes

𝒢_{c.t.} = \frac{p^{T} M^{−1}p}{2} + \frac{π_{θ}^{2}}{2μ} + e^{2θ}(f(q) + kgT θ − E).

We still must show that the new Hamiltonian recovers the desired statistical-mechanical property; this result is summarized in the following theorem, analogous to previous results. [3,18] In what follows, we take the number of degrees of freedom of the mechanical system to be N_{f} = N, although technically, this number should generally be reduced by three in the statistical mechanics since we typically work in an ensemble constrained by the conservation of linear momentum.

Theorem 1. Assuming ergodicity and with g = N − 1, the mechanical formulation with Hamiltonian 𝒢_{c.t.} samples from the canonical (N,V,T) ensemble.

Proof: We compute the partition function corresponding to 𝒢 along the zero-energy surface:

Z′ = \frac{1}{N!} ∫ dpdqdπ_{θ}dθδ[𝒢_{c.t.}],

for vanishing extended energy. Since the change of variables (q, p, σ, π_{σ}) → (q, p, θ, π_{θ}) is symplectic, hence area conserving, we can directly rewrite the integral in the Nosé extended variables,

Z′ = \frac{1}{N!} ∫ dpdqdπ_{σ}dσδ[σ^{2}(ℋ _{c.t.} − E)],

and the result now follows [18] under the ergodicity assumption, decrementing g by one due to the extra power of σ. □

Theorem 1 indicates that the averages of any particular thermodynamic quantity along the microcanonical dynamics of 𝒢_{c.t.} are equivalent to canonical constant temperature averages associated to (2).

A. Dynamics in the Modified Potential

The modified or ‘isothermal’ potential energy is (with gk_{B}T replaced by an arbitrary positive parameter γ):

f̃(q, θ) = e^{2θ}(γθ + f(q) − E).

The stationary points occur at points (q̄, θ̄), where q̄ is a critical point of f and θ̄ = (E − f(q̄))/γ − 1/2. The stability of these equilibria are determined by the stability of q̄ for (2). A center of the original mechanical system is thus mapped to a corresponding center of the constant temperature dynamics defined by 𝒢_{c.t.}. The modified potential retains the shape of the of the original potential in constant θ slices (since it consists of a shift and a scaling for any θ). On the positive θ side, the motion is uniformly bounded, while potentially lengthy excursions in the negative θ direction are possible.

Looking in the q direction, we have f̃(q, 0) = f(q) − E, and the limits

lim_{θ→+∞} f̃(q, θ) = +∞, lim_{θ→−∞} f̃(q, θ) = 0,

For any value of q, f̃ achieves a unique global minimum in the θ direction, hence the modified potential can be viewed as a “valley” with sides rising in the theta direction. The isothermal (γ = 1, E = 1) potential energy surface corresponding to the harmonic oscillator (f(q) = q^{2}/2) is illustrated in Figure 1. Note that the motion takes place in the regime determined by f̃(q, θ) ≤ 0, since we have shifted the energy to zero as part of the Poincaré transformation.

IV. INTEGRATION IN THE TP ENSEMBLE

Nosé [20] proposed a unified extended Hamiltonian formulation for the constant temperature and pressure (TP) ensemble We were unable to solve the required equations for the separating transformation for this Hamiltonian. Instead, we propose the alternative TP formulation:

ℋ _{t.p.}^{∗} = \frac{p̂^{T} M^{−1}p̂}{2σ^{2}V^{2/3}} + \frac{π_{σ}^{2}}{2μV^{2/3}} + \frac{π_{V}^{2}}{2νσ^{2}}

+f(V^{1/3}q̂) + P V − \frac{k_{B}T ln V}{3} + gk_{B}T ln σ,

where P is the external pressure. A separating transformation for ℋ _{t.p.}^{∗} is constructed as follows. Apply the canonical coordinate transformation of Section III, replacing σ by e^{θ}, and resulting in

ℋ _{t.p.}^{∗} = \frac{p̂^{T} M^{−1}p̂}{2e^{2θ}V^{2/3}} + \frac{π_{θ}^{2}}{2μe^{2θ}V^{2/3}} + \frac{π_{V}^{2}}{2νe^{2θ}}

+f(V^{1/3}q̂) + P V − \frac{k_{B}T ln V}{3} + gk_{B}T θ.

Next, construct a map V = ϕ(ψ), π_{V} = ϕ′(V)^{−1}π_{ψ} by solving

k(ϕ′)^{2} = ϕ^{2/3}.

A particular solution is given by ϕ = ψ^{3/2}, with k = 9/4. Finally make a Poincaré timetransformation, rescaling time by ψ exp(2θ) to obtain the mechanical form:

𝒢_{t.p.}^{∗} = \frac{p̂^{T} M^{−1}p̂}{2} + \frac{π_{θ}^{2}}{2μ} + \frac{2π_{ψ}^{2}}{9ν}

+ψe^{2θ}(f(ψ^{1/2}q̂) + gk_{B}T θ + P ψ^{3/2} − \frac{k_{B}T ln ψ}{2} − E).

We next investigate the relationship between the partition function of extended microcanonical dynamics with the Hamiltonian 𝒢_{t.p.}^{∗} and the partition function for the TP ensemble. Write 𝒢_{t.p.}^{∗} = 𝒢_{t.p.}^{∗}(q̂, p̂, σ, π_{σ}, V, π_{V}); the partition function is

Z = \frac{1}{N!} ∫ δ𝒢_{t.p.}^{∗},

where the integral is taken over the entire admissable extended phase space. Factoring out the time transformation, rewriting in the original control variables, and making the substitutions π̃_{V} = π_{V} /σ, π̃_{σ} = π_{σ}/V^{1/3}, p̃ = p̂/(σV^{1/3}), q̃ = q̂V^{1/3}, we obtain

Z = \frac{1}{N!} ∫ σ^{N+1}V^{1/3}δ [V^{2/3}σ^{2} (H(p̃, q̃) + \frac{π̃_{σ}^{2}}{2μ} + \frac{π̃_{V}^{2}}{2ν} + gk_{B}T ln σ + P V − \frac{k_{B}T ln V}{3} − E)] .

Following Nosé, we now apply the formula, δh(σ) = δ(σ − σ_{∗})/h′(σ_{∗}), which holds provided h has a unique simple zero at σ_{∗}. Thus we find, after a short computation,

Z = \frac{1}{N!} ∫ \frac{σ^{N+1}V^{1/3}}{V^{2/3}σ_{∗}gk_{B}T} δ[σ − σ_{∗}],

which becomes, after integrating over σ,

Z = \frac{1}{N!} ∫ \frac{σ_{∗}^{N}V^{−1/3}}{gk_{B}T}

with

gk_{B}T ln σ_{∗} = E − H(p̃, q̃) − \frac{π̃_{σ}^{2}}{2μ} − \frac{π̃_{V}^{2}}{2ν} − P V + \frac{k_{B}T ln V}{3} .

For g = N, the factor of V^{−1/3} cancels with the last term in the exponent, and equivalence to the TP partition function follows after integrating over the momenta π̃_{V} and π̃_{σ}.

V. SYMPLECTIC INTEGRATION

The implementation of a symplectic method for constant temperature is straightforward once the mechanical formulation is known. From this point on, we direct our comments to Nosé dynamics; numerical considerations for the TP ensemble are expected to be similar.

Denote the time variable in the original mechanical system (2) by t. The flow of Nosé’s extended Hamiltonian then is parameterized by τ with dt/dτ = σ^{−1}. The Nosé-Poincare method [18] and the Nosé-Hoover formulation [4] rescale the time variable according to dτ/dτ̂ = σ, thus recovering the original time variable, since dt/dτ̂ = (dt/dτ)(dτ/dτ̂) = 1. The new formulation (𝒢_{c.t.}) works in yet a different timescale, parameterized by τ̆ with dτ/dτ̆ = σ^{2}, or dt/dτ̆ = σ. This allows recovery of stationary averages 𝒜 (f, t_{∗}) = (∫_{0}^{t_{∗}} f(t)dt)/t_{∗} by making an appropriate change of variables in the integral and adjusting the time interval to the fictive time. For the computation of time-correllated statistics such as diffusion constants or velocity autocorelation functions, it is more convenient to have the solution at fixed real-time steps, requiring the solution of dt/dτ̆ = σ, numerically in tandem with the flow, and interpolating the approximate solution to the fixed timestep mesh.

A. Experiment with a 1D Lennard-Jones Chain

A Störmer-Verlet method was implemented for mechanical Nosé dynamics and applied to a system of N unit-mass particles moving in 1-d with periodic boundary conditions and interacting in a nearest-neighbor Lennard-Jones potential ϕ(r) = r^{−12} − 2r^{−6}. The system was given unit density (length N). An initial condition was chosen at low energy and temperature, the momenta were scaled up in Störmer-Verlet integration until the temperature reached a value near the target, and the system was equilibrated for a few hundred steps. We compared the mechanical formulation with the symplectic Nosé-Poincaré method. [18]

In all runs, the mechanical formulation performed at least comparably to the Nosé-Poincaré method. At lower temperatures T ≤ 1, we observed no substantial qualitative differences between the two approaches. A sample integration, showing the fluctuations in extended energy, mechanical energy and temperature for T = 1 is shown in Fig. 2.

As the target temperature was increased beyond T = 2, the mechanical-Störmer-Verlet scheme showed increasing improvement, particularly in terms of the extended energy error but also, to a somewhat lesser extent, in terms of the convergence of average temperature. Improvements in efficiency (measured in the ratio of work to temperature error) varied from 20 to 30%. On the other hand, at higher temperature, we found that the increased fluctuations of the thermostat made computation of the average temperature more difficult in the mechanical approach, particularly at larger stepsizes, and a more accurate (4th order) interpolation was needed to recover the average.

Acknowledgement.

This work was performed while the author was a guest at the Mathematical Sciences Research Institute (MSRI) in Berkeley, CA. Research at MSRI is supported in part by NSF grant DMS-9701755. The author was also supported during his stay at MSRI by NSF Grant #DMS 9627330). Brian Laird and Robert McLachlan contributed helpful comments.

REFERENCES

[1] Andersen, H.C., Journal of Chemical Physics 72, 2384-2393 (1980)

[2] Nosé, S., Molecular Physics 52, 255-268 (1984)

[3] Nosé, S., Progress in Theoretical Physics Supplement 103, 1-46 (1991)

[4] Hoover, W.G., Physical Review A 31, 1695-1697 (1985)

[5] Toxvaerd, S., and Olsen, O., Physica Scripta T33, 98-101 (1990)

[6] Holian, B., De Groot, A., Hoover, W.G., and Hoover, C.G., Physical Review A 41, 4552-4553 (1990)

[7] Martyna, G., Klein, M., and Tuckerman, M., Journal of Chemical Physics 97, 2635-2643 (1992)

[8] Tuckermann, M.E., Berne, B.J., Martyna, G.J., and Klein, M.L., Journal of Chemical Physics 99, 2796 (1993)

[9] von Boehm, J., and Nieminen, R.M., Physical Review B 50, 6450-6452 (1994)

[10] Cho, K., Joannopoulos, J.D., and Kleinmann, L, Physical Review E 47, 3145 (1993)

[11] Martyna, G., Physical Review E 50, 3234-3236 (1994)

[12] Jang, S., and Voth, G.A., Journal of Chemical Physics 107, 9514-9526 (1997)

[13] G. Benettin and A. Giorgilli, Journal of Statistical Physics 74, 1117 (1994).

[14] Sanz-Serna, J.M., and Calvo, M.P., Numerical Hamiltonian Problems, Chapman and Hall (1994).

[15] Ryckaert, J.P., Ciccotti, G., and Berendsen, H.C., Journal of Computational Physics 23, 327-341, 1977.

[16] Rowlands, G., Journal of Computational Physics 97, 235-239 (1991)

[17] Lopez-Marcos, M.A., Sanz-Serna, J.M., and Skeel, R.D., SIAM Journal on Scientific Computing 18, 223-238 (1997).

[18] Bond, S., Laird, B. and Leimkuhler, B., Journal of Computational Physics, to appear.

[19] Zare, K. and Szebehely, V., Celestial Mechanics 11, 469-482, 1975.

[20] Nosé, S., J. Chem. Phys. 81, 511-515, 1984.

FIG. 1. The isothermal potential f̃ corresponding to the harmonic oscillator f(q) = q^{2}/2 for γ = 1, along the energy slice E = 1.

FIG. 2. Evolution of energies and temperature in the mechanical approach

13