The Canonical Ensemble via Symplectic Integrators using Nosé and Nosé-Poincaré chains .

Abstract Simulations that sample from the canonical ensemble can be generated by the addition of a single degree of freedom, provided that the system is ergodic, as described by Nosé with subsequent modifications by Hoover to allow sampling in real time. Nosé-Hoover dynamics is not ergodic for small or stiff systems and the addition of auxiliary thermostats is needed to overcome this deficiency. Nosé-Hoover dynamics, like it’s derivatives, does not have a Hamiltonian structure, precluding the use of symplectic integrators which are noted for their long term stability and structure preservation. As an alternative to Nosé-Hoover, the Hamiltonian Nosé-Poincaré method was proposed by Bond, Laird and Leimkuhler [S.D. Bond, B.B. Laird and B.J. Leimkuhler, J. Comp. Phys., 151, 114, (1999)], but the straightforward addition of thermostatting chains does not sample from the canonical ensemble. In this paper a method is proposed whereby additional thermostats can be applied to a Hamiltonian system while retaining sampling from the canonical ensemble. This technique has been used to construct thermostatting chains for the Nosé and Nosé-Poincaré methods.


I. INTRODUCTION
͑a͒ Nose dynamics and its derivatives [1][2][3][4][5][6][7] are popular schemes for implementing molecular simulations at constant temperature.In all such schemes, sampling from the canonical ensemble is conditional on the system being ergodic, and in small or stiff systems this is often not the case.To overcome the lack of ergodicity in these systems further modifications to the Nose ´-Hoover method were proposed by Martyna, Klein, and Tuckerman 8 with new thermostats added to control each previous thermostat to form a thermostatting chain.This method is successful but suffers from the limitations of the Nose ´-Hoover method: there is no Hamiltonian from which it is derived and hence symplectic methods are not applicable.More recently the real time Nose ´-Poincare ḿethod was proposed by Bond, Laird, and Leimkuhler 6 which is based on an extended Hamiltonian, but the application of thermostatting chains to the Nose ´-Poincare ´method in a straightforward manner does not result in sampling from the canonical ensemble.A generalized thermostatting technique has been developed by Leimkuhler and Laird; 9 in this scheme an auxiliary heat bath is coupled into the thermostatting variables, an example of the auxiliary heat bath is a box containing billiards.This requires the design of the auxiliary heat bath and can require long integration time for the correct sampling to occur if the bath is poorly chosen.
This paper introduces the idea of adding multiple thermostats to a Hamiltonian which has been modified by Nose ´'s method, while retaining sampling from the canonical ensemble.Here we employ a regularizing term in the Nose ´or Nose ´-Poincare ´Hamiltonian to ensure bounded integrals over the auxiliary variables.This technique is used to introduce additional terms into the thermostatting chain in both the Nose ´and Nose ´-Poincare ´Hamiltonians, it is then possible to prove analytically that they sample from the canonical ensemble.In addition the fast convergence to the canoni-cal ensemble, which is a characteristic of the Nose ´-Hoover chains, is observed.Before proceeding to describe the new thermostatting technique, we introduce the Nose ´-Poincare formulation and other elements on which our method is based.

A. Nose ´and Nose ´-Hoover schemes
͑b͒ Nose ´'s method 1 adds a thermostatting variable to the equations of motion to act as a heat bath.Given a Hamiltonian system where H(q,p) is the energy of an N-body system, qϭ(q 1 ,q 2 ,...,q N ) and pϭ(p 1 ,p 2 ,...,p N ) are the positions and momenta of the N bodies, Nose ´proposed the extended Hamiltonian, where s is the new thermostatting variable, p s its corresponding momentum, T is temperature, Q the Nose ´mass, and k is the Boltzmann constant.In the equations of motion for this Hamiltonian, the time is scaled by the thermostatting variable s.
Hoover developed the idea of applying a Sundman transformation, dt/dtЈϭs, to correct the dynamics, but this destroys the Hamiltonian structure so that symplectic methods are no longer applicable.Applying the Sundman transformation, and substituting p i Јϭp i /s, tЈϭ͐ dt/s, p s Јϭp s /s, then making the substitution p ϭQ(1/s)ds/dtЈ, ϭln s the equations of motion become This form is now known as the Nose ´-Hoover thermostat. 2

B. Nose ´-Hoover chains
͑c͒ Martyna, Klein, and Tuckerman 8 proposed a method to overcome the lack of ergodicity in small or stiff systems.
Here each thermostat is controlled by another thermostat, forming a thermostat chain.In standard Nose ´-Hoover dynamics the distribution has a Gaussian dependence on the particle momenta, p, as well as the thermostat momentum, p .The Gaussian fluctuations of p are driven by the thermostat but there is nothing to drive the fluctuations of p unless further thermostats are added as described above.The modified dynamics, for M thermostats, can then be expressed as ϪkT ͪ.

