On ergodicity for multi-dimensional harmonic oscillator systems with Nose-Hoover type thermostat
aa r X i v : . [ n li n . C D ] A ug On ergodicity for multi-dimensional harmonic oscillator systemswith Nose-Hoover type thermostat
Ikuo Fukuda , Kei Moritsugu , and Yoshifumi Fukunishi Institute for Protein Research, Osaka University,3-2 Yamadaoka, Suita, Osaka 565-0871,Japan and Graduate School of Simulation Studies,University of Hyogo, Kobe 650-0047, Japan Graduate School of Medical Life Science,Yokohama City University, Yokohama 230-0045, Japan and Cellular and Molecular Biotechnology Research Institute, AIST (Dated: August 16, 2020)
Abstract
A simple proof and detailed analysis on the non-ergodicity for multidimensional harmonic os-cillator systems with Nose-Hoover type thermostat are given. The origin of the nonergodicity issymmetries in the multidimensional target physical system, and is differ from that in the Nose-Hoover thermostat with the 1-dimensional harmonic oscillator. A new simple deterministic methodto recover the ergodicity is also presented. An individual thermostat variable is attached to eachdegree of freedom, and all these variables act on a friction coefficient for each degree of freedom.This action is linear and controlled by a Nos´e mass matrix Q , which is a matrix analogue of thescalar Nos´e’s mass. Matrix Q can break the symmetry and contribute to attain the ergodicity. . INTRODUCTION The Nos´e-Hoover (NH) [16, 29] equation has been utilized as basic equations of motion(EOM) in molecular dynamics (MD), which is now an important tool to perform a realisticsimulation of a physical system [1, 17, 34]. The NH equation is an ordinary differentialequation (ODE), based on the Newtonian EOM described by physical coordinates x ∈ R n and momenta p ∈ R n . It is obtained by adding a friction force that is − ( ζ /Q ) p to theNewtonian EOM and by adding a EOM for ζ , where Q is a real parameter often called asthe Nos´e’s mass. The friction coefficient-like quantity ζ ∈ R is thus a dynamical variable,and it is introduced to control the temperature of the target physical system described by( x, p ) and maintain the value around a desired value T ex . It is shown that the physicalsystem obeys the Boltzmann-Gibbs (BG), or canonical, distribution at temperature T ex ifthe total system described by ( x, p, ζ ) satisfies the ergodic condition.1-dimensional harmonic oscillator (1HO) has been investigated, for theoretically studyingthe NH equation, as the most simple model system that describes near a physical equilibrium.In a viewpoint of dynamical system study, the NH equation with a 1HO has first been studiednumerically [32] and revealed to include both regular and chaotic motions, implying thatNH EOM with 1HO is nonergodic. The origin of the nonergodicity has been considered as alack of “complexity”, that is, small degrees of freedom of the system (which is three) and asimple form of the ODE involving only two nonlinear terms − ζ p and p . Its non-ergodicityhas been demonstrated theoretically in case of a sufficiently large Q by using KAM theoryto show the existence of invariant tori [23]. In contrast, the NH chain method, which is anextension of the original NH method via introducing multidimensional ζ ∈ R m , has longbeen considered that it gives the ergodicity even in the case of the 1HO. However, recently,Patra and Bhattacharya indicated in long-time numerical simulations that the NHC with1HO is non ergodic [31].However, 1HO is a special model in a viewpoint of the energy density of state Ω( e ) forwhich Ω( e ) =constant, which is not increasing with respect to the energy e of the physicalsystem. In this respect, n -dimensional harmonic oscillator ( n HO) with n > e ) = ce n − ( c is irrelevant to e ) shows the increasing. Despite the importance inthis respect, NH with n HO has not been much investigated. One of a few examples is thestudy by Nos´e [28], where recurrence phenomena strongly depending on initial conditions2ere found and secular periodic modes can be captured by a simple Hamiltonian.In the current study, we discuss the non ergodicity of n HO with n >
1. We specificallyshow in a simple manner that the NH EOM with n HO with n > n HOis identical. This applies not only to the original NH system but also for NH type systems,which includes a number of kinds of EOM such as the NHC equation. The identical conditionmeans that all masses and the spring constants are identical for all the degrees of freedom.With or without this condition, symmetry of the system may differ much and the dynamicsof the system can differ quantitatively [28]. We also demonstrate the non ergodicity of NHtype with n HO with n > x , . . . , x i . . . , x j , . . . , x n ) ( x , . . . , x j . . . , x i , . . . , x n ) becomes one of the symmetries ifthe system is isotropic. Even if the EOM that is considered to be ergodic for 1HO, it shouldbe non ergodic for n HO with n > n HO with n >
1. These issues will be clarified in the current study.We further propose a method to recover the ergodicity even for the idealized model of n HO with n >
1. For this, we first modify the density dynamics (DD) scheme [13], whichhas been developed to produce an arbitrary phase-space density. The DD is described by( x, p, ζ ) where ζ is an additional scalar variable as in the NH EOM and plays a role to controlthe dynamics to yield the arbitrarily given phase-space density. We then obtain a new EOM(we call this splitting DD), which utilizes a vertorized ζ and splits the role of the originaladditional variable to directly act on each degree of freedom. In general, the splitting DDcan be free from the symmetry discussed above. Next, we apply the BG density, for thephase-space density, into the splitting DD to produce the BG distribution and have a newNH form EOM, which we call splitting NH EOM. This new EOM can also be free from the3ymmetry through the splitting feature, even if K ( p ) and U ( x ), the source of the symmetry,are introduced.One of the key issues to avoid the problem originated from the symmetry and to reachthe ergodicity in the splitting NH is an extension of the original Nos´e mass parameter Q .Nos´e mass is not a scalar but a matrix. This extension can be viewed as natural in that thedensity of ζ is based on a quadratic form. Such a matrix approach would be simple, andonly a few discussions have been done [33]. The simplicity of this idea might be the reasonthat it is has not paid much attention. However, we show that this idea gives a higher costperformance than supposed to avoid the nonergodicity.The splitting NH should not be restricted in the applications to the n HO model system.As well as short-rage vibrational interactions, many part of long range interactions of parti-cles in classical physical system can be approximated by harmonic interactions around theirequilibriums, which can then be mimicked by n HO interactions. More directly, e.g., repre-sentations of the physical system as a set of harmonic oscillators in normal mode study welldescribe the feature of a biomolecular system. It is discussed that perturbations by a chemi-cal reaction via enzyme or by a docking of a medicinal molecule excite normal modes [3, 36].These excited normal modes allow an oscillator approximation such as elastic network model[35]. Conversely, intramolecular vibrational energy can be transferred from a given normalmode [27].Thus, the ability to enhance the phase-space sampling and recover the ergodicity for the n HO with n > n HO with n = 2 , n HO model with n >
I. EQUATIONS OF MOTION
Our target dynamical system can be represented as the following ODE:˙ x = M − p ∈ R n , ˙ p = F ( x ) + λ ( ω ) p ∈ R n , ˙ ζ = Λ( ω ) ∈ R m , (1)where x ≡ ( x , . . . , x n ) ∈ D ⊂ R n and p ≡ ( p , . . . , p n ) ∈ R n are atomic coordinates andmomenta of a physical system of n degrees of freedom; F represents a force, which is a C function defined on a domain D ; M represents the mass parameters, which is a symmetric,positive-definite square matrix of size n over R . Along with these quantities associatedwith Newtonian EOM, ζ ∈ R m is an additional dynamical variable, relating to a notion offrictional coefficient or thermostat. Thus the phase space is Ω := D × R n × R m ⊂ R N with N ≡ n + m , and the phase-space point is represented as ω ≡ ( x, p, ζ ) ∈ Ω. To the physicalsystem, a C function λ : Ω → R provides − λ ( ω ), which can be viewed as a dynamicalfrictional “coefficient” and essentially depends on the additional variable ζ ∈ R m . Thetime development of ζ is described by Λ : Ω → R m , which is of class C . The functions λ and Λ may contain potential energy U ( x ) ∈ R , wherein F = −∇ U , and kinetic energy K ( p ) ≡ ( p | M − p )= P ni,j =1 M − ij p i p j of the physical system .We give several examples that fall into the form of EOM (1). Example 1 λ ( ω ) ≡ ζ /Q ∈ R and Λ( ω ) ≡ K ( p ) − nk B T ex with m ≡ , where Q > is aparameter (often called as Nos´e’s mass) give the the NH equation [16, 29]: ˙ x = M − p ∈ R n , ˙ p = F ( x ) − ( ζ /Q ) p ∈ R n , ˙ ζ = 2 K ( p ) − nk B T ex ∈ R . (2) This is introduced to control the physical system temperature K ( p ) /nk B ( k B is Boltzmann’sconstant) into a target temperature T ex > and can yield the canonical ensemble. Example 2 λ ( ω ) ≡ ζ /Q ∈ R and Λ( ω ) ≡ ( G ( ω ) − ζ ζ /Q , . . . , G j ( ω ) − ζ j ζ j +1 /Q j +1 , . . . , G m ( ω )) , with G ( ω ) ≡ K ( p ) − nk B T ex ∈ R ,G j ( ω ) ≡ ζ j − /Q j − − k B T ex ∈ R , j = 2 , . . . , m, here Q . . . , Q m > are parameters, give the NH chain (NHC) equation [26]: ˙ x = M − p ∈ R n , ˙ p = F ( x ) − ( ζ /Q ) p ∈ R n , ˙ ζ j = G j ( ω ) − ζ j ζ j +1 /Q j +1 ∈ R , j = 1 , . . . , m − , ˙ ζ m = G m ( ω ) ∈ R . This can be viewed as an extended form of the NH EOM.
Example 3
The kinetic moments method [15, 19], which can be represented by ˙ x = M − p ∈ R n , ˙ p = F ( x ) − ( ζ /Q + ˆ K ( p ) ζ /Q ) p ∈ R n , ˙ ζ = ˆ K ( p ) − ∈ R , ˙ ζ = ˆ K ( p )( ˆ K ( p ) − ( n + 2) /n ) ∈ R , where ˆ K ( p ) ≡ K ( p ) /nk B T ex , becomes an example of (1). Example 4
The generalized Gaussian moment thermostatting method [25] represented by ˙ x = M − p ∈ R n , ˙ p = F ( x ) − m P j =1 j P k =1 a k − (2 K ( p )) k − ( k B T ex ) j − k ζ j /Q j ! p ∈ R n , ˙ ζ j = a j − (2 K ( p )) j − n ( k B T ex ) j ∈ R , j = 1 , . . . , m, with a j ≡ Q jk =1 ( n + 2 k ) − , becomes also an example. There can be found more examples in thermostat methods; see e.g., Refs. [2, 4, 5, 7,14, 18, 20–22, 30, 33] for details on thermostat methods and their development. Otherexamples include e.g., the coupled NH equations of motion, which is introduced to fluctuatethe temperature of the heat bath for the physical system [9, 11].
III. NON ERGODICITY FOR ISOTROPIC OSCILLATOR SYSTEM
The target ODE (1) can be represented as˙ ω = X ( ω ) , (4)6here X becomes a C vector field defined on a domain Ω of R N . Assuming the completenessof X , we let { T t } ≡ { T t : Ω → Ω | t ∈ R } be the flow generated by the field X . Considerthe case where we have an invariant measure µ of the flow { T t } : ∀ t ∈ R , ∀ A ∈ L Ω N , µ ( T − t A ) = µ ( A ) , where L Ω N ≡ L N ∩ Ω with L N being the Lebesgue measurable sets on R N . We assume that0 < µ (Ω) < ∞ and that µ ∼ l N , (5)i.e., µ and the Lebesgue measure l N of R N are absolutely continuous each other. For example,on Examples 1 and 2, we have an invariant measure defined by L Ω N → [0 , ∞ ) , A µ ( A ) := Z A ρdl N (6)with a (strictly positive and measurable) density ρ : Ω → R , satisfying the Liouville equa-tion [6, 10] div ρX = 0 . (7)A subset A ∈ L Ω N is said to be an invariant set if T − t ( A ) = A for all t ∈ R . The ergodicityfor the measure space (Ω , L Ω N , µ ) with the flow { T t } holds if any invariant set A is trivial,i.e., µ ( A ) = 0 or µ (Ω \ A ) = 0. In other word, if we have an invariant set A such that µ ( A ) > µ (Ω \ A ) > , (8)then the ergodicity does not hold.We show { T t } is not ergodic for a harmonic oscillator system with n >
1. The conditionfor n > n = 1 is demonstrated in Ref. [23]. Here, an isotropic condition in a harmonicoscillator system in Eq. (1) is described by M = m n and F ( x ) = − kx ∈ R n , (9)where n is the unit matrix of size n , and m and k are strictly positive parameters (repre-senting a mass and spring constant, respectively); i.e., we have˙ x = m − p ∈ R n , ˙ p = − kx + λ ( ω ) p ∈ R n , ˙ ζ = Λ( ω ) ∈ R m . (10)7efining a map γ : Ω → T ( R n ) ∼ = R n by γ ( ω ) := x ∧ p = 12 ( x ⊗ p − p ⊗ x ) , (11)we show Lemma 5 Υ ± ij ≡ { ω ∈ Ω | γ ij ( ω ) ≷ } is an invariant space for the flow of Eq. (10) for i, j = 1 , . . . , n . Proof.
For any solution of ODE (10), ϕ : R ⊃ J → Ω , t ϕ ( t ) ≡ ( x ( t ) , p ( t ) , ζ ( t )) , wehave D ( γ ◦ ϕ )( t ) = Dx ( t ) ∧ p ( t ) + x ( t ) ∧ Dp ( t )= m − p ( t ) ∧ p ( t ) + x ( t ) ∧ ( − kx ( t ) + λ ( ϕ ( t )) p ( t ))= x ( t ) ∧ λ ( ϕ ( t )) p ( t )= λ ( ϕ ( t ))( γ ◦ ϕ )( t )for all t ∈ J , which is an open interval (that may be R ) involving 0. Thus( γ ◦ ϕ )( t ) = exp (cid:18)Z t λ ( ϕ ( s )) ds (cid:19) ( γ ◦ ϕ )(0) ∈ R n for any t ∈ J , implying that ( γ ◦ ϕ )(0) = 0 reads as ( γ ◦ ϕ )( t ) = 0 for all t . HenceΥ ij ≡ { ω ∈ Ω | γ ij ( ω ) = 0 } is an invariant space for i, j = 1 , . . . , n . For any i and j , the continuity of γ ij indicates thatΥ + ij = { γ ij > } and Υ − ij = { γ ij < } are also invariant.Thus, we get Proposition 6
The flow { T t } of ODE (10) is not ergodic with respect to (Ω , L Ω N , µ ) . Proof.
Choose any i , j ∈ { , . . . , n } such that i = j (recall n > ± ij = ∅ .From assumption (5) and the fact that Υ ± ij becomes a nonempty open set of R N , we see µ (Υ + ij ) > µ (Ω \ Υ + ij ) ≥ µ (Υ − ij ) > n >
1, by showing the existence of an invariant set Υ + ij that has a desired property (for8his purpose, finding one such a set is relevant but the existence of many subsets such as ∪ i,j Υ ± ij and ∩ i,j Υ ± ij is less important). However, this discussion lacks explanations why andhow the invariant set Υ ± ij arises. In Appendix, we show that symmetry group O ( N ) acts onODE (10) and its invariant set splits the total phase space Ω producing invariant sets thatcorrespond to Υ ± ij . It should be noted that finding just one solution (or countably manysolutions) confined in a certain subset B does not necessarily indicate the non ergodicity,since B may be a null set.The above discussion on the isotropic harmonic oscillator case can be generalized to aharmonic oscillator with F ( x ) = − K x ∈ R n with K ∈ End R n under the condition that M − K is symmetric, positive definite, and degenerate, i.e., there exist eigenvalues suchthat λ i = λ j for i = j . Proposition 6 in this case is proven by observing that ˜Υ ij ≡{ ω ∈ Ω | ˜ γ ij ( ω ) = 0 } becomes an invariant space, where ˜ γ ( ω ) ≡ V Gx ∧ G M − p with V ≡ G M − K G − being the diagonal matrix for a certain G ∈ GL( n ), and by observing that˜Υ ± ij ≡ { ω ∈ Ω | ˜ γ ij ( ω ) ≷ } becomes a nonempty open invariant set. IV. SPLITTING DENSITY DYNAMICS
To attain the ergodicity for symmetric systems such as the n HO, we propose a generalizedversion of the density dynamics (DD) [13]. The original version of the DD is defined by˙ x i = D p i Θ( ω ) ∈ R , i = 1 , . . . , n, ˙ p i = − D x i Θ( ω ) − D ζ Θ( ω ) p i ∈ R , i = 1 , . . . , n, ˙ ζ = n P i =1 D p i Θ( ω ) p i − nβ − ∈ R , (12)with Θ = − β − ln ρ , where ρ : Ω → R is an arbitrarily given density function, i.e., ρ is afunction that is of class C , strictly positive, and integrable. ζ ∈ R is a dynamical variableand β > ρ , and it involves the NH EOM (2), viz., the NH is recovered ifwe set ρ = ρ NH , where ρ NH ( ω ) ≡ exp (cid:20) − β (cid:18) U ( x ) + K ( p ) + 12 Q ζ (cid:19)(cid:21) (13)with β = 1 /k B T ex .Our generalization for Eq. (12) is based on (i) an extension of the additional scalarvariable ζ ∈ R to a vector variable ζ ≡ ( ζ , . . . , ζ n ) ∈ R n , and (ii) ζ i ’s EOM that is a9atural decomposition of the third equation of (12). Namely, a generalized DD, which wecall splitting density dynamics, is˙ x i = D p i Θ( ω ) ∈ R , i = 1 , . . . , n, ˙ p i = − D x i Θ( ω ) − D ζ i Θ( ω ) p i ∈ R , i = 1 , . . . , n, ˙ ζ i = D p i Θ( ω ) p i − β − ∈ R , i = 1 , . . . , n. (14)As is easily confirmed that the Liouville equation (7) holds for any density function ρ , thisnew EOM can be replaced with Eq. (12). That is, for any P -integrable function g on phasespace Ω, ¯ g := ∃ lim τ →∞ τ Z τ g ( T t ( ω )) dt = Z Ω gρ dl N (cid:30)Z Ω ρdl N =: h g i ∈ R (15)holds with respect to a P -almost every initial point ω , if the flow { T t } is ergodic with respectto the measure P ≡ ρdl N .Some remarks are made [13]. First, fixed points for X, which can be obstructions to theergodicity, do not exist, as long as β >
0. Second, div X = 0 holds (otherwise ρ becomes aninvariant function and should not be almost everywhere constant, so the system does notbecome ergodic), as long as n X i =1 D ζ i ρ = 0 (16)(not identically zero). This is because div X = − P ni =1 D ζ i Θ. Condition (16) is valid inmany cases (see the case later).
V. SPLITTING NOS´E-HOOVER METHOD UTILIZING NOS´E-MASS MATRIX
The meaning of the generalization of (i) and (ii) in Section IV will be clearer when weconsider the NH limit. Here, the NH limit is obtained if we set ρ = ˜ ρ BG , where˜ ρ BG ( ω ) ≡ exp [ − β ( U ( x ) + K ( p ) + K z ( ζ ))] (17)with K z ( ζ ) ≡
12 ( ζ | Q − ζ ) = 12 n X i,j =1 Q − ij ζ i ζ j (18)10eing a quantity corresponding to the kinetic energy for ζ , and β = 1 /k B T ex . The differencefrom Eq. (13) is to utilize, as well as the vectorized ζ ∈ R n , a matrix form of Q , which isa natural extension of the original scalar Nos´e’s mass Q (recovered when n = 1, of course).Specifically, Q ≡ (Q ij ) ∈ End R n should be symmetric and positive definite: we should set itsymmetric without loss of generality, considering that the kinetic energy is a quadratic form;and the positive-definite condition is a natural extension of Q > ρ . Now, applying Eq. (17), the splitting DD (14)turns out to be ˙ x i = ( M − p ) i ∈ R , i = 1 , . . . , n, (19a)˙ p i = F i ( x ) − τ i ( ζ ) p i ∈ R , i = 1 , . . . , n, (19b)˙ ζ i = 2 K i ( p ) − k B T ex ∈ R , i = 1 , . . . , n, (19c)where τ i ( ζ ) ≡ − k B T ex D ζ i ln ˜ ρ BG ( ω )= D i K z ( ζ )= ( Q − ζ ) i = n X j =1 Q − ij ζ j (20)and 2 K i ( p ) ≡ ( M − p ) i p i . (21)Thus, the dynamical frictional “coefficient” τ i ( ζ ) depends on not only one component for ζ (as in the NH and NHC) but also all components ζ , . . . , ζ n or at least two components ζ k , ζ l when we choose Q as a non-diagonal matrix. This is a motivation of above (i), and this“mixing” of ζ components will play a part for avoiding the nonergodicity, as detailed below.Briefly speaking, the fact that the contribution of τ i ( ζ ) to ˙ p i can be different for each i iseffective to break down the isotropic symmetry if it exists in the system. Equation (19c)intends the law of equipartition, i.e., the expected equilibrium condition, ˙ ζ i ∼
0, shouldcontribute to the exact relationship of the law 2 K i = k B T ex for every degree of freedom i ,11hich is validated by 2 K i = h K i i≡ Z Ω ( M − p ) i p i ˜ ρ BG ( ω ) dl N ( ω ) (cid:30)Z Ω ˜ ρ BG dl N = k B T ex (22)where the first equation is owing to (15) under the ergodic condition.Some remarks are made. First note that the EOM (19) is still physically natural in thesense that the first equation is exactly same as that in the Newtonian EOM, and the i thcomponent of the frictional force in Eq. (19b) is proportional to p i and take a form − τ i ( ζ ) p i using a scalar quantity τ i ( ζ ) ∈ R , which conforms to the conventional form for the classicaldynamics treatment. Second, EOM (19) can be viewed as a generalization of the originalNH. This is because, by setting Q − = Q − ≡ Q − · · · · · · ∈ End R n (23)with a scalar input Q − > ζ ≡ P ni =1 ζ i ∈ R , which plays the same role in the original NH variable. In this sense, we willcall Eq. (19) splitting Nos´e-Hoover EOM. Third, as can also be seen from Eq. (22), ODE (19)generates the BG distribution at temperature T ex for physical quantities under the ergodiccondition. This is clearly seen by separately rewriting ˜ ρ BG such that˜ ρ BG ( ω ) = ρ BG ( x, p ) ρ Z ( ζ ) ,ρ BG ( x, p ) ≡ exp (cid:20) − k B T ex ( U ( x ) + K ( p )) (cid:21) ,ρ Z ( ζ ) ≡ exp (cid:20) − k B T ex K z ( ζ ) (cid:21) , Eq. (15) implies the relation ¯ f = h f i BG with respect to a physical quantity f : D × R n → R f ≡ lim τ →∞ τ Z τ f ( x ( t ) , p ( t )) dt = Z Ω f ( x, p ) ˜ ρ BG ( ω ) dl N (cid:30)Z Ω ˜ ρ BG dl N ∈ R = R D × R n f ( x, p ) ρ BG ( x, p ) dxdp R D × R n ρ BG ( x, p ) dxdp =: h f i BG . (24)Finally, note that Condition (16) is valid for ρ = ˜ ρ BG since P ni =1 D ζ i Θ ( ζ ) = P ni,j =1 Q − ij ζ j .Note that matrix form of Q has been utilized before in literature. Samoletov et al. [33]have used it in their development of configurational thermostats, which are thermostatequations in configuration space. Their matrix Q is of size of 3, which corresponds to theadditional 3-vector introduced for their purpose in order to control x instead of p . Theyalso utilized a diagonal form (uncoupled case) for Q in considering physical necessities,though admitted the possibility of coupled case. In these respects, their approach andours are different. Leimkuhler et al. [24] has considered a NH type EOM with introducingrandom noise to improve ergodicity. The frictional term they treated is of form of λ ( ω ) = ζ n + M S ( t, ζ ) ∈ End R n , where ζ is a scalar as in the original NH, and S ( t, ζ ) is ananti-symmetric matrix depending on ζ . Thus, it is different from our term (20). Theydemonstrated the ergodicity for their stochastic dynamics with harmonic oscillators, relatingas a counterpart to our statement of the nonergodicity for ODE. A. On the choice of matrix Q
We describe how the matrix Q defines the the distribution of ζ and how we should set Q for effectively realize the ergodicity. Q determines ζ ’s distribution In contrast to the (marginal) distribution of ( x, p ) described by the RHS of Eq. (24),which is the BG distribution, the distribution of ζ is characterized by the matrix Q and is13escribed by P ζ ≡ P π − ζ : B n → R such that B P ( π − ζ ( B ))= Z D × R n × B ˜ ρ BG dl N (cid:30)Z Ω ˜ ρ BG dl N = Z B ρ Z ( ζ ) dl n ( ζ ) (cid:30)Z R n ρ Z dl n = N z Z B exp (cid:20) − k B T ex ( ζ | Q − ζ ) (cid:21) dζ , where N z ≡ [(2 πk B T ex ) n det Q ] − / . Namely, ζ ∈ R n is distributed ellipsoidally around theorigin. Note that, instead of directly using P ζ , it is often connivent to use the distributionof principal component y ≡ O − ζ ∈ R n for which Q − is diagonalized as O − Q − O =diag( λ , . . . , λ n ) with λ i being a strictly positive eigenvalue of Q − . The distribution of y isgiven by P Y ≡ P ( G − ◦ π ζ ) − : B n → R where B P ( π − ζ ( G ( B )))= N z Z B exp " − k B T ex n X i =1 λ i y i dy, (25)which is the joint distribution of 1-dimensional Gaussian distributions exp h − k B T ex λ i y i i dy i , i = 1 , . . . , n . Note that the distribution of ζ , or y = O − ζ , is not used in obtaining physicalinformation, such as the long-time average of physical variable in Eq. (24), but the explicitform of P Y can be utilized to monitor the convergence of the distribution generated by theflow and numerically judge the ergodicity.
2. We determine Q Here we discuss how we determine Q or Q − . Its overall amplitude can be set by a scalefactor as in the case of the scalar Q , as in the original NH [30]. Thus we should determinethe difference between the matrix elements it in a finer manner. Our criteria for setting thematrix Q − are as follows:(i) it is symmetric and positive definite;(ii) its eigenvalues and eigen vectors are explicitly obtained;14iii) it should not be diagonal;(iv) its diagonal components are nevertheless sufficiently larger than off-diagonal compo-nents;(v) randomness can be easily introduced in the elements.The reason of these requirements is as follows: (i) has been already assumed and thenecessity has also been discussed. (ii) is required to explicitly obtain the distribution of ζ or y ≡ G − ζ . (iii) is needed to enhance the mixing of different components ζ , . . . , and ζ n through the friction term − τ i ( ζ ) p i = − (cid:0) Q − i ζ + · · · + Q − ii ζ i + · · · Q − in ζ n (cid:1) p i (26)in Eq. (19b). Otherwise, Eq. (19b) turns out to be the same form as that of the originalNH, leading to the nonergodicity in the case of the n HO as discussed in Section III. (iv) ifQ − ii is small, then the contribution of ζ i derived by Eq. (19c) will not be much assessed, sothat the equipartition (22) will not be enhanced at least in a relatively short time scale. (v)is needed to break isotropy or symmetries in the target physical system. It is also usefulto emphasize the difference between the splitting NH and the original NH, where the lattercan be characterized as a uniform matrix Q − seen in (23).Using the fact that a symmetric matrix W ∈ End R n is positive definite if and only if ∃ O ∈ O ( n ), ∃ d , . . . , d n > W = O diag( d , . . . , d n ) T O , and using a representation of thegroup O ( n ), we propose the following procedures (1)–(4) for setting Q − :(1) Choose values randomly for θ k,j ∈ ]0 , ǫ [ with 0 < ǫ ≪ π for 1 ≤ j < k ≤ n ,(2) define O := h n h n − · · · h for which h k := r ( θ k, ) r ( θ k, ) · · · r k − ( θ k,k − ), where r i ( θ ) ≡ i − u ( θ )
00 0 1 n − i − ∈ End R n with i being the unit matrix of size i and u ( θ ) ≡ cos θ sin θ − sin θ cos θ ; for convenience, any q ∈ O ( n ), such as q ≡ diag(1 , . . . , , −
1) or interchange matrices q ≡ J ij can be insertedamong the products of h k in defining O , 153) set d i = 1 + δ i with − δ < δ i < δ < i = 1 , . . . , n , where δ i = δ j for i = j , andthen, using a scale factor λ , put D ≡ λ diag( d , . . . , d n ) ≡ : diag( λ , . . . , λ n ), and finally,(4) define Q − := O D T O. (27)Then, condition (i) holds, and (ii) is clear since the eigen values of Q − are λ , . . . , λ n and O i ≡ T ( O i , . . . , O ni ) ∈ R n is obtained to be the eigen vector corresponding to λ i for i = 1 , . . . , n . Condition (iii) will hold, or reset some values of θ k,j if needed. A small ǫ thatbecomes the threshold of θ k,j is useful to contribute to the purpose (iv) in that each r i ( θ ) isnear the identity matrix. Randomness can be introduced through the n ( n − / θ k,j for the sake of (v). Example 7
For n = 2 , we have O = h = r ( θ , ) = u ( θ , ) . For n = 3 , we have O = h h = r ( θ , ) r ( θ , ) r ( θ , )= cos θ , sin θ , − sin θ , cos θ ,
00 0 1 θ , sin θ , − sin θ , cos θ , cos θ , sin θ , − sin θ , cos θ ,
00 0 1 . Instead, we can use e.g., O = J r ( θ , ) J J r ( θ , ) J J r ( θ , ) J , viz., O = cos θ , − sin θ , θ , cos θ ,
00 0 1 cos θ , θ , − sin θ , θ , θ , − sin θ , θ , cos θ , , (28) which indicates the composition of rotations in R around the x,y,z-axis with angles θ , , θ , ,and θ , , respectively. VI. NUMERICS
We numerically tested our considerations, nonergodic property for the conventionalschemes and the ergodic property for the current scheme, using the isotropic harmonicoscillator system defined by (9) with n >
1. We set both the mass m and spring constant k to be 1, and put k B T ex = 1 (all quantities were treated as dimensionless). Numericalintegrations of ODEs were done by the explicit second order scheme described in Ref. [8] for160 time steps with a unit time of h = 10 − , and the numerical errors were checked to bewithin a tolerance within the scheme of the extended system [12].We show that conventional method employing an EOM of the form of Eq. (10) fail inthe ergodic sampling, as stated in Proposition 6. As a conventional method, we have usedthe NHC method (Example 2) with the chain length m = 2 and masses Q = Q = 1.A first case we show is that with n = 2, where we used the initial value x (0) = (0 , ,p (0) = (1 , , and ζ (0) = (0 , x , p ) and theirmarginal distributions. The trajectory shows a hall in a vicinity of the origin, and thisclearly affects the distribution of x . The distribution of p is also weird. Due to the aspecial setting of the initial condition, x ( t ) = x ( t ) and p ( t ) = p ( t ) for all time t , sothat the trajectory of ( x , p ) and the distributions are totally the same as that for ( x , p ) , respectively. In terms of the symmetry, this special initial condition obeys a symmetry ofthe interchange, S = ∈ O (2), so that the initial value ω ≡ ( x (0) , p (0) , ζ (0)) falls inan invariant set A = A R = A O (2) × , defined in Eq. (30) and utilized in the decomposition (29)(see Appendix). Thus, the solution always falls in the invariant set A , indicating the fact that x ( t ) = x ( t ) and p ( t ) = p ( t ) for all t . This initial condition, however, seems too specialand very severe. Thus, we also treated another condition such that x (0) = x (0) = 0 ,p (0) = 1 , p (0) = 2 (with ζ (0) = (0 , S = − ∈ O (2), so that thissolution is also confined in the invariant space, ω ( t ) ∈ A for all t . Although it does not holdthat x ( t ) = x ( t ) and p ( t ) = p ( t ) for all t , this confinement severely effects the motion,and the trajectories of x, p and their distributions exhibited similar nonergodic behavior asthat in Figure 1 (not shown).In the second case, we changed the initial condition only and set as x (0) = 1 , x (0) = 0 ,p (0) = 0 , p (0) = 0 .
01, which does not have any symmetry in O (2) × . The initial value ω is in A [+] because γ ( ω ) = ( x (0) p (0) − p (0) x (0)) > A but confined in A [+] , whichis “large”. However, as shown in Figure 2, the distributions are far from the theoreticalGaussian distributions and the trajectories are biased, suggesting a certain structure. Thethird case we studied is the case with x (0) = 1 , x (0) = 0 , p (0) = 0 , p (0) = −
2, which17lso has no symmetry in O (2) × and ω ∈ A [ − ] (because γ ( ω ) < x, p ) trajectories andunignorable errors were admitted in the distributions, it might be sufficient in practicalsimulations. However, a clear numerical evidence for nonergodicity is a definiteness of thesignature of γ ( ω ( t )). It should be a null occurrence that γ ( ω ( t )) = 0 for all t , if theflow is ergodic. Furthermore, there should not be the case where γ ( ω ( t )) > t or γ ( ω ( t )) < t . Otherwise, it breaks the ergodicity and contradicts the BGdistribution. In fact, its average should be zero under the BG distribution: ¯ γ ij = h γ ij i = h γ ij i BG = 0 if the flow is ergodic with respect to exp [ − βU ( x ) + K ( p )] ρ Z ( ζ ) dω for anysmooth, positive, integrable ρ Z (as long as R D x k exp [ − βU ( x )] dx are finite for k = i and j ). Figure 3 shows γ ( ω ( t )) for the three cases above. The first case (Fig. 3a) correspondsto the null case γ ( ω ( t )) = 0 for all t , and the second (Fig. 3b) and third (Fig. 3c) casescorrespond to γ ( ω ( t )) > γ ( ω ( t )) < t , respectively. These results show thatthe conventional method sampled the phase space in a nonergodic manner. Note also thatthe magnitude of γ ( ω ( t )) in the third case (Fig. 3c) is larger than the second case (Fig. 3b).This result may be the reason why the third case shows relatively good sampling; namely,trajectories staying near the invariant set A shows bad sampling, whereas trajectories thatcan be away from A relatively show good (but not exact) sampling. These staying featuresnear A may suggest a kind of stability of the invariant set A .We tested the splitting NH EOM (19), currently provided scheme, using the same har-monic oscillator system as above. A first example is the case with n = 2, where the initialvalue is the same as the most stiff case used above, viz., x (0) = (0 , , p (0) = (1 , ,ζ (0) = (0 , Q − in the manner stated in Section V A 2, where δ = 0, δ = 0 . λ = 10, and θ , = 0 . x , x , p , p , ζ , and ζ are shown in Figure 4. They are well sampled in the phase space. The distributions agreedthe theoretical distribution, and the errors were sufficiently small, for which we have usedvariable y , instead of ζ , as indicated in Eq. (25). We also observed that γ ( ω ( t )) does notindicate the positive/negative definiteness, as in the conventional method, and rapidly con-verged to the theoretical value 0. We had similar results for other initial conditions. Next,we show the results for other setting of Q − , where δ = 0 . θ , = 0 .
8, while the otherconditions are the same as above. This is a setting where ellipsoid distributed feature for ζ
18s emphasized. We observe in Figure 5 that ζ were distributed ellipsoidally around the originand sampled correctively, as indicated in the theoretical contours and the distributions for y . The next example for the splitting NH EOM (19) is the case with n = 3. In the proceduresfor setting Q − , we put θ , = θ , = θ , = 0 . δ = − . δ = 0, δ = 0 .
2, and λ = 10,and utilized Eq. (28). Initial values were x i (0) = 0 , p i (0) = 1 , ζ i (0) = 0 for all i = 1 , , n = 2). As shown in Figure 6, the scatter plots indicate the ergodicsampling, and the distribution for each variable x i , p i , ζ i for i = 1 , , θ k,j were not set randomly. On the Basis of these results, weconclude that the current method accurately corresponds to the ergodicity. Appendix
We say that a linear symmetry S : R n → R n acts on ODE (1) if it satisfies the following: Definition 8 S ∈ End R n preserves the functions λ and Λ and the domain D such that λ ( S ( x ) , S ( p ) , ζ ) = λ ( x, p, ζ ) and Λ ( S ( x ) , S ( p ) , ζ ) = Λ ( x, p, ζ ) hold for all ( x, p, ζ ) ∈ Ω and S ( D ) ⊂ D . Commutativities also hold: M − ◦ S = S ◦ M − and F ◦ S = S ◦ F . We denote by S the set of all S ∈ End R n that acts on ODE (1). Lemma 9
For any S ∈ S , we have S ( x ( t )) = x ( t ) and S ( p ( t )) = p ( t ) for all t in an interval J ⊂ R , if ϕ : J → Ω , t ( x ( t ) , p ( t ) , ζ ( t )) is a solution of ODE (1) with an initial conditionsatisfying S ( x (0)) = x (0) and S ( p (0)) = p (0) . Proof.
It follows from Definition 8 that ˆ ϕ : R ⊃ J → Ω , t d ( S ( x ( t )) , S ( p ( t )) , ζ ( t )) alsobecomes a solution of the C ODE and ˆ ϕ (0) = ϕ (0) holds. Thus, the uniqueness of theinitial value problem ensures ˆ ϕ = ϕ , so that S ( x ( t )) = x ( t ) and S ( p ( t )) = p ( t ) for all t ∈ J .We thus have for every S ∈ S an invariant set,Ω S := Γ Ds × Γ s × R m , with Γ s ≡ { p ∈ R n | S ( p ) = p } Ds ≡ Γ s ∩ D ≡ { x ∈ D | S ( x ) = x } , indicating that the symmetry ( S ( x ) , S ( p )) = ( x, p )is kept in the dynamics or compatible with the ODE. To show the nonergodic condition (8),we take an approach that is to find an invariant set A whether it meets condition (8) itselfor it separates the total phase space into three invariant sets,Ω = A ⊔ A [+] ⊔ A [ − ] , (29)wherein A [+] meets condition (8). For this, A should be ”large” (for otherwise situation inchoosing A = Ω S , there is the extremely small case Γ s = ∅ or a case of a low dimensionalsubspace). Our target for A is thus an invariant set that are summed up these Ω S in acertain manner: A R ≡ [ S ∈R Ω S = [ S ∈R (Γ Ds × Γ s ) × R m , (30)where R is a certain subset of S such that it is sufficiently large but not too large. Forthe latter condition, for example, we should remove the case where S is the identity id R n ,otherwise A S becomes ”too large” ( S = id R n provides Ω id R n = Ω and so yields Ω \ A R = ∅ ,which does not contribute to the nonergodic condition (8) for A ≡ A R ).We will show that (29) holds with A [ ± ] ≡ Υ ± ij if A ≡ A R , in a special case of the isotropic n HO. Here Υ ± ij are defined in Lemma 5, and this fact can explain the route why Υ ± ij arise.That is, they arise as a complementary set to a sum, in a certain range R , of the invariantset Ω S based on the symmetry S that acts on the ODE. To show the issue, we restrict thecondition such that the dependence of x, p in the functions λ and Λ is only through thepotential and kinetic energies; viz., Condition 10
There exist C functions ˜ λ, ˜Λ : R × R × R m → R such that λ ( x, p, ζ ) =˜ λ ( U ( x ) , K ( p ) , ζ ) and Λ ( x, p, ζ ) = ˜Λ ( U ( x ) , K ( p ) , ζ ) for all ( x, p, ζ ) ∈ ΩThis is not a special condition and are satisfied by Examples 1–4. For the harmonicoscillator system described by Eq. (9) under condition 10, orthogonal transforms of R n actually act on ODE (1): Lemma 11
S ⊃ O ( n ) ≡ { S ∈ End R n | T SS = id R n } holds for the harmonic oscillatorsystem. roof. Take any S ∈ O ( n ). Since U ( S ( x )) = U ( x ) and K ( S ( p )) = K ( p ) hold forany ( x, p ), and since M − and F become diagonal, the conditions in definition 8 are valid,indicating O ( n ) ⊂ S .We then put R ≡ O ( n ) \{ id R n } =: O ( n ) × , viz., we sum up Ω S to make A R for all S that is a non-identical orthogonal transform of R n . We also restrict our consideration for n = 2 for ease, wherein the discussion would beextended to a larger n . The following proposition shows how Υ ± arises from R . Proposition 12
It holds that A R = Υ , and the decomposition (29) holds with A = A R and A [ ± ] = Υ ± . Proof.
By using the fact that O (2) is bijectively parametrized by S and signature, S × {± } → O (2) , ( θ, σ ) cos θ − σ sin θ sin θ σ cos θ =: s σθ , and by solving an eigenvalue problem s σθ ( p ) = p , we seeΓ s +1 θ = R if θ = 0 { , } otherwise , Γ s − θ = R × { } if θ = 0 { } × R if θ = π (cid:8) ( p , p ) ∈ R | p = − cos θ sin θ p (cid:9) otherwise . Note that Γ s − θ is a line through origin of R with a gradient k = k θ for θ ∈ S \{ , π } , where k θ can take any value in R × , and that Γ s − = R × { } and Γ s − π = { } × R are also lineswith gradients 0 and ∞ , respectively. Thus Γ s − θ = L k , a line through origin of R with agradient k ∈ ( −∞ , ∞ ]. Applying R = O (2) × = { s +1 θ | θ ∈ S × } ∪ { s − θ | θ ∈ S } , we hence21et [ S ∈R (Γ Ds × Γ s )= [ θ ∈ S × (Γ s +1 θ × Γ s +1 θ ) ∪ [ θ ∈ S (Γ s − θ × Γ s − θ )= [ θ ∈ S (Γ s − θ × Γ s − θ )= [ k ∈ ( −∞ , ∞ ] ( L k × L k )= { ( x , kx , p , kp ) | x , p , k ∈ R }∪ { } × R × { } × R = { ( x , x , p , p ) ∈ R | x p − x p = 0 } . Thus A R = S S ∈R (Γ Ds × Γ s ) × R m = { ( x, p ) ∈ R | x p − x p = 0 } × R m = { ω ∈ Ω | γ ( ω ) = 0 } = Υ . Therefore, decomposition (29) holds as Ω = Υ ⊔ Υ +12 ⊔ Υ − .Note that the explicit form of Γ s in the proof directly indicates nontrivial examples toexplain that R should be sufficiently large. For example, if we take R as a one point set,then A R = S S ∈R (Γ Ds × Γ s ) × R m does not separate the phase space into the three spacesas described in (29): if we set R = { s +1 θ } with θ = 0 or R = { s − θ } , then S S ∈R (Γ Ds × Γ s ) =Γ s +1 θ × Γ s +1 θ = { , , , } ⊂ R or S S ∈R (Γ Ds × Γ s ) = Γ s − θ × Γ s − θ =line × line ⊂ R , respectively,clearly induces no separation. So does for e.g. any finite set R = { s ± θ , · · · , s ± θ q } . ACKNOWLEDGMENTS
This work was supported by a Grant-in-Aid for Scientific Research (C) (17K05143 and20K11854) from JSPS and the “Development of innovative drug discovery technologies formiddle-sized molecules” from Japan Agency for Medical Research and development, AMED.We thank Profs. Haruki Nakamura and Akinori Kidera for their continuous encouragement. [1] Allen M. and Tildesley D.
Computer Simulation of Liquids . New York: Oxford; 2002.[2] Collins P., Ezra G.S., and Wiggins S., Phase Space Structure and Dynamics for the Hamilto-nian Isokinetic Thermostat, J. Chem. Phys., 2010, vol. 133, no. 1, pp. 014105.
3] Dobbins S.E., Lesk V.I., and Sternberg M.J., Insights into Protein Flexibility: The Rela-tionship between Normal Modes and Conformational Change Upon Protein-Protein Docking,Proc. Natl. Acad. Sci. U S A, 2008, vol. 105, no. 30, pp. 10390-10395.[4] Ezra G.S., Reversible Measure-Preserving Integrators for Non-Hamiltonian Systems, J. Chem.Phys., 2006, vol. 125, no. 3, pp. 34104.[5] Dettmann C. and Morriss G., Hamiltonian Formulation of the Gaussian Isokinetic Thermostat,Phys. Rev. E, 1996, vol. 54, no. 3, pp. 2495.[6] Fukuda I., Comment on ”Preserving the Boltzmann Ensemble in Replica-Exchange MolecularDynamics” [J. Chem. Phys. 129, 164112 (2008)], J. Chem. Phys., 2010, vol. 132, no. 12, pp.127101.[7] Fukuda I., Coupled Nos´e-Hoover Lattice: A Set of the Nos´e-Hoover Equations with DifferentTemperatures, Phys. Lett. A, 2016, vol. 380, no. pp. 2465-2474.[8] Fukuda I., Symmetric, Explicit Numerical Integrator for Molecular Dynamics Equations ofMotion with a Generalized Friction, Journal of Mathematical Physics, 2019, vol. 60, no. 4,pp. 042903.[9] Fukuda I. and Moritsugu K., Coupled Nos´e-Hoover Equations of Motions without Time Scal-ing, J. Phys. A: Math. Theor., 2017, vol. 50, no. 1, pp. 015002.[10] Fukuda I. and Moritsugu K., Coupled Nos´e-Hoover Equations of Motion to Implement aFluctuating Heat-Bath Temperature, Phys. Rev. E, 2016, vol. 93, no. 3, pp. 033306.[11] Fukuda I. and Moritsugu K., Double Density Dynamics: Realizing a Joint Distribution of aPhysical System and a Parameter System, J. Phys. A: Math. Theor., 2015, vol. 48, no. 45,pp. 455001.[12] Fukuda I. and Nakamura H., Construction of an Extended Invariant for an Arbitrary OrdinaryDifferential Equation with Its Development in a Numerical Integration Algorithm, Phys. Rev.E, 2006, vol. 73, no. 2, pp. 026703.[13] Fukuda I. and Nakamura H., Tsallis Dynamics Using the Nos´e-Hoover Approach, Phys. Rev.E, 2002, vol. 65, no. 2, pp. 026105.[14] H¨uenberger P.H.,
Thermostat Algorithms for Molecular Dynamics Simulations . In: AdvancedComputer Simulation: Approaches for Soft Matter Sciences I. Holm C, Kremer K (editors).Berlin: Springer; 2005. pp. 105-149.[15] Harish M.S. and Patra P.K., Temperature and Its Control in Molecular Dynamics Simulations, Computational Statistical Mechanics . Amsterdam: Elsevier; 1991.[18] Hoover W.G.
Molecular Dynamics . Berlin: Springer-Verlag; 1986.[19] Hoover W.G. and Holian B.L., Kinetic Moments Method for the Canonical Ensemble Distri-bution, Phys. Lett. A, 1996, vol. 211, no. 5, pp. 253-257.[20] Jarzynski C., Nonequilibrium Equality for Free Energy Differences, Phys. Rev. Lett., 1997,vol. 78, no. 14, pp. 2690.[21] Jepps O.G. and Rondoni L., Deterministic Thermostats, Theories of Nonequilibrium Systemsand Parallels with the Ergodic Condition, J. Phys. A: Math. Theor., 2010, vol. 43, no. 13, pp.133001.[22] Krajˇn´ak V., Ezra G.S., and Wiggins S., Roaming at Constant Kinetic Energy: Chesnavich’sModel and the Hamiltonian Isokinetic Thermostat, Regular and Chaotic Dynamics, 2019, vol.24, no. 6, pp. 615-627.[23] Legoll F., Luskin M., and Moeckel R., Non-Ergodicity of the Nos´e-Hoover ThermostattedHarmonic Oscillator, Archive for Rational Mechanics and Analysis, 2006, vol. 184, no. 3, pp.449-463.[24] Leimkuhler B., Noorizadeh E., and Theil F., A Gentle Stochastic Thermostat for MolecularDynamics, J. Stat. Phys., 2009, vol. 135, no. 2, pp. 261-277.[25] Liu Y. and Tuckerman M.E., Generalized Gaussian Moment Thermostatting: A New Contin-uous Dynamical Approach to the Canonical Ensemble, J. Chem. Phys., 2000, vol. 112, no. 4,pp. 1685-1700.[26] Martyna G.J., Klein M.L., and Tuckerman M., Nos´e–Hoover Chains: The Canonical EnsembleVia Continuous Dynamics, J. Chem. Phys., 1992, vol. 97, no. 4, pp. 2635.[27] Moritsugu K., Miyashita O., and Kidera A., Vibrational Energy Transfer in a ProteinMolecule, Phys. Rev. Lett., 2000, vol. 85, no. 18, pp. 3970-3973.[28] Nos´e S., Dynamical Behavior of a Thermostated Isotropic Harmonic Oscillator, Phys. Rev. E,1993, vol. 47, no. 1, pp. 164-177.[29] Nos´e S., A Unified Formulation of the Constant Temperature Molecular Dynamics Methods,J. Chem. Phys., 1984, vol. 81, no. 1, pp. 511-519.
30] Nos´e S.,
Constant Temperature Molecular-Dynamics Methods , Progress Theor. Phys. Suppl.1991, pp. 103.[31] Patra P.K. and Bhattacharya B., Nonergodicity of the Nose-Hoover Chain Thermostat inComputationally Achievable Time, Phys. Rev. E, 2014, vol. 90, no. 4, pp. 043304.[32] Posch H., Hoover W., and Vesely F., Canonical Dynamics of the Nos´e Oscillator: Stability,Order, and Chaos, Phys. Rev. A, 1986, vol. 33, no. 6, pp. 4253-4265.[33] Samoletov A.A., Dettmann C.P., and Chaplain M.A.J., Thermostats for “Slow” Configura-tional Modes, J. Stat. Phys., 2007, vol. 128, no. 6, pp. 1321-1336.[34] Schlick T. Molecular Modeling and Simulation: An Interdisciplinary Guide. New York:Springer Science & Business Media; 2006.[35] Tirion M.M., Large Amplitude Elastic Motions in Proteins from a Single-Parameter, AtomicAnalysis, Phys. Rev. Lett., 1996, vol. 77, no. 9, pp. 1905-1908.[36] Zheng W. and Thirumalai D., Coupling between Normal Modes Drives Protein Conforma-tional Dynamics: Illustrations Using Allosteric Transitions in Myosin Ii, Biophys. J., 2009,vol. 96, no. 6, pp. 2128-2137. igure Captions Fig. 1. Simulation results obtained by a conventional thermostat method (NHC with thechain length 2) for the 2HO using an initial condition x (0) = x (0) = 0 , p (0) = p (0) = 1:(a) trajectory (scatter plot) of ( x , p ) and marginal distributions for (b) x and (c) p , wherecorresponding theoretical distributions and the discrepancies are also shown. The resultsfor x and p are exactly the same as that for x and p (see text).Fig. 2. Simulation results obtained by the NHC for the 2HO using x (0) = 1 , x (0) = 0 ,p (0) = 0 , p (0) = 0 .
01: (a) marginal distributions for (a) x and (b) p , and trajectory of( ζ , ζ ).Fig. 3. Trajectories of γ ( ω ) obtained by the NHC for the 2HO using initial conditionsof (a) that in Fig. 1, (b) that in Fig. 2, and (c) x (0) = 1 , x (0) = 0 , p (0) = 0 , p (0) = − x (0) = x (0) = 0 , p (0) = p (0) = 1: (a) trajectories (scatter plot) of( x , p ) and (b) ( x , p ), and the marginal distributions for (c) x and x and (d) p and p (theoretical distributions and the discrepancies are also shown); (e) trajectories of ( ζ , ζ )and (f) their marginal distributions represented in the principal components y and y ; (g)trajectories of γ ( ω ) and its time average.Fig. 5. Simulation results obtained by the splitting NH for the 2HO, using the differentsetting of the Nos´e mass matrix Q than that in Fig. 4: (a) trajectory of ( ζ , ζ ) and thecontours of the theoretical distribution (rotated ellipsoids); the marginal distributions forprincipal components (b) y and (c) y .Fig. 6. Simulation results obtained by the splitting NH for the 3HO: trajectories anddistributions for (a) x i , (b) p i , and (c) y i ( i = 1 ,,
01: (a) marginal distributions for (a) x and (b) p , and trajectory of( ζ , ζ ).Fig. 3. Trajectories of γ ( ω ) obtained by the NHC for the 2HO using initial conditionsof (a) that in Fig. 1, (b) that in Fig. 2, and (c) x (0) = 1 , x (0) = 0 , p (0) = 0 , p (0) = − x (0) = x (0) = 0 , p (0) = p (0) = 1: (a) trajectories (scatter plot) of( x , p ) and (b) ( x , p ), and the marginal distributions for (c) x and x and (d) p and p (theoretical distributions and the discrepancies are also shown); (e) trajectories of ( ζ , ζ )and (f) their marginal distributions represented in the principal components y and y ; (g)trajectories of γ ( ω ) and its time average.Fig. 5. Simulation results obtained by the splitting NH for the 2HO, using the differentsetting of the Nos´e mass matrix Q than that in Fig. 4: (a) trajectory of ( ζ , ζ ) and thecontours of the theoretical distribution (rotated ellipsoids); the marginal distributions forprincipal components (b) y and (c) y .Fig. 6. Simulation results obtained by the splitting NH for the 3HO: trajectories anddistributions for (a) x i , (b) p i , and (c) y i ( i = 1 ,, ,,