Bulgac-Kusnezov-Nosé-Hoover thermostats
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b October 10, 2018
Bulgac-Kusnezov-Nos´e-Hoover thermostats
Alessandro Sergi ∗ School of Physics, University of KwaZulu-Natal,Pietermaritzburg, Private Bag X01 Scottsville,3209 Pietermaritzburg, South Africa
Gregory S. Ezra † Department of Chemistry and Chemical Biology, Baker Laboratory,Cornell University, Ithaca, New York 14853, USA
Abstract
In this paper we formulate Bulgac-Kusnezov constant temperature dynamics in phase spaceby means of non-Hamiltonian brackets. Two generalized versions of the dynamics are similarlydefined: one where the Bulgac-Kusnezov demons are globally controlled by means of a singleadditional Nos´e variable, and another where each demon is coupled to an independent Nos´e-Hooverthermostat. Numerically stable and efficient measure-preserving time-reversible algorithms arederived in a systematic way for each case. The chaotic properties of the different phase space flowsare numerically illustrated through the paradigmatic example of the one-dimensional harmonicoscillator. It is found that, while the simple Bulgac-Kusnezov thermostat is apparently not ergodic,both of the Nos´e-Hoover controlled dynamics sample the canonical distribution correctly. ∗ Electronic address: [email protected] † Electronic address: [email protected] . INTRODUCTION In condensed matter studies, there are many situations in which molecular dynamics sim-ulation at constant-temperature [1–3] is needed. For example, this occurs when magneticsystems are modelled in terms of classical spins [4–7]. Deterministic methods [8–10], basedon non-Hamiltonian dynamics [11–19], can sample the canonical distribution provided thatthe motion in the phase space of the relevant degrees of freedom is ergodic [1, 3]. However,classical spin systems are usually formulated in terms of non-canonical variables [20, 21],without a kinetic energy expressed through momenta in phase space, so that Nos´e dynamicscannot be applied directly. To tackle this problem, Bulgac and Kusnezov (BK) introduceda deterministic constant-temperature dynamics [22–24] which can be applied to spins. Anumber of numerical approaches to integration of spin dynamics can be found in the lit-erature [25–28]. However, BK dynamics, as any other deterministic canonical phase spaceflow, is able to correctly sample the canonical distribution only if the motion in phase spaceis ergodic on the timescale of the simulation. In general, this condition is very difficult tocheck for statistical systems with many degrees of freedom, while it is known that, despiteits simplicity, the one-dimensional harmonic oscillator provides a difficult and importantchallenge for deterministic thermostatting methods [9, 29–31].In this paper, we accomplish two goals. First, by reformulating BK dynamics throughnon-Hamiltonian brackets [14, 15] in phase space, we introduce two generalized versionsof the BK time evolution which are able to sample the canonical distribution for a stiffharmonic system. Second, using a recently introduced approach based on the geometryof non-Hamiltonian phase space [19], we are able to derive stable and efficient measure-preserving and time-reversible algorithms in a systematic way for all the phase space flowstreated here.The BK phase space flow introduces temperature control by means of fictitious coor-dinates (and their associated momenta in an extended phase space) traditionally called‘demons’. Our generalizations of the BK dynamics are obtained by controlling the BKdemons themselves by means of additional Nos´e-type variables [8]. In one case, the BKdemons are controlled globally by means of a single additional Nos´e-Hoover thermostat [8, 9].In the following this will be referred to as (BKNH) Bulgac-Kusnezov-Nos´e-Hoover dynam-ics. In the second case, each demon is coupled to an independent Nos´e-Hoover thermostat.2his will be called the Bulgac-Kusnezov-Nos´e-Hoover chain (BKNHC), and corresponds to‘massive’ NH thermostatting of the demon variables [32].The ability to derive numerically stable measure-preserving time-reversible algorithms[19] for Nos´e controlled BK dynamics is very encouranging for future applications to ther-mostatted spin systems.This paper is organized as follows. In Sec. II we briefly sketch the unified formalism fornon-Hamiltonian phase space flows and measure-preserving integration. The BK dynamicsis formulated in phase space and a measure-preserving integration algorithm is derived inSec. III. The BKNH and BKNH-chain thermostats are treated in Secs. IV and V respectively.Numerical results for the one-dimensional harmonic oscillator using these thermostats arepresented and discussed in Sec. VI. Section VII reports our conclusions.In addition we include several Appendices. A useful operator formula is derived in Ap-pendix A, while invariant measures for the BK, BKNH, and BKNHC phase space flows arederived in Appendices B, C, and D, respectively.3
I. NON-HAMILTONIAN BRACKETS AND MEASURE-PRESERVING ALGO-RITHMS
Consider an arbitrary system admitting a time-independent (extended) Hamiltonian ex-pressed in terms of the phase space coordinates x i , i = 1 , . . . , N . In this case, the Hamil-tonian can be interpreted as the conserved energy of the system.Upon introducing an antisymmetric tensor field (generalized Poisson tensor [21, 33]) inphase space, B ( x ) = − B T ( x ), one can define non-Hamiltonian brackets [14–16] as { a, b } = n X i,j =1 ∂a∂x i B ij ∂b∂x j , (1)where a = a ( x ) and b = b ( x ) are two arbitrary phase space functions. The bracket definedin Eq. (1) is classified as non-Hamiltonian [14–16] since, in general, it does not obey theJacobi relation, i.e., in general the Jacobiator J 6 = 0, where [21] J = { a, { b, c }} + { b, { c, a }} + { c, { a, b }} , (2)with c = c ( x ) arbitrary phase space function (in addition to the functions a and b , previouslyintroduced). If J 6 = 0, the tensor B ij is said to define an ‘almost-Poisson’ structure [34].(Such systems have also been called ‘pseudo-Hamiltonian’ [33].)An energy-conserving and in general non-Hamiltonian phase space flow is then definedby the vector field ˙ x i = { x i , H } = N X j =1 B ij ∂H∂x j , (3)where conservation of H ( x ) follows directly from the antisymmetry of B ij .It has previously been shown how equilibrium statistical mechanics can be comprehen-sively formulated within this framework [16]. It is also possible to recast the above formalismand the corresponding statistical mechanics in the language of differential forms [17, 18]. Ifthe matrix B is invertible (this is true for all the cases considered here), with inverse Ω ij ,we can define the 2-form [35] Ω = Ω ij dx i ∧ dx j . (4)The dynamics of Eq. (3) is then Hamiltonian if and only if the form (4) is closed , i.e., haszero exterior derivative, d Ω = 0 [35]. This condition is independent of the particular systemof coordinates used to describe the dynamics.4he structure of Eq. (3) can be taken as the starting point for derivation of efficienttime-reversible integration algorithms that also preserve the appropriate measure in phasespace [19]. Measure-preserving algorithms can be derived upon introducing a splitting ofthe Hamiltonian H = n s X α =1 H α (5)which in turn induces a splitting of the Liouville operator associated with the non-Hamiltonian bracket in Eq. (1), L α x i = { x i , H α } = N X j =1 B ij ∂H α ∂x j . (6)When the phase space flow has a non-zero compressibility κ = N X i,j =1 ∂ B ij ∂x i ∂H∂x j (7)the statistical mechanics must be formulated in terms of a modified phase space measure[12–18] ω = e − w ( x ) ω (8)where ω = dx ∧ dx ∧ . . . ∧ dx N (9)is the standard phase space volume element (volume form [35]) and the statistical weight w ( x ) is defined by dwdt = κ ( x ) . (10)It has been shown that, provided the condition ∂∂x j (cid:2) e − w ( x ) B ij (cid:3) = 0 , i = 1 , . . . N (11)is satisfied, then L α ω = 0 for every α , (12)so that the volume element ω is invariant under each of the L α [19]. The condition (11) issatisfied for all the cases considered below, so that, exploiting the decomposition in Eq. (6),5lgorithms derived by means of a symmetric Trotter factorization of the Liouville propagator:exp[ τ L ] = n s − Y α =1 exp h τ L α i exp exp [ τ L n x ] × n s − Y β =1 exp h τ L n s − β i (13)are not only time-reversible but also measure-preserving.6 II. PHASE SPACE FORMULATION OF THE BK THERMOSTAT
A phase space formulation of the BK thermostat can be achieved upon introducing theHamiltonian H BK = p m + V ( q ) + K ( p ζ ) m ζ + K ( p ξ ) m ξ + k B T ( ζ + ξ ) (14a)= H ( q, p ) + K ( p ζ ) m ζ + K ( p ξ ) m ξ + k B T ( ζ + ξ ) , (14b)where ( q, p ) are the physical degrees of freedom (coordinates and momenta), with mass m , to be simulated at constant temperature T , while ζ and ξ are the BK ‘demons’, withcorresponding inertial parameters m ζ and m ξ , and associated momenta ( p ζ , p ξ ) [22–24]. K and K provide the kinetic energy of demon variables, and for the moment are left arbitrary.Upon defining the phase space point as x = ( q, ζ , ξ, p, p ζ , p ξ ) = ( x , x , x , x , x , x ), onecan introduce an antisymmetric BK tensor field as B BK = − G ∂G ∂p
00 0 0 0 0 ∂G ∂q − − G − ∂G ∂p G G − ∂G ∂q (15)where G and G are functions of system variables ( p, q ) only.Substituting B BK and H BK into Eq. (3), we obtain the energy-conserving equations˙ q = ∂H∂p − G ( q, p ) m ξ ∂K ∂p ξ (16a)˙ ζ = 1 m ζ ∂G ∂p ∂K ∂p ζ (16b)˙ ξ = 1 m ξ ∂G ∂q ∂K ∂p ξ (16c)˙ p = − ∂H∂q − G ( q, p ) m ζ ∂K ∂p ζ (16d)˙ p ζ = G ∂H∂p − k B T ∂G ∂p (16e)˙ p ξ = G ∂H∂q − k B T ∂G ∂q . (16f)The associated invariant measure for the BK flow is discussed in Appendix B.7 . Algorithm for BK Dynamics In order to derive a measure preserving algorithms, the first step, following Eq. (5), is tointroduce a splitting of H BK : H BK1 = V ( q ) (17a) H BK2 = p m (17b) H BK3 = k B T ζ (17c) H BK4 = k B T ξ (17d) H BK5 = K ( p ζ ) m ζ (17e) H BK6 = K ( p ξ ) m ξ . (17f)A measure-preserving splitting of the Liouville operator then follows from Eq. (6): L BK1 = − ∂V∂q ∂∂p + G ∂V∂q ∂∂p ξ (18a) L BK2 = pm ∂∂q + G pm ∂∂p ζ (18b) L BK3 = − k B T ∂G ∂p ∂∂p ζ (18c) L BK4 = − k B T ∂G ∂q ∂∂p ξ (18d) L BK5 = 1 m ζ ∂G ∂p ∂K ∂p ζ ∂∂ζ − G m ζ ∂K ∂p ζ ∂∂p (18e) L BK6 = − G m ξ ∂K ∂p ξ ∂∂q + 1 m ξ ∂G ∂q ∂K ∂p ξ ∂∂ξ . (18f)Upon choosing a symmetric Trotter factorization of the BL Liouville operator based on thedecomposition L BK = X α =1 L BK α (19)a measure-preserving algorithm can be produced in full generality.In practice, a choice of K , K , G , G must be made in order obtain explicit formulas.8n this paper, we make the following simple choices: G = p (20a) G = q (20b) K = p ζ K = p ξ . (20d)In terms of Eqs (20a–20d), the antisymmetric BK tensor becomes˜ B BK = − q − − p − p q − , (21)and the Hamiltonian reads˜ H BK = H ( q, p ) + p ζ m ζ + p ξ m ξ + k B T ( ζ + ξ ) . (22)The split Liouville operators now simplify as follows:˜ L BK1 = − ∂V∂q ∂∂p + q ∂V∂q ∂∂p ξ (23a)˜ L BK2 = pm ∂∂q + p m ∂∂p ζ (23b)˜ L BK3 = − k B T ∂∂p ζ (23c)˜ L BK4 = − k B T ∂∂p ξ (23d)˜ L BK5 = p ζ m ζ ∂∂ζ − p ζ m ζ p ∂∂p + p ζ m ζ ∂∂p η (23e)˜ L BK6 = − p ξ m ξ q ∂∂q + p ξ m ξ ∂∂ξ + p ξ m ξ ∂∂p χ (23f)For the purposes of defining an efficient integration algorithm, we combine commuting9iouville operators as follows: L BK A ≡ ˜ L BK1 + ˜ L BK4 = F ( q ) ∂∂p + F p ξ ∂∂p ξ (24a) L BK B ≡ ˜ L BK2 + ˜ L BK3 = pm ∂∂q + F p ζ ∂∂p ζ (24b) L BK C ≡ ˜ L BK5 + ˜ L BK6 = − p ζ m ζ p ∂∂p − p ξ m ξ q ∂∂q + p ζ m ζ ∂∂ζ + p ξ m ξ ∂∂ξ (24c)where F ( q ) = − ∂V /∂q (25a) F p ξ = q ∂V∂q − k B T (25b) F p ζ = p m − k B T . (25c)Defining U BK α ( τ ) = exp h τ ˜ L BK α i , (26)where α = A, B, C , one possible reversible measure-preserving integration algorithm for theBK thermostat is then U ( τ ) BK = U BK B (cid:16) τ (cid:17) U BK C (cid:16) τ (cid:17) U BK B (cid:16) τ (cid:17) × U BK A ( τ ) × U BK B (cid:16) τ (cid:17) U BK C (cid:16) τ (cid:17) U BK B (cid:16) τ (cid:17) . (27)Using the so-called direct translation technique [36] we can expand the above symmetricbreak-up of the Liouville operator into a pseudo-code form, ready to be implemented on thecomputer: • q → q + τ pm p ζ → p ζ + τ F p ζ : U BK B (cid:0) τ (cid:1) p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ : U BK C (cid:0) τ (cid:1) • q → q + τ pm p ζ → p ζ + τ F p ζ : U BK B (cid:0) τ (cid:1) • p → p + τ Fp ξ → p ξ + τ F p ξ : U BK A ( τ ) • q → q + τ pm η → η + τ p η m η p ζ → p ζ + τ F p ζ : U BK B (cid:0) τ (cid:1) • p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ : U BK C (cid:0) τ (cid:1) • q → q + τ pm p ζ → p ζ + τ F p ζ : U BK B (cid:0) τ (cid:1) V. BULGAC-KUSNEZOV-NOS´E-HOOVER DYNAMICS
The BKNH Hamiltonian H BKNH = H ( q, p ) + K ( p ζ ) m ζ + K ( p ξ ) m ξ + p η m η + k B T ( ζ + ξ ) + 2 k B T η (28)is simply the BK Hamiltonian augmented by the Nos´e variables ( η, p η ) with mass m η . Withthe antisymmetric BKNH tensor B BKNH = − G
00 0 0 0 0 ∂G ∂p ∂G ∂q
00 0 0 0 0 0 0 1 − − G − ∂G ∂p G − p ζ G − ∂G ∂q − p ξ − p ζ p ξ . (29)we obtain from Eq. (3) equations of motion for the phase space variables x =( q, ζ , ξ, η, p, p ζ , p ξ , p η ) = ( x , x , x , x , x , x , x , x ):˙ q = ∂H∂p − G ( q, p ) m ξ ∂K ∂p ξ (30a)˙ ζ = 1 m ζ ∂G ∂p ∂K ∂p ζ (30b)˙ ξ = 1 m ξ ∂G ∂q ∂K ∂p ξ (30c)˙ η = p η m η (30d)˙ p = − ∂H∂q − G ( q, p ) m ζ ∂K ∂p ζ (30e)˙ p ζ = G ∂H∂p − k B T ∂G ∂p − p ζ p η m η (30f)˙ p ξ = G ∂H∂q − k B T ∂G ∂q − p ξ p η m η (30g)˙ p η = p ζ m ζ ∂K ∂p ζ + p ξ m ξ ∂K ∂p ξ − k B T . (30h)Here, a single Nos´e variable is coupled to both of the BK demons ζ and ξ . The associatedinvariant measure is discussed in Appendix C.12 . Algorithm for BKNH dynamics The Hamiltonian can be split as H BKNH1 = V ( q ) (31a) H BKNH2 = p m (31b) H BKNH3 = k B T ζ (31c) H BKNH4 = k B T ξ (31d) H BKNH5 = K ( p ζ ) m ζ (31e) H BKNH6 = K ( p ξ ) m ξ (31f) H BKNH7 = p η m η (31g) H BKNH8 = 2 k B T η (31h)The measure-preserving splitting [19] of the Liouville operator L α = B BKNH ij ∂H BKNH α ∂x j ∂∂x i (32)yields L BKNH1 = − ∂V∂q ∂∂p + G ∂V∂q ∂∂p ξ (33a) L BKNH2 = pm ∂∂q + G pm ∂∂p ζ (33b) L BKNH3 = − k B T ∂G ∂p ∂∂p ζ (33c) L BKNH4 = − k B T ∂G ∂q ∂∂p ξ (33d) L BKNH5 = 1 m ζ ∂G ∂p ∂K ∂p ζ ∂∂ζ − G m ζ ∂K ∂p ζ ∂∂p + p ζ m ζ ∂K ∂p ζ ∂∂p η (33e) L BKNH6 = − G m ξ ∂K ∂p ξ ∂∂q + 1 m ξ ∂G ∂q ∂K ∂p ξ ∂∂ξ + p ξ m ξ ∂K ∂p ξ ∂∂p η (33f) L BKNH7 = p η m η ∂∂η − p η m η p ζ ∂∂p ζ − p η m η p ξ ∂∂p ξ (33g) L BKNH8 = − k B T ∂∂p η . (33h)13t this stage, we leave the general formulation and adopt the particular choice of K , K , G , and G given in Eq. (20). The antisymmetric BKNH tensor becomes˜ B BKNH = − q
00 0 0 0 0 1 0 00 0 0 0 0 0 1 00 0 0 0 0 0 0 1 − − p − p − p ζ q − − p ξ − p ζ p ξ , (34)and the Hamiltonian simplifies to˜ H BKNH = H ( q, p ) + p ζ m ζ + p ξ m ξ + p η m η + k B T ( ζ + ξ ) + 2 B T η . (35)The split Liouville operators are now˜ L BKNH1 = − ∂V∂q ∂∂p + q ∂V∂q ∂∂p ξ (36a)˜ L BKNH2 = pm ∂∂q + p m ∂∂p ζ (36b)˜ L BKNH3 = − k B T ∂∂p ζ (36c)˜ L BKNH4 = − k B T ∂∂p ξ (36d)˜ L BKNH5 = p ζ m ζ ∂∂ζ − p ζ m ζ p ∂∂p + p ζ m ζ ∂∂p η (36e)˜ L BKNH6 = − p ξ m ξ q ∂∂q + p ξ m ξ ∂∂ξ + p ξ m ξ ∂∂p η (36f)˜ L BKNH7 = p η m η ∂∂η − p η m η p ζ ∂∂p ζ − p η m η p ξ ∂∂p ξ (36g)˜ L BKNH8 = − k B T ∂∂p η . (36h)For the purposes of defining an efficient integration algorithm, we combine commuting Li-14uville operators as follows: L BKNH A ≡ ˜ L BKNH1 + ˜ L BKNH4 + ˜ L BKNH7 = F ( q ) ∂∂p + p η m η ∂∂η − p η m η p ζ ∂∂p ζ + (cid:18) − p χ m χ p ξ + F p ξ (cid:19) ∂∂p ξ (37a)˜ L BKNH B ≡ ˜ L BKNH2 + ˜ L BKNH3 = pm ∂∂q + F p ζ ∂∂p ζ (37b) L BKNH C ≡ ˜ L BKNH5 + ˜ L BKNH6 + ˜ L BKNH8 = − p ζ m ζ p ∂∂p − p ξ m ξ q ∂∂q + p ζ m ζ ∂∂ζ + p ξ m ξ ∂∂ξ + F p η ∂p η , (37c)where F ( q ) = − ∂V∂q (38a) F p ξ = q ∂V∂q − k B T (38b) F p ζ = p m − k B T (38c) F p η = p ζ m ζ + p ξ m ξ − k B T . (38d)In L A there appears an operator with the form L i = (cid:18) − p k m k p i + F p i (cid:19) ∂∂p i , (39)where ( k, i ) = ( χ, ξ ) for L A . The action of the propagator associated with this operator on p i is derived in Appendix A, and is given by e τL i p i = p i e − τ pkmk + τ F p i e − τ pk mk (cid:18) τ p k m k (cid:19) − sinh (cid:20) τ p k m k (cid:21) . (40)The apparently singular function (cid:18) τ p k m k (cid:19) − sinh (cid:20) τ p k m k (cid:21) (41)is in fact well behaved as p k →
0, and can be expanded in a Maclaurin series to suitablyhigh order [37]. In our implementation we used an eighth order expansion.The propagators for the BKNH dynamics can now be defined as U BKNH α ( τ ) = exp h τ ˜ L BKNH α i , (42)15here α = A, B, C . One possible reversible measure-preserving integration algorithm forthe BKNH thermostat can then be derived from the following Trotter factorization: U ( τ ) BKNH = U BKNH B (cid:16) τ (cid:17) U BKNH C (cid:16) τ (cid:17) U BKNH B (cid:16) τ (cid:17) × U BKNH A ( τ ) × U BKNH B (cid:16) τ (cid:17) U BKNH C (cid:16) τ (cid:17) U BKNH B (cid:16) τ (cid:17) . (43)The direct translation technique gives the following pseudo-code: • q → q + τ pm p ζ → p ζ + τ F p ζ : U BKNH B (cid:0) τ (cid:1) • p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ p η → p η + τ F p ζ : U BKNH C (cid:0) τ (cid:1) • q → q + τ pm p ζ → p ζ + τ F p ζ : U BKNH B (cid:0) τ (cid:1) • p → p + τ F ( q ) p ξ → p ξ + τ F p ξ η → η + τ p η m η p ζ → p ζ exp h − τ p η m η i : U BKNH A ( τ ) • q → q + τ pm p ζ → p ζ + τ F p ζ : U BKNH B (cid:0) τ (cid:1) • p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ p η → p η + τ F p η : U BKNH C (cid:0) τ (cid:1) q → q + τ pm p ζ → p ζ + τ F p ζ : U BKNH B (cid:0) τ (cid:1) . BULGAC-KUSNEZOV-NOS´E-HOOVER CHAIN For simplicity, we explicitly treat only the case in which the p ζ and p ξ demons are eachcoupled to a standard NH thermostat (length one). It would be straightforward to coupleeach of the demons to NH chains [32], and the general case can be easily inferred from whatfollows. Define the Hamiltonian H BKNHC = H ( q, p ) + K ( p ζ ) m ζ + K ( p ξ ) m ξ + p η m η + p χ m χ + k B T ( ζ + ξ + η + χ ) . (44)Upon defining the phase space point x = ( q, ζ , ξ, η, χ, p, p ζ , p ξ , p η , p χ ) =( x , x , x , x , x , x , x , x , x , x ) and the antisymmetric BKNHC tensor B BKNHC = − G ∂G ∂p ∂G ∂q − − G − ∂G ∂p G − p ζ G − ∂G ∂q − p ξ − p ζ − p ξ , (45)associated non-Hamiltonian equations of motion are˙ x i = B BKNHC ij ∂H BKNHC ∂x j (46)with i = 1 , . . . ,
10. 18 . Algorithm for BKNHC chain dynamics
Splitting the BKNHC chain Hamiltonian as H BKNHC1 = V ( q ) (47a) H BKNHC2 = p m (47b) H BKNHC3 = k B T ζ (47c) H BKNHC4 = k B T ξ (47d) H BKNHC5 = K ( p ζ ) m ζ (47e) H BKNHC6 = K ( p ξ ) m ξ (47f) H BKNHC7 = p η m η (47g) H BKNHC8 = k B T η (47h) H BKNHC9 = p χ m χ (47i) H BKNHC10 = k B T χ , (47j)we obtain the corresponding measure-preserving splitting of the Liouville operator L α = B BKNHC ij ∂H BKNHC α ∂x j ∂∂x i . (48)At this stage we go directly to Eqs (20). The antisymmetric Nos´e-Hoover-Bulgac-Kusnezov tensor becomes˜ B BKNHC = − q − − p − p − p ζ q − − p ξ − p ζ − p ξ , (49)19he Hamiltonian˜ H BKNHC = H ( q, p ) + p ζ m ζ + p ξ m ξ + p η m η + p χ m χ + k B T ( ζ + ξ + η + χ ) (50)and associated Liouville operators˜ L BKNHC1 = − ∂V∂q ∂∂p + q ∂V∂q ∂∂p ξ (51a)˜ L BKNHC2 = pm ∂∂q + p m ∂∂p ζ (51b)˜ L BKNHC3 = − k B T ∂∂p ζ (51c)˜ L BKNHC4 = − k B T ∂∂p ξ (51d)˜ L BKNHC5 = p ζ m ζ ∂∂ζ − p ζ m ζ p ∂∂p + p ζ m ζ ∂∂p η (51e)˜ L BKNHC6 = − p ξ m ξ q ∂∂q + p ξ m ξ ∂∂ξ + p ξ m ξ ∂∂p χ (51f)˜ L BKNHC7 = p η m η ∂∂η − p η m η p ζ ∂∂p ζ (51g)˜ L BKNHC8 = − k B T ∂∂p η (51h)˜ L BKNHC9 = p χ m χ ∂∂χ − p χ m χ p ξ ∂∂p ξ (51i)˜ L BKNHC10 = − k B T ∂∂p χ . (51j)We combine commuting Liouville operators as follows: L BKNHC A ≡ ˜ L BKNHC1 + ˜ L BKNHC4 + ˜ L BKNHC9 = F ( q ) ∂∂p + p χ m χ ∂∂χ + (cid:18) − p χ m χ p ξ + F p ξ (cid:19) ∂∂p ξ (52a) L BKNHC B ≡ ˜ L BKNHC2 + ˜ L BKNHC3 + ˜ L BKNHC7 = pm ∂∂q + p η m η ∂∂η + (cid:18) − p η m η p ζ + F p ζ (cid:19) ∂∂p ζ (52b) L BKNHC C ≡ ˜ L BKNHC5 + ˜ L BKNHC6 + ˜ L BKNHC8 + ˜ L BKNHC10 = − p ζ m ζ p ∂∂p − p ξ m ξ q ∂∂q + p ζ m ζ ∂∂ζ + p ξ m ξ ∂∂ξ + F p η ∂p η + F p χ ∂p χ , (52c)20here F ( q ) = − ∂V∂q (53a) F p ξ = q ∂V∂q − k B T (53b) F p ζ = p m − k B T (53c) F p η = p ζ m ζ − k B T (53d) F p χ = p ξ m ξ − k B T . (53e)Both in L BKNHC A and L BKNHC B there appears an operator with the form L i = (cid:18) − p k m k p i + F p i (cid:19) ∂∂p i , (54)where ( k, i ) = ( χ, ξ ) for L A and ( k, i ) = ( η, ζ ) for L B . Again following the derivation inAppendix A, we find e τL i p i = p i e − τ pkmk + τ F p i e − τ pk mk (cid:18) τ p k m k (cid:19) − sinh (cid:20) τ p k m k (cid:21) . (55)The function (cid:16) τ p k m k (cid:17) − sinh h τ p k m k i is treated through an eighth order expansion [37].The propagators U BKNHC α ( τ ) = exp h τ ˜ L BKNHC α i (56)with α = A, B, C can now be introduced. One possible reversible measure-preserving inte-gration algorithm for the BKNHC chain thermostat is then U ( τ ) BKNHC = U BKNHC B (cid:16) τ (cid:17) U BKNHC C (cid:16) τ (cid:17) U BKNHC B (cid:16) τ (cid:17) × U BKNHC A ( τ ) × U BKNHC B (cid:16) τ (cid:17) U BKNHC C (cid:16) τ (cid:17) U BKNHC B (cid:16) τ (cid:17) . (57)In pseudo-code form, we have the resulting integration algorithm: • q → q + τ pm η → η + τ p η m η p ζ → p ζ e − τ pηmη + τ F p ζ e − τ pη mη (cid:16) τ p η m η (cid:17) − sinh h τ p η m η i : U BKNHC B (cid:0) τ (cid:1) p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ p η → p η + τ F p ζ p χ → p χ + τ F p χ : U BKNHC C (cid:0) τ (cid:1) • q → q + τ pm η → η + τ p η m η p ζ → p ζ e − τ pηmη + τ F p ζ e − τ pη mη (cid:16) τ p η m η (cid:17) − sinh h τ p η m η i : U BKNHC B (cid:0) τ (cid:1) • p → p + τ Fχ → χ + τ p χ m χ p ξ → p ξ e − τ pχmχ + τ F p ξ e − τ pχ mχ (cid:16) τ p χ m χ (cid:17) − sinh h τ p χ m χ i : U BKNHC A ( τ ) • q → q + τ pm η → η + τ p η m η p ζ → p ζ e − τ pηmη + τ F p ζ e − τ pη mη (cid:16) τ p η m η (cid:17) − sinh h τ p η m η i : U BKNHC B (cid:0) τ (cid:1) • p → p exp h − τ p ζ m ζ i q → q exp h − τ p ξ m ξ i ζ → ζ + τ p ζ m ζ ξ → ξ + τ p ξ m ξ p η → p η + τ F p η p χ → p χ + τ F p χ : U BKNHC C (cid:0) τ (cid:1) • q → q + τ pm η → η + τ p η m η p ζ → p ζ e − τ pηmη + τ F p ζ e − τ pη mη (cid:16) τ p η m η (cid:17) − sinh h τ p η m η i : U BKNHC B (cid:0) τ (cid:1) I. NUMERICAL RESULTS
In its simplicity, the dynamics of a harmonic mode in one dimension is a paradigmaticexample for checking the chaotic (ergodic) properties of constant-temperature phase spaceflows and the correct sampling of the canonical distribution. It is well known that it isnecessary to generalize basic Nos´e-Hoover dynamics [1, 8, 9] to thermostats such as theNos´e-Hoover chain [32, 37] in order to produce correct constant-temperature averages forsystems such as the harmonic oscillator.Some time ago, BK dynamics was devised to provide a deterministic thermostat forsystems such as classical spins [23, 24]. To ensure efficient thermostatting, BK found itnecessary to introduce several ‘demons’ per thermostatted degree of freedom, where eachdemon was taken to have a different and in principle complicated coupling to the systemdegree of freedom [23, 24]. In the present work, we keep the form of the system-thermostatcoupling as simple as possible, in order to facilitate the formulation of explicit, reversibleand measure-preserving integrators [19]. It is then of interest to investigate the ability of ourBK-type thermostats to produce the correct canonical sampling in the case of the harmonicoscillator. Interest in harmonic modes is also justified by the possibility of devising modelsof condensed matter systems in terms of coupled spins and harmonic modes, as already donein quantum dynamics with so-called spin-boson models [38]. We therefore investigate theperformance of our integration schemes on the simple one-dimensional harmonic oscillator.For the particular calculations reported here, the oscillator angular frequency, all massesand k B T were taken to be unity. The time step in all cases was τ = 0 . time steps, starting from the same initial conditions: harmonic oscillatorcoordinate q = 0 .
3, all other phase space variables zero at t = 0.The measure-preserving algorithms derived here result in stable numerical integration forall the three cases treated: BK, BKNH, and BKNHC chain dynamics. Figure 1 shows thethree extended Hamiltonians (normalized by their respective initial time value) versus time.All three Hamiltonians are numerically conserved by the corresponding measure-preservingalgorithm to very high accuracy (which is maintained in all the three cases).However, the basic BK phase space flow is not capable of producing the correct canonicalsampling for a harmonic mode. This can be easily checked since the canonical distributionfunction of the harmonic oscillator is isotropic in phase space and its radial dependence23an be calculated exactly. Details of this way of visualizing the phase space sampling havealready been given in [14, 15]. Figure 2, displaying the comparison between the theoreticaland the calculated value of the radial probability in phase space, clearly shows that theBK dynamics is not able to produce canonical sampling. A look at the inset of Fig. 2,showing the phase space distribution for the harmonic mode, also immediately shows thatthe dynamics is not ergodic.The same analysis has been carried out for BKNH and BKNHC phase space flows, andthese are displayed in Fig 3 and Fig 4, respectively. Within numerical errors, both BKNHand BKNHC thermostats are able to produce the correct canonical distribution for the stiffharmonic modes.Introduction of a single, global Nos´e-type variable in the BKNH thermostat effectivelyintroduces additional coupling between the two demon variables. The effectiveness of theBKNH thermostat is consistent with our findings (results not discussed here) that introduc-tion of explicit coupling between demons in BK thermostat dynamics also leads to efficientthermostatting of the harmonic oscillator. 24 II. CONCLUSIONS
We have formulated Bulgac-Kusnezov [23, 24], Nos´e-Hoover controlled Bulgac-Kusnezov,and Bulgac-Kusnezov-Nos´e-Hoover chain thermostats in phase space by means of non-Hamiltonian brackets [14, 15]. We have derived time-reversible measure-preserving algo-rithms [19] for these three cases and showed that additional control by a single Nos´e-Hooverthermostat or independent Nos´e-Hoover thermostats is necessary to produce the correctcanonical distribution for a stiff harmonic mode.Measure-preserving dynamics of the kind discussed here is associated with equilibriumsimulations (where, for example, there is a single temperature parameter T ). Stationaryphase space distributions associated with non-equilibrium situations are much more com-plicated than the smooth equilbrium densities analyzed in the present paper [11, 39, 40].Nonequilibrium simulations of heat flow could be carried out by extending the present ap-proach to multimode systems (e.g., a chain of oscillators) coupled to BK-type demons withassociated NH thermostats corresponding to two different temperatures [41–43].The techniques presented here for derivation and implementation of thermostats havebeen shown to be efficient and versatile. We anticipate that analogous approaches can beusefully applied to systems of classical spins coupled to both harmonic and anharmonicmodes. 25 ppendix A: Operator formula We wish to determine the action of the propagator associated with the Liouville operatorEq. (39). This is equivalent to solving the evolution equation (recall i = k ) dp i dt = (cid:18) − p k m k p i + F p i (cid:19) (A1)from t = 0 to t= τ . Integrating, we have − m k p k ln (cid:18) − p k m k p i + F p i (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) τ = τ (A2)giving p i ( τ ) ≡ exp (cid:20) τ (cid:18) − p k m k p i + F p i (cid:19) ∂∂p i (cid:21) p i (A3a)= p i e − τp k /m k + m k p k F p i (cid:0) − e − τp k /m k (cid:1) (A3b)= p i e − τp k /m k + τ F p i e − τ pk mk sinh h τ p k m k i τ p k m k . (A3c) Appendix B: Invariant Measure of the BK phase space flows
The phase space compressibility of the phase space BK thermostat is κ BK = ∂ B BK ij ∂x i ∂H BK ∂x i = − m ζ ∂G ∂p ∂K ∂p ζ − m ξ ∂G ∂q ∂K ∂p ξ (B1)Upon introducing the function H BKT = H + K m ζ + K m ξ (B2)one can easily find that κ BK = 1 k B T dH
BKT dt (B3)so that the invariant measure in phase space reads dµ = dx exp (cid:20) − Z t dtκ BK (cid:21) (B4a)= dx exp (cid:2) − βH BKT (cid:3) (B4b)= dx exp[ − βH BK ] exp[ ζ + ξ ] (B4c)as desired. 26 ppendix C: Invariant Measure of the BKNH phase space flows The phase space compressibility of the NH controlled Bulgac-Kusnezov thermostat is κ BKNH = ∂ B BKNH ij ∂x i ∂H BKNH ∂x i = − m ζ ∂G ∂p ∂K ∂p ζ − m ξ ∂G ∂q ∂K ∂p ξ − p η m η (C1)Upon introducing the function H BKNHT = H + K m ζ + K m ξ + p η m η (C2)we have κ BKNH = 1 k B T dH
BKT dt (C3)so that the invariant measure in phase space is dµ = dx exp (cid:20) − Z t dtκ BKNH (cid:21) (C4a)= dx exp (cid:2) − βH BKNHT (cid:3) (C4b)= dx exp[ − βH BKNH ] exp[ ζ + ξ + 2 η ] . (C4c) Appendix D: Invariant Measure of the BKNHC chain phase space flows
The phase space compressibility of the Nos´e-Hoover-Bulgac-Kusnezov chain is κ BKNHC = ∂ B BKNHC ij ∂x i ∂H BKNHC ∂x i = − m ζ ∂G ∂p ∂K ∂p ζ − m ξ ∂G ∂q ∂K ∂p ξ − p η m η − p χ m χ (D1)Upon introducing the function H BKNHCT = H + K m ζ + K m ξ + p η m η + p χ m χ (D2)we have κ BKNHC = 1 k B T dH
BKT dt (D3)so that the invariant measure in phase space reads dµ = dx exp (cid:20) − Z t dtκ BKNHC (cid:21) (D4a)= dx exp (cid:2) − βH BKNHCT (cid:3) (D4b)= dx exp[ − βH BKNHC ] exp[ ζ + ξ + η + χ ] . (D4c)27
1] S. Nos´e, Prog. Theo. Phys. Suppl. , 1 (1991).[2] D. Frenkel and B. Smit,
Understanding Molecular Simulation: From Algorithms to Applica-tions (Academic, New York, 2001), 2nd ed.[3] B. Leimkuhler and S. Reich,
Simulating Hamiltonian Dynamics (Cambridge University Press,Cambridge, 2004).[4] K. Chen and D. P. Landau, Phys. Rev. B , 3266 (1994).[5] A. Bunker, K. Chen, and D. P. Landau, Phys. Rev. B , 9259 (1996).[6] H. G. Evertz and D. P. Landau, Phys. Rev. B , 12302 (1996).[7] B. V. Costa, J. E. R. Costa, and D. P. Landau, J. Appl. Phys. , 5746 (1997).[8] S. Nos´e, J. Chem. Phys. , 511 (1984).[9] W. G. Hoover, J. Chem. Phys. , 1695 (1985).[10] J. Jellinek, J. Phys. Chem. , 3163 (1988).[11] D. J. Evans and G. P. Morriss, Statistical Mechanics of Nonequilibrium Liquids (Academic,New York, 1990).[12] M. E. Tuckerman, C. J. Mundy, and G. J. Martyna, Europhys. Lett. , 149 (1999).[13] M. E. Tuckerman, Y. Liu, G. Ciccotti, and G. J. Martyna, J. Chem. Phys. , 1678 (2001).[14] A. Sergi and M. Ferrario, Phys. Rev. E , Art. No. 056125 (2001).[15] A. Sergi, Phys. Rev. E , Art. No. 021101 (2003).[16] A. Sergi and P. V. Giaquinta, J. Stat. Mech. , P02013 (2007).[17] G. S. Ezra, J. Math. Chem. , 339 (2002).[18] G. S. Ezra, J. Math. Chem. , 29 (2004).[19] G. S. Ezra, J. Chem. Phys. , Art. No. 034104 (2006).[20] A. Bulgac and D. Kusnezov, Ann. Phys. , 187 (1990).[21] J. E. Marsden and T. S. Ratiu, Introduction to Mechanics and Symmetry (Springer-Verlag,New York, 1999).[22] A. Bulgac and D. Kusnezov, Phys. Rev. A (8), 5045 (1990).[23] D. Kusnezov, A. Bulgac, and W. Bauer, Ann. Phys. , 155 (1990).[24] D. Kuznezov and A. Bulgac, Ann. Phys. , 180 (1992).[25] J. Frank, W. Z. Huang, and B. Leimkuhler, J. Comp. Phys. , 160 (1997).
26] T. Arponen and B. Leimkuhler, BIT Numerical Math. , 403 (2004).[27] S. H. Tsai, M. Krech, and D. P. Landau, Braz. J. Phys. , 384 (2004).[28] R. I. McLachlan and D. R. J. O’Neale, J. Phys. A , L447 (2006).[29] W. G. Hoover, C. G. Hoover, and D. J. Isbister, Phys. Rev. E , 026209 (2001).[30] F. Legoll, M. Luskin, and R. Moeckel, Arch. Rat. Mech. Anal. , 449 (2007).[31] F. Legoll, M. Luskin, and R. Moeckel, Nonlinearity , 1673 (2009).[32] G. J. Martyna, M. L. Klein, and M. E. Tuckerman, J. Chem. Phys. , 2635 (1992).[33] J. D. Ramshaw, J. Non-Equil. Thermo. , 33 (1991).[34] A. C. da Silva and A. Weinstein, Geometric Models for Noncommutative Algebra (AMS, NewYork, 1999).[35] B. Schutz,
Geometrical methods of mathematical physics (Cambridge University Press, Cam-bridge, 1980).[36] M. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem. Phys. , 1990 (1992).[37] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. Klein, Mol. Phys. , 1117 (1996).[38] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. Fisher, A. Garg, and W. Zwerger, Rev.Mod. Phys. , 1 (1987).[39] W. G. Hoover, J. Chem. Phys. , 4164 (1998).[40] J. R. Dorfman, An Introduction to Chaos in Nonequilibrium Statistical Mechanics (CambridgeUniversity Press, Cambridge, 1999).[41] C. J. Mundy, S. Balasubramanian, K. Bagchi, M. E. Tuckerman, G. J. Martyna, and M. L.Klein,
Nonequilibrium Molecular Dynamics (Wiley-VCH, New York, 2000), vol. 14 of
Reviewsin Computational Chemistry , pp. 291–397.[42] W. G. Hoover, Molecular Simulation , 13 (2007).[43] W. G. Hoover and C. G. Hoover, J. Chem. Phys. , Art. No. 164113 (2007). H t H BK H BKNH H BKNHC