͑11͒
These equations can be shown to produce the correct phasespace distributions.

C. Nose ´-Poincare ´method
͑d͒ Neither the Nose ´-Hoover or the Nose ´-Hoover chain methods have a corresponding Hamiltonian which means that symplectic integrators, with their associated long term stability and structure preserving characteristics, are not applicable.The real time Nose ´-Poincare ´method was introduced ͑along with a symplectic integrator͒ by Bond, Laird, and Leimkuhler. 6The reformulation for a Hamiltonian system with energy H(q, p) is ϩNkT ln sϪH 0 ͪ s.

͑12͒
Here N is the number of degrees of freedom of the real system, and H 0 is chosen such that the Nose ´-Poincare Hamiltonian, H NP , is zero when evaluated at the initial conditions.
The Nose ´-Poincare ´method has been extended to NPT ensemble simulation 10 and shown there to be an efficient method.However, it can be shown that the application of thermostatting chains to the Nose ´-Poincare method in a straightforward manner 6 does not sample from the canonical ensemble.

II. MULTIPLE THERMOSTATS
͑e͒ It is possible to introduce additional thermostats into the Nose ´and Nose ´-Poincare ´methods while retaining both their Hamiltonian structure and sampling from the canonical ensemble.This can be illustrated in a more general setting by rewriting the Nose ´method, ͑1͒, to include the momenta of the thermostatting variable with the system momenta such that p ˆϭ( p 1 ,p 2 ,...,p N ,p Nϩ1 ) to give ,p Nϩ1 ͪ .
A second thermostat can be added as follows: where f 2 (s 2 ) is a real valued function, g is a scalar and the new thermostat is applied to M of the momenta, the thermostatted set being ͕p i 1 Ј ,...,p i M Ј ͖ and the nonthermostatted set being ͕p j 1 Ј ,...,p j Nϩ1ϪM Ј ͖ for some integers i 1 ,...,i M , j 1 ,..., j Nϩ1ϪM .Note that the thermostatted set may include any of the system momenta and the thermostatting variable momenta.The partition function for this method, for energy E, is defined as

͑13͒
We can substitute p i Јϭp i /s 1 , 1рiрN, p Nϩ1 Ј ϭp Nϩ1 , the volume element then becomes dp ˆϭs 1 N dp ˆЈ, where p ˆЈ is defined as above.There is no upper limit in momentum space so we can change the order of integration of dp ˆЈ and ds 1 giving Using the equivalence relation for ␦, ␦(g(x))ϭ␦(xϪx 0 )/͉gЈ(x)͉, where x 0 is the zero of g(x)ϭ0, for xϭs 1 , and noting that ͉s 1 ͉ϭs 1 since ln s 1 is not defined for s 1 Ͻ0, we get We can substitute p k Љϭp k Ј/s 2 , k͕i 1 ,...,i M ͖; p k Љϭp k Ј , k͕ j 1 ,..., j Nϩ1ϪM ͖ the volume element then becomes dp ˆЈϭs 2

M dp ˆЉ
where p ˆЉϭ(p 1 Љ ,p 2 Љ ,...,p Nϩ1 Љ ).There is no upper limit in momentum space so we can change the order of integration of dp ˆЉ and ds 2 giving Zϭ 1 If we arrange that gϭM and that Integrating over both thermostat momenta, p Nϩ1 Љ and p Nϩ2 , gives where A similar proof can be applied to the Nose ´-Poincare ḿethod.This process can be repeated to add more thermostats, with the possibility at each stage of thermostatting the previous thermostat's momenta in addition to any of the other momenta.

III. NOSE ´CHAINS AND NOSE ´-POINCARE ´CHAINS
͑f͒ In this section the application of multiple thermostats to both the Nose ´and Nose ´-Poincare ´extended Hamiltonians is considered, in order to generate thermostatting chains.

A. Nose ´chains
͑g͒ Thermostatting chains, consisting of M thermostats, can be added to the Nose ´equation ͑1͒, with some additional terms as follows: where the auxiliary functions ͕ f i (s i )͖ are real valued and satisfy

͑15͒
This extended system produces a canonical ensemble of N ϩM degrees of freedom.The partition function for this ensemble, for energy E, is defined as

͑16͒
We can substitute pЈϭp/s 1 , the volume element then becomes dpϭs 1 N dpЈ.There is no upper limit in momentum space so we can change the order of integration of dpЈ and ds 1 giving the integral over s 1 as Using the equivalence relation for ␦, ␦(g(x))ϭ␦(xϪx 0 )/͉gЈ(x)͉, where x 0 is the zero of g(x)ϭ0, for xϭs 1 , and noting that ͉s 1 ͉ϭs 1 since ln s 1 is not defined for s 1 Ͻ0, we get

͑20͒
Substituting p s 1 Ј ϭp s 1 /s 2 , changing the order of integration, and integrating ͑20͒ over s 2 we get

͑21͒
where K 2 is defined in ͑15͒.Repeating this for s 3 ,...,s M gives

͑22͒
Changing the order of integration and integrating ͑22͒ over all p s i Ј and p s M gives

͑24͒
This means that constant energy dynamics of the extended Hamiltonian H NC correspond to constant temperature dynamics of H(q,p/s 1 ).

B. Nose ´-Poincare ´chains
͑h͒ In a similar manner, thermostatting chains, consisting of M thermostats, can be added to the Nose ´-Poincare ´equation ͑12͒, with some additional terms as follows: where the auxiliary functions ͕ f i (s i )͖ are real valued and satisfy Eq. ͑15͒, and H 0 is Eq.͑14͒ evaluated at the initial conditions.
This extended system produces a canonical ensemble of NϩM degrees of freedom.The partition function for this ensemble is defined as

͑25͒
Substituting pЈϭp/s 1 , the volume element then becomes dpϭs 1 N dpЈ.There is no upper limit in momentum space so we can change the order of integration of dpЈ and ds 1 giving the integral over s 1 as where H(q,pЈ,p ˆs), p ˆs , and F i are defined in Eqs.͑17͒-͑19͒.Using the equivalence relation for ␦, ␦(g(x))ϭ␦(x Ϫx 0 )/͉gЈ(x 0 )͉, where x 0 is the zero of g(x)ϭ0, for xϭs 1 , to get The remaining thermostatting variables can be integrated out as above in Eqs.͑21͒-͑24͒ to get the partition function

C. Auxiliary function
͑i͒ For the Nose ´chains to work correctly an auxiliary function, f i (s i ), must be chosen not only to satisfy Eq. ͑15͒ but to provide a suitable modification to the thermostats.One such choice is where C i , the auxiliary function coefficient, is a constant.The value a i is chosen as the required average value of s i , generally 1, as the additional term will operate as a negative feedback loop to minimize (a i Ϫs i ), as can be seen from the equations of motion.For a Hamiltonian of the form the equations of motion for s i and p s i in the equivalent Nose ´chain system will be If C i is sufficiently small, if s i increases above a i then p s i will decrease, eventually decreasing s i .Conversely, if s i decreases below a i then p s i will increase, eventually increasing s i .

D. Estimation of the auxiliary function coefficient
͑j͒ The value of C i , iу2 can be estimated by considering the equation of motion for the momenta of one of the thermostats, s i , Then the changes in s i are driven by the changes in p s iϪ1 .The purpose of the auxiliary function is to limit the excursions of s i , which can be achieved if ds i /dp s i is a maximum at s i ϭa i .The negative feedback loops arising in Nose ´dynamics drive ͗p ˙si ͘ to zero, in the above equation, over a sufficiently long integration time.For the purpose of estimating the value of C i , we will assume that p ˙si is small.Then, from Eq. ͑28͒,

͑29͒
differentiating with respect to s i , Differentiating again to get the turning points, and substituting s i ϭa i ,

͑31͒
Setting d 2 p s iϪ1 /ds i 2 ϭ0 and solving for C i gives . ͑32͒ Evaluating d 3 p s iϪ1 /ds i 3 at this point gives a positive value, indicating a minimum for dp s iϪ1 /ds i or a maximum for ds i /dp s iϪ1 as required.
Experimental data from tests with the harmonic oscillator, with a i ϭ1, show that Nose ´chains will not work with C i Ͼ1/8kT, but will work for all C i Ͻ1/8kT.However, with very small values of C i the additional thermostats become ineffective as s i is restricted to a value close to 1.

A. Hamiltonian splitting method
͑k͒ The numerical methods used for the following experiments are based on the following general Hamiltonian: The Nose ´-chains method derived from this, with M thermostats based on the auxiliary function in Eq. ͑27͒ is then giving a Nose ´-Poincare ´-chains method, where H 0 is chosen as the initial value of H NC .The equations of motion are where sϭ(s 1 ,...,s M ), p s ϭ(p s 1 ,...,p s M ), and pϭp ˜/s 1 .The thermostats have introduced an implicit coupling into the equations of motion, but an explicit method can be formulated by splitting the Hamiltonian and corresponding Liouville operator.For an odd number of thermostats, M, this can be reduced to three Hamiltonians by employing even-odd splitting of the extended variables.Then if, we have Using a symmetric splitting of the Liouville operator to get a symplectic and time reversible method, iL H ϭ͕.,H͖ϭ͕.,H 1 ͖ϩ͕.,H 2 ͖ϩ͕.,H 3 ͖ϭiL This splitting introduces an error of order ⌬t 3 at each step in terms of the solution operator, giving a second order method, The dynamics for H 1 and H 3 can be solved in a straightforward manner as each s i and p s i are decoupled, leaving H 2 to be solved either analytically or by using the generalized leapfrog algorithm. 6,11 Harmonic oscillator ͑l͒ The harmonic oscillator is generally regarded as one of the hardest models to thermostat and as such is a good test for these methods.The Hamiltonian for the test system is giving the Nose ´-Poincare ´-chains method, To illustrate the importance of the correct selection of thermostatting parameters, and the improvements obtained using thermostatting chains, further experiments were carried out.When very small values are used for the C i , the thermostats are forced to be close to 1, preventing them from operating and producing distributions that would normally be expected from the standard Nose ´or Nose ´-Poincare ´methods, without thermostatting chains, for the Harmonic oscillator.With the parameters the same as the above experiment except ⌬tϭ0.01,C 2 ϭ0.0008,C 3 ϭ0.0004,C 4 ϭ0.0002,C 5 ϭ0.0001 gave the results in Fig. 2.
As we can see from these examples, it is possible to choose thermostat masses to implement an efficient canonical sampling with the chain technique, while using symplectic integrators.

C. Estimating thermostat masses
͑m͒ Values for the thermostatting masses, Q j , can be estimated by simplifying the equations of motion to decouple the thermostats, the resulting equations can then be linearized to evaluate their behavior near to a point of equilibrium.When calculating the masses for Nose ´-Hoover chains 8 it is assumed that adjacent thermostats are slow in comparison to the thermostat of interest, average values can then be used for s jϪ1 , s jϩ1 , p s jϪ1 , p s jϩ1 .Under these conditions, the analysis for the Nose ´-Poincare ´chains is similar to that provided by Nose ´4 provided that we also consider the variation in s 1 to be slow, in order that we can replace it with it's average value, for all but the first thermostat.For Nose ´-Poincare ´chains the equations of motion for s j , p s j for 1 Ͻ jрM are

͑34͒
Rearranging and differentiating ͑33͒, then substituting into ͑34͒, We will consider a fluctuation ␦s j of s j around an average ͗s j ͘, s j ϭ͗s j ͘ϩ␦s j .

͑36͒
Linearizing ͑35͒ to obtain an equation for ␦s, If the change in s j is much faster than the rest of the system, then the change of the momentum can be ignored as the constant temperature is maintained by s j , then, since a j is chosen as ͗s j ͘, as discussed in Sec.III C͑i͒.Sub- stituting ͑38͒ into ͑37͒, expanding the left-hand side and substituting ͗s j ͘ϭa j , ͗s jϩ1 ͘ϭa jϩ1 , s 1 Ϸ͗s 1 ͘, we get ␦s ¨jϭϪ

͑39͒
giving a self-oscillation frequency, w j , of .

͑40͒
Since we normally choose a j ϭ1, j 1, ͑40͒ reduces to the more general form .

͑41͒
For the remaining thermostat's variables, s 1 and p s 1 , the equations of motion are

͑43͒
Following a similar procedure to that above gives a selfoscillation frequency, where a 2 ϭ1, of which is of the form familiar from Nose ´'s paper. 4

D. Optimum thermostat masses
͑n͒ To evaluate the relationship between the selfoscillation frequencies and the optimum choice of the Nose ánd auxiliary masses, experiments were carried out to assess the deviation from the required distribution with varying masses.The experiments were based on the harmonic oscillator model with frequency 1.0, auxiliary function coefficient C 2 ϭ0.08, a 2 ϭ1, and having two thermostats, the original Nose ´thermostat and one auxiliary thermostat.The initial conditions were chosen such that the average value of the Nose ´variable s 1 was 1.0 and the results were taken after 5 000 000 steps, with a step size of 0.005.
In the first experiment the auxiliary thermostat mass was chosen to have a self-oscillation frequency equal to that of the harmonic oscillator and the Nose ´thermostat was varied over a range of values, producing the results in Figs. 3 and 4.
Here Q 1 has been normalized so that 1.0 is the value given by Eq. ͑44͒ and ⌬ D represents the mean square difference between the actual and theoretical distributions.These indicate that the optimum choice of Nose ´mass is near, but less than, its self-oscillation frequency, which is consistent with the results obtained for the auxiliary heat bath method, 9 but with good results for smaller values.
For the second experiment the Nose ´mass was fixed at half that of its self-oscillation frequency and the mass of the auxiliary thermostat was varied, giving the results in Fig. 5.
Here Q 2 has been normalized so that 1.0 is the value given by Eq. ͑41͒ and ⌬ D is defined as above.From these results the optimum value for the auxiliary mass is around its selfoscillation frequency, but good results are obtained over a large range of values.