**Numerical Integration Techniques Based on a Geometric View and Application to Molecular Dynamics Simulations**

Ikuo Fukuda<sup>1</sup> and Séverine Queyroy2

<sup>1</sup>*Computational Science Research Program, RIKEN (The Institute of Physical and Chemical Research)* <sup>2</sup>*Laboratoire chimie Provence, chimie théorique, UMR 6264, Aix-Marseille Universités-CNRS* <sup>1</sup>*Japan* <sup>2</sup>*France*

#### **1. Introduction**

42 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy

Yang, J., Hu, W., Tang, J., Xu, M. (2008). Long-Time Scale Molecular Dynamics Study of

Yang, X., Liu, L., Zhang, Q., Zhai, P. (2011). Molecular Dynamics Study of the Structural and

Yang, Z., Yang, X., Xu, Z. (2008). Molecular Dynamics Simulation of the Melting Behavior of Pt-Au Nanoparticles with Core-Shell Structure. *J. Phys. Chem. C* 112: 4937-4947. Yang, Z., Yang, X., Xu, Z., Liu, S. (2009). Structural evolution of Pt-Au nanoalloys during

Yao, C.-h., Song, B., Cao, P.-l. (2004). Structure of Al19 cluster: A full-potential LMTO

Yildirim, E. K., Atis, M., Guvenc, Z. B. (2007). Molecular dynamics simulation of melting behaviour of small gold clusters: AuN (N=12-14). *Phys. Scr.* 75: 111-118. Yukna, J. & Wang, L. (2007). Molecular Dynamics Studies of the Coalescence of Siliver

Zeng, Q., Jiang, X., Yu, A., Lu, G. (2007). Growth mechanisms of silver nanoparticles: a

Zhang, L., Zhang, C.-B., Qi, Y. (2009). Molecular-dynamics investigation of structural transformations of a Cu201cluster in its melting process. *Physica B* 404: 205-209. Zhang, W., Ge, Q., Wang, L. (2003). Structure effects on the energetic, electronic, and magnetic properties of palladium nanoparticles. *J. Chem. Phys.* 118: 5793-5780. Zhang, W., Ran, X., Zhao, H., Wang, L. (2004a). The non-metallicity of molybdenium

Zhang, W., Xiao, L., Hirata, Y., Pawluk, T., Wang, L. (2004b). The simple cubic structure of Ir clusters and the element effect on cluster structures. *Chem. Phys. Lett.* 383: 67-71. Zhang, W., Zhang, F., Zhu, Z. (2006). Molecular dynamics study on the melting phase transition of aluminum clusters with around 55 atoms. *Phys. Rev. B* 74: 033412. Zhang, W., Zhao, H., Wang, L. (2004c). The simple cubic structure of ruthenium clusters. *J.* 

Zhao, S. J., Wang, S. Q., Yang, Z. Q., Ye, H. Q. (2001). Coalescence of three silver nanoclusters: a molecular dynamics study. *J. Phys.: Condens. Matter* 13: 8061-8069.

molecular-dynamics study. *Phys. Rev. B* 70: 195431.

molecular dynamics study. *Nanotechnology* 18: 035708.

Clusters. *J. Phys. Chem. C* 111: 13337-13347.

clusters. *J. Chem. Phys.* 121: 7717-7724.

*Phys. Chem. B* 108: 2140-2147.

2074-2078.

489-492.

*Phys.* 11: 6249-6255.

Diffusion Dynamics of Small Cu Clusters on Cu(111) Surface. *J. Phys. Chem. C* 112:

Mechanical Properties of Skutterudite CoSb3: Surface Effect. *J. Electronic Mat.* 40:

heating process: comparison of random and core-shell orderings. *Phys. Chem. Chem.* 

In this chapter we address numerical integration techniques of ordinary differential equation (ODE), especially that for molecular dynamics (MD) simulation. Since most of the fundamental equations of motion in MD are represented by nonlinear ODEs with many degrees of freedom, numerical integration becomes essential to solve the equations for analyzing the properties of a target physical system. To enhance the molecular simulation performance, we demonstrate two techniques for numerically integrating the ODE. The first object we present is an invariant function, viz., a conserved quantity along a solution, of a given ODE. The second one is a numerical integrator itself, which numerically solves the ODE by capturing certain geometric properties of the ODE.

In our proposed procedure (Fukuda & Nakamura, 2006), for an ODE defined on an *N*-dimensional phase space Ω we construct an extended phase space Ω of *N* + 1 dimension, by introducing an additional degree of freedom. Then, on Ω we constitute a new ODE, which has an invariant but retains every solution of the original ODE, and we construct efficient integrators for the extended ODE. Advantageous features of our proposal are the simplicity and the applicability to a wide class of ODEs beyond the Hamiltonian equations. In fact, in MD methods, non-Hamiltonian equations are often used (Hoover, 1991); they have been developed (Hoover & Holian, 1996), e.g., to provide more robustness than conventional one or a rapid convergence to a targeted statistical thermodynamic ensemble, or to define a new ensemble itself. Considering such a development, new equations must be designed successively in future studies, and the simplicity and the applicability for the current techniques will be useful also in such a circumstance.

Specifically, by the first technique, the invariant can be simply constructed for any (smooth) ODE, including non-Hamiltonian equation. It can thus be easily used to examine the accuracy of numerical integration of the ODE by monitoring the invariant value, as done in Hamiltonian system by using the Hamiltonian value. Namely, the integration accuracy check can be done in a system that does not have or may not have an invariant.

where a point in Ω� is expressed as *ω*� = (*ω*, *v*) [**R**<sup>×</sup> denotes non-zero numbers; *Di* the partial differentiation with respect to *ω<sup>i</sup>* (all quantities are supposed to be sufficiently smooth)]. Here,

<sup>45</sup> Numerical Integration Techniques Based on a

defines an invariant (dimensionless); i.e., for an arbitrary solution *φ*� of Eq. (3), the value of

The way of defining the extended equation is not unique and Eq. (3b) is constructed with the aim that we will make a volume-preserving integrator, as demonstrated in the following subsection. In fact, alternative schemes are discussed in the literature (Fukuda & Nakamura, 2006), and by applying the scheme to several MD equations the conserved quantity individually defined in each equation can be reproduced in a uniform, generalized manner. Concerning the technical issue, we assume that, for simplicity, the Liouville equation,

holds for a certain density function *ρ*, and put *B*≡ − ln *ρ* as the choice of *B*. Thus, it follows

The above conception is illustrated in figure 1: A new ODE (2) in an extended phase space Ω� is defined such that (i) a projection of any solution *φ*� onto the original phase space Ω is a solution of the original ODE (1), and (ii) a function (invariant) *L* exists, whose value is constant

, and the similar holds for *ω<sup>t</sup>* ≡ *φ*(*t*)). The above (i) is a trivial matter, since Eq. (3a) in the extended ODE is the original ODE itself, which does not couple to the new variable *v*. By the suitable generation of the field (viz., *X*e) in the extended dimension, the above (ii) is attained,

> ( ) *X* Z

( ) *X* Z*<sup>t</sup>* c c

*<sup>t</sup> <sup>X</sup>* ( )

*<sup>t</sup> X* ( ) Z

c c c

Z*t*

Fig. 1. A schematic figure to illustrate the basic concept for an invariant and extended phase

As was stated, a projection of any solution with initial value (*ω*0, *v*0) for the extended ODE (3) onto Ω is also a solution with initial value *ω*<sup>0</sup> for the original ODE (3a). Conversely, any solution of the original ODE can be lifted to the solution of the extended ODE. However, this correspondence is not one-to-one, i.e., many lifts are possible according to the choice of

Z

Z

Z

Zc *L*(*ω*, *v*) = *B*(*ω*) + ln |*v*| (4)

div(*ρX*) = 0, (5)

(*t*) is a point of the solution at time *t* with initial value

: Level set *<sup>c</sup>*

( ) for all *Lc t* Z*t* c

L

: : Original phase space

*B* is any phase-space function. Then we can show that

Geometric View and Application to Molecular Dynamics Simulations

) = −div*X*(*ω*)*v*.

*<sup>t</sup>* ≡ *φ*�

(*t*)) is constant for all time *t*.

from Eq. (3b) that *X*e(*ω*�

along any solution (in the figure, *ω*�

as straightforwardly confirmed.

: Added space

*L*(*φ*�

*ω*�

space

The second technique is to constitute the integrator that is widely applicable to many ODEs, including those for non-Hamiltonian systems. We require higher accuracy to avoid the accumulation of the numerical errors, as the simulation time increases. In addition, we should seriously address the issue of the computational cost: with the advancement of computational architecture, larger and larger systems can be treated and the computational cost grows consequently. In fact, the number of the interaction evaluations, which characterize the total cost in MD calculations, is considerably grown with increasing the number of degrees of freedom of the system, *n*; e.g. the number of the evaluation is of order of *n*<sup>2</sup> in a typical classical system. We provide a route to easily construct efficient integrators for many kinds of ODEs. Specifically, on the extended space we present integrators that are explicit, symmetric, and phase space volume-preserving.

Geometric features of integrator, including the volume preserving property (Zai-jiu, 1994; Kang & Zai-jiu, 1995, Quispel, 1995; Okunbor, 1995; Quispel & Dyt, 1998), give rise to stable simulation, as if the symplectic integrator does for Hamiltonian system (Ruth, 1983; Yoshida, 1990; McLachlan & Atela, 1992; Sanz-Serna, 1992). In fact, volume-preserving property is a generalization of the symplectic property; this corresponds to the fact that a divergence-free system, which is realized in the extended system by our protocol, is a generalization of a Hamiltonian system. For relevant geometric concepts in the analysis of molecular dynamics equations, see e.g., the work of Ezra (Ezra, 2006) and the references therein.

In section 2, we demonstrate the details of the integration method and explain the geometric view lying in our technique. To effectively present such a view in MD study, the implementation of algorithms is explicitly described in reference to the Nosé-Hoover (NH) equation (Nosé, 1984; Hoover, 1985), which is one of the representative ODEs in MD studies. In section 3 we discuss the applicability of our method and clarify the issues that should be considered for a further development of the method.

#### **2. Geometric concept in the integration techniques**

#### **2.1 Invariant on extended space: Fiber bundle structure**

Beginning with a review of our work (Fukuda & Nakamura, 2006) concerning the construction of the invariant, we give a new interpretation of this matter to clarify the geometric view under consideration. The main idea is to suitably generate a vector field in the extended dimension so that the extended ODE has an invariant.

For any ODE in phase space Ω (domain of **R***N*),

$$
\dot{\omega} = X(\omega),
\tag{1}
$$

$$\text{viiz.}\,d\omega\_i/dt = \text{X}\_i(\omega\_1, \dots, \omega\_N) \text{ for } i = 1, \dots, N\_\prime \text{ consider an extended ODE}$$

$$
\dot{\omega}' = X'(\omega'),
\tag{2}
$$

which is defined, on an extended phase space Ω� ≡ Ω × **R**×, by

$$
\dot{\omega} = X(\omega),
\tag{3a}
$$

$$\dot{\boldsymbol{v}} = \boldsymbol{X\_{e}}(\boldsymbol{\omega}') \equiv -\sum\_{i=1}^{N} \boldsymbol{X\_{i}}(\boldsymbol{\omega}) \, \boldsymbol{D\_{i}} \boldsymbol{B}(\boldsymbol{\omega}) \, \boldsymbol{v\_{\star}} \tag{3b}$$

2 Will-be-set-by-IN-TECH

Hamiltonian system by using the Hamiltonian value. Namely, the integration accuracy check

The second technique is to constitute the integrator that is widely applicable to many ODEs, including those for non-Hamiltonian systems. We require higher accuracy to avoid the accumulation of the numerical errors, as the simulation time increases. In addition, we should seriously address the issue of the computational cost: with the advancement of computational architecture, larger and larger systems can be treated and the computational cost grows consequently. In fact, the number of the interaction evaluations, which characterize the total cost in MD calculations, is considerably grown with increasing the number of degrees of freedom of the system, *n*; e.g. the number of the evaluation is of order of *n*<sup>2</sup> in a typical classical system. We provide a route to easily construct efficient integrators for many kinds of ODEs. Specifically, on the extended space we present integrators that are explicit, symmetric,

Geometric features of integrator, including the volume preserving property (Zai-jiu, 1994; Kang & Zai-jiu, 1995, Quispel, 1995; Okunbor, 1995; Quispel & Dyt, 1998), give rise to stable simulation, as if the symplectic integrator does for Hamiltonian system (Ruth, 1983; Yoshida, 1990; McLachlan & Atela, 1992; Sanz-Serna, 1992). In fact, volume-preserving property is a generalization of the symplectic property; this corresponds to the fact that a divergence-free system, which is realized in the extended system by our protocol, is a generalization of a Hamiltonian system. For relevant geometric concepts in the analysis of molecular dynamics

In section 2, we demonstrate the details of the integration method and explain the geometric view lying in our technique. To effectively present such a view in MD study, the implementation of algorithms is explicitly described in reference to the Nosé-Hoover (NH) equation (Nosé, 1984; Hoover, 1985), which is one of the representative ODEs in MD studies. In section 3 we discuss the applicability of our method and clarify the issues that should be

Beginning with a review of our work (Fukuda & Nakamura, 2006) concerning the construction of the invariant, we give a new interpretation of this matter to clarify the geometric view under consideration. The main idea is to suitably generate a vector field in the extended dimension

*ω*˙ � = *X*�

) ≡ −

(*ω*�

*N* ∑ *i*=1

*ω*˙ = *X*(*ω*), (1)

*ω*˙ = *X*(*ω*), (3a)

), (2)

*Xi*(*ω*) *DiB*(*ω*) *v*, (3b)

can be done in a system that does not have or may not have an invariant.

equations, see e.g., the work of Ezra (Ezra, 2006) and the references therein.

viz., *dωi*/*dt* = *Xi*(*ω*1, ...., *ωN*) for *i* = 1, ..., *N*, consider an extended ODE

which is defined, on an extended phase space Ω� ≡ Ω × **R**×, by

*v*˙ = *X*e(*ω*�

considered for a further development of the method.

so that the extended ODE has an invariant. For any ODE in phase space Ω (domain of **R***N*),

**2. Geometric concept in the integration techniques 2.1 Invariant on extended space: Fiber bundle structure**

and phase space volume-preserving.

where a point in Ω� is expressed as *ω*� = (*ω*, *v*) [**R**<sup>×</sup> denotes non-zero numbers; *Di* the partial differentiation with respect to *ω<sup>i</sup>* (all quantities are supposed to be sufficiently smooth)]. Here, *B* is any phase-space function. Then we can show that

$$L(\omega, v) = B(\omega) + \ln|v| \tag{4}$$

defines an invariant (dimensionless); i.e., for an arbitrary solution *φ*� of Eq. (3), the value of *L*(*φ*� (*t*)) is constant for all time *t*.

The way of defining the extended equation is not unique and Eq. (3b) is constructed with the aim that we will make a volume-preserving integrator, as demonstrated in the following subsection. In fact, alternative schemes are discussed in the literature (Fukuda & Nakamura, 2006), and by applying the scheme to several MD equations the conserved quantity individually defined in each equation can be reproduced in a uniform, generalized manner. Concerning the technical issue, we assume that, for simplicity, the Liouville equation,

$$\text{div}(\rho X) = 0,\tag{5}$$

holds for a certain density function *ρ*, and put *B*≡ − ln *ρ* as the choice of *B*. Thus, it follows from Eq. (3b) that *X*e(*ω*� ) = −div*X*(*ω*)*v*.

The above conception is illustrated in figure 1: A new ODE (2) in an extended phase space Ω� is defined such that (i) a projection of any solution *φ*� onto the original phase space Ω is a solution of the original ODE (1), and (ii) a function (invariant) *L* exists, whose value is constant along any solution (in the figure, *ω*� *<sup>t</sup>* ≡ *φ*� (*t*) is a point of the solution at time *t* with initial value *ω*� , and the similar holds for *ω<sup>t</sup>* ≡ *φ*(*t*)). The above (i) is a trivial matter, since Eq. (3a) in the extended ODE is the original ODE itself, which does not couple to the new variable *v*. By the suitable generation of the field (viz., *X*e) in the extended dimension, the above (ii) is attained, as straightforwardly confirmed.

Fig. 1. A schematic figure to illustrate the basic concept for an invariant and extended phase space

As was stated, a projection of any solution with initial value (*ω*0, *v*0) for the extended ODE (3) onto Ω is also a solution with initial value *ω*<sup>0</sup> for the original ODE (3a). Conversely, any solution of the original ODE can be lifted to the solution of the extended ODE. However, this correspondence is not one-to-one, i.e., many lifts are possible according to the choice of *v*0. This correspondence can be naturally understood via a geometric concept, fiber bundle (e.g., Husemoller, 1966), which plays an important role in physics (see Choquet-Bruhat et al., 1982), including gauge theory, Berry's phase, and the quantum Hall effect, as well as in mathematics. In this context, our target is a product (trivial) vector bundle over Ω with fiber **R**, i.e., (*E* ≡ Ω × **R**, *π*, Ω), where we have considered that the total space is Ω × **R** rather than <sup>Ω</sup> <sup>×</sup> **<sup>R</sup>**×. A level set (i.e., a generalization of energy surface) is represented by <sup>L</sup>*<sup>c</sup>* <sup>=</sup> <sup>L</sup><sup>+</sup> *<sup>c</sup>* ∪ L<sup>−</sup> *c* , where

$$\mathcal{L}\_{\mathfrak{c}}^{\pm} \equiv \left\{ \omega' \in \Omega' \mid L(\omega') = c, v \gtrless 0 \right\}. \tag{6}$$

:c :c

method, estimate the succeeding value by using the several values previously sought with

To explain the procedures in a specific manner, we shall use the NH equation (Nosé, 1984; Hoover, 1985; Nosé, 1991) as an example of ODE. This is because the NH equation is frequently used in MD simulation to control the temperature, and furthermore it gives a foundation of various realizations of the Boltzmann-Gibbs dynamics (see e.g., Hoover, 1991;

*<sup>p</sup>* · **<sup>M</sup>**<sup>−</sup>1, −∇*U*(*x*) <sup>−</sup> (*ζ*/*Q*)*p*, 2*K*(*p*) <sup>−</sup> *nk*B*<sup>T</sup>*

where *ω* ≡ (*x*, *p*, *ζ*) and **M** ≡diag(*m*1, ..., *mn*); *x* ≡ (*x*1,..., *xn*), *p* ≡ (*p*1,..., *pn*), *U*(*x*), and

respectively, of a physical system; *k*<sup>B</sup> is Boltzmann's constant, *ζ* a real variable to control the temperature of the physical system to targeted temperature *T*, and *Q* a positive parameter

*ω*˙ = *X*NH(*ω*) (11a)

*<sup>i</sup>* /2*mi* represent the coordinates, momenta, potential energy, and kinetic energy,

<sup>−</sup>[*U*(*x*) + *<sup>K</sup>*(*p*) + *<sup>ζ</sup>*2/2*Q*]/*k*B*<sup>T</sup>*

/*k*B*T* + ln |*v*|. (14)

. (15)

*v*˙ = − div *X*(*ω*)*v* = (*nζ*/*Q*)*v*, (13)

, (11b)

(12)

Fig. 2. A schematic picture to show that one-step map integrator is a mapping of the

Nosé, 1991; Fukuda & Nakamura, 2004). The NH equation can be represented by

<sup>1</sup> *T* { <*<sup>h</sup>*

<sup>47</sup> Numerical Integration Techniques Based on a

u

u

u

Geometric View and Application to Molecular Dynamics Simulations

Zc

(extended) phase space onto itself.

certain interpolation techniques.

*<sup>K</sup>*(*p*) <sup>≡</sup> <sup>∑</sup>*<sup>n</sup>*

*<sup>i</sup>*=<sup>1</sup> *<sup>p</sup>*<sup>2</sup>

and invariant (4) becomes

associated with *ζ*. The function

≡ 

*<sup>ρ</sup>* (*ω*) <sup>≡</sup> exp

) =

(i) For a given ODE (1), consider a solvable decomposition

*L*(*ω*�

Now, the five procedures are described as follows:

is a density for the Liouville equation. Thus, extended equation (3b) becomes

*U*(*x*) + *K*(*p*) + *ζ*2/2*Q*

*<sup>X</sup>* <sup>=</sup> *<sup>X</sup>*[1] <sup>+</sup> *<sup>X</sup>*[2] <sup>+</sup> ··· <sup>+</sup> *<sup>X</sup>*[*L*]

u

In each piece, i.e., one component of a level set, the invariant function takes a constant value, *c*. It should be noted that each piece becomes a (image of) global cross section from a base space Ω, where the cross section is a map defined by

$$s\_c^{\pm}: \Omega \to \mathbb{E}, \ \omega \mapsto (\omega, \pm \exp(-\mathcal{B}(\omega) + \mathfrak{c})).\tag{7}$$

This indicates that L<sup>±</sup> *<sup>c</sup>* is equivalent, in a topological sense, to the original space Ω. This situation is different in the Hamiltonian system in which a level set {(*x*, *<sup>p</sup>*) <sup>∈</sup> **<sup>R</sup>**2*<sup>n</sup>* <sup>|</sup> *<sup>H</sup>*(*x*, *<sup>p</sup>*) = *<sup>c</sup>*} often becomes compact and is not necessarily equivalent to **<sup>R</sup>**2*n*−<sup>1</sup> nor **<sup>R</sup>**2*n*. The choice of *v*0, under an arbitrary fixed *ω*0, corresponds to the choice of the level set, i.e., the choice of the parameter *c* and the signature, which also characterizes the section map. The time development in the extended space is always on the (image of the) cross section.

#### **2.2 Integrator on extended space: Volume-preserving discrete dynamical-system structure**

The second technique is on the numerical integrator. We explain a general idea for constructing the current integrator through the following five procedures:

(i) Considering a solvable decomposition for the ODE originally given,

(ii) Making a solvable decomposition for the extended ODE,

(iii) Turning attention from the solutions of the (decomposed) ODE to phase-space maps in the extended space,

(iv) Combining the individual maps to get a first order integrator,

(v) Combining the maps furthermore to get a higher order integrator.

Our view is clear: the integrator Ψ*<sup>h</sup>* is a (one-step) invertible map (see figure 2) parametrized by a time step *h*. In metaphorical terms, we have the most fundamental geometric view as

$$\text{Integration} \equiv \{ \overline{T}\_{\mathfrak{m}} \equiv \Psi\_{\mathfrak{h}} \circ \stackrel{\mathfrak{m}}{\cdots} \circ \Psi\_{\mathfrak{h}} \text{ (iteration of the map)} : \Omega' \to \Omega' \}\_{\mathfrak{m} \text{ integers}} \tag{8}$$

$$\equiv \text{Discrete dynamical system (generated by map)}.\tag{9}$$

Namely, the iterations of the map <sup>Ψ</sup>*<sup>h</sup>* constitute the integration [*T*<sup>0</sup> <sup>≡</sup>identity, *<sup>T</sup>*−*<sup>m</sup>* <sup>≡</sup> (*Tm*)−<sup>1</sup> can be defined]. This iteration corresponds to an approximation to

$$\begin{aligned} \text{Exact flow} & \equiv \{ T\_l : \Omega' \to \Omega' \} \_{l \text{: real numbers}} \\ & \equiv \text{Continuous dynamical system (generated by ODE)}. \end{aligned} \tag{10}$$

Of course, this view is applicable for a one-step map integrator and not necessarily valid for an arbitrary integrator; e.g., linear multistep methods, including the predictor-corrector 4 Will-be-set-by-IN-TECH

*v*0. This correspondence can be naturally understood via a geometric concept, fiber bundle (e.g., Husemoller, 1966), which plays an important role in physics (see Choquet-Bruhat et al., 1982), including gauge theory, Berry's phase, and the quantum Hall effect, as well as in mathematics. In this context, our target is a product (trivial) vector bundle over Ω with fiber **R**, i.e., (*E* ≡ Ω × **R**, *π*, Ω), where we have considered that the total space is Ω × **R** rather than <sup>Ω</sup> <sup>×</sup> **<sup>R</sup>**×. A level set (i.e., a generalization of energy surface) is represented by <sup>L</sup>*<sup>c</sup>* <sup>=</sup> <sup>L</sup><sup>+</sup>

*ω*� ∈ Ω� | *L*(*ω*�

In each piece, i.e., one component of a level set, the invariant function takes a constant value, *c*. It should be noted that each piece becomes a (image of) global cross section from a base

situation is different in the Hamiltonian system in which a level set {(*x*, *<sup>p</sup>*) <sup>∈</sup> **<sup>R</sup>**2*<sup>n</sup>* <sup>|</sup> *<sup>H</sup>*(*x*, *<sup>p</sup>*) = *<sup>c</sup>*} often becomes compact and is not necessarily equivalent to **<sup>R</sup>**2*n*−<sup>1</sup> nor **<sup>R</sup>**2*n*. The choice of *v*0, under an arbitrary fixed *ω*0, corresponds to the choice of the level set, i.e., the choice of the parameter *c* and the signature, which also characterizes the section map. The time

**2.2 Integrator on extended space: Volume-preserving discrete dynamical-system structure** The second technique is on the numerical integrator. We explain a general idea for

(iii) Turning attention from the solutions of the (decomposed) ODE to phase-space maps in

Our view is clear: the integrator Ψ*<sup>h</sup>* is a (one-step) invertible map (see figure 2) parametrized by a time step *h*. In metaphorical terms, we have the most fundamental geometric view as

Namely, the iterations of the map <sup>Ψ</sup>*<sup>h</sup>* constitute the integration [*T*<sup>0</sup> <sup>≡</sup>identity, *<sup>T</sup>*−*<sup>m</sup>* <sup>≡</sup> (*Tm*)−<sup>1</sup>

Of course, this view is applicable for a one-step map integrator and not necessarily valid for an arbitrary integrator; e.g., linear multistep methods, including the predictor-corrector

}*t*: real numbers

···◦ Ψ*<sup>h</sup>* (iteration of the map) : Ω� → Ω�

≡ Discrete dynamical system (generated by map). (9)

≡ Continuous dynamical system (generated by ODE). (10)

development in the extended space is always on the (image of the) cross section.

constructing the current integrator through the following five procedures: (i) Considering a solvable decomposition for the ODE originally given,

(ii) Making a solvable decomposition for the extended ODE,

Integration ≡ {*Tm* <sup>≡</sup> <sup>Ψ</sup>*<sup>h</sup>* ◦ *<sup>m</sup>*

(iv) Combining the individual maps to get a first order integrator, (v) Combining the maps furthermore to get a higher order integrator.

can be defined]. This iteration corresponds to an approximation to

Exact flow ≡ {*Tt* : Ω� → Ω�

) = *c*, *v* ≷ 0

*<sup>c</sup>* : Ω → *E*, *ω* �→ (*ω*, ± exp(−*B*(*ω*) + *c*)). (7)

*<sup>c</sup>* is equivalent, in a topological sense, to the original space Ω. This

L<sup>±</sup> *<sup>c</sup>* <sup>≡</sup>

space Ω, where the cross section is a map defined by

*s* ±

where

This indicates that L<sup>±</sup>

the extended space,

*<sup>c</sup>* ∪ L<sup>−</sup> *c* ,

. (6)

}*m*: integers (8)

Fig. 2. A schematic picture to show that one-step map integrator is a mapping of the (extended) phase space onto itself.

method, estimate the succeeding value by using the several values previously sought with certain interpolation techniques.

To explain the procedures in a specific manner, we shall use the NH equation (Nosé, 1984; Hoover, 1985; Nosé, 1991) as an example of ODE. This is because the NH equation is frequently used in MD simulation to control the temperature, and furthermore it gives a foundation of various realizations of the Boltzmann-Gibbs dynamics (see e.g., Hoover, 1991; Nosé, 1991; Fukuda & Nakamura, 2004). The NH equation can be represented by

$$
\dot{\omega} = X\_{\rm NH}(\omega) \tag{11a}
$$

$$\mathbf{x} \equiv \left( p \cdot \mathbf{M}^{-1} , -\nabla \mathcal{U}(\mathbf{x}) - (\mathbb{\zeta}/\mathbb{Q}) p \,\prime \mathcal{R} \mathbf{K}(p) - n k\_{\mathbb{B}} T \right) , \tag{11b}$$

where *ω* ≡ (*x*, *p*, *ζ*) and **M** ≡diag(*m*1, ..., *mn*); *x* ≡ (*x*1,..., *xn*), *p* ≡ (*p*1,..., *pn*), *U*(*x*), and *<sup>K</sup>*(*p*) <sup>≡</sup> <sup>∑</sup>*<sup>n</sup> <sup>i</sup>*=<sup>1</sup> *<sup>p</sup>*<sup>2</sup> *<sup>i</sup>* /2*mi* represent the coordinates, momenta, potential energy, and kinetic energy, respectively, of a physical system; *k*<sup>B</sup> is Boltzmann's constant, *ζ* a real variable to control the temperature of the physical system to targeted temperature *T*, and *Q* a positive parameter associated with *ζ*. The function

$$\rho\left(\omega\right) \equiv \exp\left(-\left[\mathcal{U}(\mathbf{x}) + \mathcal{K}(p) + \zeta^2 / 2Q\right] / k\_\mathcal{B} T\right) \tag{12}$$

is a density for the Liouville equation. Thus, extended equation (3b) becomes

$$
\dot{v} = -\operatorname{div} X(\omega) v = (n\zeta/Q)v,\tag{13}
$$

and invariant (4) becomes

$$L(\omega') = \left[\mathcal{U}(x) + K(p) + \zeta^2 / 2Q\right] / k\_\mathcal{B}T + \ln|v|.\tag{14}$$

Now, the five procedures are described as follows:

(i) For a given ODE (1), consider a solvable decomposition

$$X = X^{[1]} + X^{[2]} + \cdots + X^{[L]}.\tag{15}$$

That is, decompose *X* such that each ODE, *ω*˙ = *X*[*j*] (*ω*), can be solved explicitly for all time; viz., decompose the field *X* into "easier" field components, prior to the procedure (see (iv)) in which the individual maps generated by the individual field components are combined.

For the NH equation, a solvable decomposition is given (but not uniquely) by *X*NH = ∑4 *<sup>j</sup>*=<sup>1</sup> *<sup>X</sup>*[*j*] defined as:

$$X^{[1]}(\omega) \equiv \left(p \cdot \mathbf{M}^{-1}, 0, 0\right),\tag{16a}$$

where *ω*� ≡ (*x*, *p*, *ζ*, *v*); and a solution with an initial value *ω*�

Geometric View and Application to Molecular Dynamics Simulations

(*t*) =

(*t*) =

Consider any decomposed component *j*. Even when we rewrite Φ[*j*]

(extended) phase-space point to be impartial. Namely Φ[*j*]

For the NH equation, it follows from Eq. (25) that e.g., Φ[1]

*<sup>p</sup>* <sup>→</sup> <sup>e</sup>−(*ζ*/*Q*)*<sup>t</sup>*

<sup>Φ</sup>*<sup>t</sup>* <sup>≡</sup> <sup>Φ</sup>[1]

Φ*t*(*ω*�

with respect to initial value by fixing time, we get a phase-space map Φ[*j*]

Φ[*j*] *<sup>t</sup>* (*ω*�

Volume-preserving property is ensured by Eq. (20); i.e., for any area *A* ⊂ Ω�

Vol(Φ[*j*]

*<sup>x</sup>* <sup>→</sup> *tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>* by <sup>Φ</sup>[1]

*<sup>p</sup>* → −*t*∇*U*(*x*) + *<sup>p</sup>* by <sup>Φ</sup>[2]

*<sup>t</sup>* ◦ <sup>Φ</sup>[2]

)··· for point *<sup>ω</sup>*�

) − *Tt*(*ω*�

*<sup>ζ</sup>* <sup>→</sup> (2*K*(*p*) <sup>−</sup> *nk*B*T*)*<sup>t</sup>* <sup>+</sup> *<sup>ζ</sup>* by <sup>Φ</sup>[3]

*<sup>x</sup>*0, e<sup>−</sup> *<sup>ζ</sup>*0*<sup>t</sup>*

*tp*<sup>0</sup> · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>*0, *<sup>p</sup>*0, *<sup>ζ</sup>*0, *<sup>v</sup>*<sup>0</sup>

<sup>49</sup> Numerical Integration Techniques Based on a

*<sup>Q</sup> <sup>p</sup>*0, *<sup>ζ</sup>*0, e *<sup>n</sup>ζ*0*<sup>t</sup>*

) <sup>≡</sup> <sup>Φ</sup>[*j*]*ω*�

where we have dropped the suffix "0" in the initial value because we consider all the

to the point reached by the solution with initial value *ω*� for ODE (23) after a period *t*.

*<sup>t</sup>* as

*<sup>p</sup>*, *<sup>v</sup>* <sup>→</sup> <sup>e</sup>(*nζ*/*Q*)*<sup>t</sup>*

*<sup>t</sup>* ◦···◦ <sup>Φ</sup>[*L*]

Through just this procedure we get a first-order integrator for the extended ODE (19). Here, "first-order" means that the map approximates the exact flow of the ODE locally with an order

) = *O*(*t*

.

*<sup>Q</sup> v*<sup>0</sup> 

<sup>0</sup> (*t*) still denotes the solution as long as we focus on the variation with

; here, the variable that changes is only *x*, and in this way we can express

*v* by Φ[4]

*<sup>t</sup>* : Ω� → Ω�

<sup>2</sup>) as *<sup>t</sup>* <sup>→</sup> 0. In general, map <sup>Φ</sup>*<sup>t</sup>* is said to be order *<sup>p</sup>* if

(*t*) = (*x*0, −*t*∇*U*(*x*0) + *p*0, *ζ*0, *v*0), (25b)

(*t*) = (*x*0, *p*0,(2*K*(*p*0) − *nk*B*T*)*t* + *ζ*0, *v*0), (25c)

Φ[1]

Φ[2]

Φ[3]

Φ[4]

respect to time *t* by fixing initial value *ω*�

and the volume of the mapped area are equal:

(iv) Combine the individual maps as

) − *Tt*(*ω*�

) = *O*(*t*

the essential changes provided by each map Φ[*j*]

(iii) Turn attention from solutions to phase-space maps.

is given, respectively, by

the initial value, Φ[*j*]*ω*�

*tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>*, *<sup>p</sup>*, *<sup>ζ</sup>*, *<sup>v</sup>*

giving Φ[1] *t* Φ[2] *t* ··· <sup>Φ</sup>[*L*] *<sup>t</sup>* (*ω*�

1: Φ*t*(*ω*�

of *t*

<sup>0</sup> ≡ (*x*0, *p*0, *ζ*0, *v*0) for ODE (23)

, (25a)

. (25d)

*<sup>t</sup>* : Ω� → Ω� as,

*<sup>t</sup>* maps every point *ω*� in Ω�

*<sup>t</sup>* maps *ω*� = (*x*, *p*, *ζ*, *v*) to

*<sup>t</sup>* , (28a)

*<sup>t</sup>* , (28b)

*<sup>p</sup>*+1) as *<sup>t</sup>* <sup>→</sup> 0. (30)

*<sup>t</sup>* , (28c)

*<sup>t</sup>* . (28d)

, (29)

<sup>0</sup> (*t*) to express

, its volume

(*t*) as Φ[*j*]*ω*�

(*t*), (26)

<sup>0</sup>. On the other hand, if we focus on the variation

*<sup>t</sup>* (*A*)) = Vol(*A*). (27)

$$X^{[2]}(\omega) \equiv \left(0, -\nabla l I(\mathbf{x}), 0\right),\tag{16b}$$

$$X^{[3]}(\omega) \equiv (0, 0, 2\mathcal{K}(p) - nk\_{\mathbb{B}}T) \, , \tag{16c}$$

$$X^{[4]}(\omega) \equiv \left(0, -(\zeta/\mathbb{Q})p, 0\right). \tag{16d}$$

This is indeed the solvable decomposition as is easily seen; e.g., *ω*˙ = *X*[1] (*ω*) is

$$
\dot{\mathbf{x}} = \boldsymbol{p} \cdot \mathbf{M}^{-1} \; , \; \dot{\boldsymbol{p}} = \mathbf{0} \; , \; \dot{\boldsymbol{\zeta}} = \mathbf{0} \; , \tag{17}
$$

so we have its solution *<sup>φ</sup>*[1] with an initial value *<sup>ω</sup>*<sup>0</sup> = (*x*0, *<sup>p</sup>*0, *<sup>ζ</sup>*0) in an explicit form for all time *t*, as

$$\boldsymbol{\phi}^{[1]}(t) = \left(t\boldsymbol{p}\_0 \cdot \mathbf{M}^{-1} + \mathbf{x}\_{0\prime} \ p\_{0\prime} \boldsymbol{\zeta}\_0\right). \tag{18}$$

(ii) Consider a solvable decomposition for, in turn, the extended ODE (2), which is now represented by

$$
\dot{\omega} = X(\omega),
\tag{19a}
$$

$$
\psi = -\text{div} X(\omega) v. \tag{19b}
$$

To obtain an extended phase-space volume-preserving integrator we should devise a decomposition *X*� = ∑ *X*�[*j*] that yields, for each *j*, the divergence-free property,

$$\text{div}X^{\prime[j]} = 0.\tag{20}$$

Fortunately, this can be done automatically; we have a desired decomposition:

$$X' = X^{\prime[1]} + X^{\prime[2]} + \cdots + X^{\prime[L]},\tag{21}$$

$$X^{\prime[j]}(\omega^{\prime}) = (X^{[j]}(\omega), -\text{div} X^{[j]}(\omega)v) \tag{22}$$

Although, in a strict sense, regarding a solution Φ[*j*] of

$$
\dot{\omega}' = X^{\prime[\dot{j}]}(\omega'),
\tag{23}
$$

its *v*-component is generally given via a form that is integrated with respect to time, this form can be evaluated explicitly in many cases as already discussed (Fukuda & Nakamura, 2006).

In fact, for the NH equation, according to Eqs. (16) and (22) we have

$$X^{\prime [1]}(\omega^{\prime}) = \left(p \cdot \mathbf{M}^{-1}, 0, 0, 0\right),\tag{24a}$$

$$X^{\prime[2]}(\omega^{\prime}) = (0, -\nabla \mathcal{U}(\mathbf{x}), 0, 0), \tag{24b}$$

$$X^{\prime [3]}(\omega^{\prime}) = \begin{pmatrix} 0, 0, 2K(p) - nk\_{\mathbb{B}}T, 0 \end{pmatrix},\tag{24c}$$

$$X^{\prime [4]}(\omega^{\prime}) = \left(0, -\frac{\zeta}{\mathcal{Q}}p, 0, \frac{n\zeta}{\mathcal{Q}}v\right),\tag{24d}$$

6 Will-be-set-by-IN-TECH

viz., decompose the field *X* into "easier" field components, prior to the procedure (see (iv)) in which the individual maps generated by the individual field components are combined.

For the NH equation, a solvable decomposition is given (but not uniquely) by *X*NH =

*<sup>p</sup>* · **<sup>M</sup>**<sup>−</sup>1, 0, 0

(*ω*) ≡ (0, −∇*U*(*x*), 0), (16b)

(*ω*) ≡ (0, 0, 2*K*(*p*) − *nk*B*T*), (16c)

(*ω*) ≡ (0, −(*ζ*/*Q*)*p*, 0). (16d)

*ω*˙ = *X*(*ω*), (19a) *v*˙ = −div*X*(*ω*)*v*. (19b)

div*X*�[*j*] = 0. (20)

(*ω*), can be solved explicitly for all time;

, (16a)

(*ω*) is

. (18)

, (21)

, (24a)

, (24d)

(*ω*)*v*) (*j* = 1, . . . , *L*). (22)

), (23)

*ζ* = 0, (17)

That is, decompose *X* such that each ODE, *ω*˙ = *X*[*j*]

*X*[1]

*X*[2]

*X*[3]

*X*[4]

*φ*[1] (*t*) = 

*X*�[*j*] (*ω*�

Although, in a strict sense, regarding a solution Φ[*j*] of

(*ω*) ≡

This is indeed the solvable decomposition as is easily seen; e.g., *ω*˙ = *X*[1]

*<sup>x</sup>*˙ <sup>=</sup> *<sup>p</sup>* · **<sup>M</sup>**<sup>−</sup>1, *<sup>p</sup>*˙ <sup>=</sup> 0, ˙

so we have its solution *<sup>φ</sup>*[1] with an initial value *<sup>ω</sup>*<sup>0</sup> = (*x*0, *<sup>p</sup>*0, *<sup>ζ</sup>*0) in an explicit form for all

(ii) Consider a solvable decomposition for, in turn, the extended ODE (2), which is now

To obtain an extended phase-space volume-preserving integrator we should devise a

decomposition *X*� = ∑ *X*�[*j*] that yields, for each *j*, the divergence-free property,

Fortunately, this can be done automatically; we have a desired decomposition:

)=(*X*[*j*]

In fact, for the NH equation, according to Eqs. (16) and (22) we have

*X*�[1] (*ω*� ) = 

*X*�[2] (*ω*�

*X*�[3] (*ω*�

*X*�[4] (*ω*� ) = 0, <sup>−</sup> *<sup>ζ</sup>*

*<sup>X</sup>*� <sup>=</sup> *<sup>X</sup>*�[1] <sup>+</sup> *<sup>X</sup>*�[2] <sup>+</sup> ··· <sup>+</sup> *<sup>X</sup>*�[*L*]

(*ω*), <sup>−</sup>div*X*[*j*]

*ω*˙ � = *X*�[*j*]

its *v*-component is generally given via a form that is integrated with respect to time, this form can be evaluated explicitly in many cases as already discussed (Fukuda & Nakamura, 2006).

(*ω*�

*<sup>p</sup>* · **<sup>M</sup>**<sup>−</sup>1, 0, 0, 0

*<sup>Q</sup> <sup>p</sup>*, 0, *<sup>n</sup><sup>ζ</sup> Q v* 

) = (0, −∇*U*(*x*),0,0), (24b)

) = (0, 0, 2*K*(*p*) − *nk*B*T*, 0), (24c)

*tp*<sup>0</sup> · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>*0, *<sup>p</sup>*0, *<sup>ζ</sup>*<sup>0</sup>

∑4

time *t*, as

represented by

*<sup>j</sup>*=<sup>1</sup> *<sup>X</sup>*[*j*] defined as:

where *ω*� ≡ (*x*, *p*, *ζ*, *v*); and a solution with an initial value *ω*� <sup>0</sup> ≡ (*x*0, *p*0, *ζ*0, *v*0) for ODE (23) is given, respectively, by

$$\Phi^{[1]}(t) = \left(tp\_0 \cdot \mathbf{M}^{-1} + \mathbf{x}\_{0\prime} \, p\_0 \, \zeta\_{0\prime} \, v\_0\right),\tag{25a}$$

$$\Phi^{[2]}(t) = \left(\mathbf{x}\_{0\prime} - t\nabla U(\mathbf{x}\_0) + p\_{0\prime}\mathcal{L}\_{0\prime}v\_0\right),\tag{25b}$$

$$\Phi^{[3]}(t) = \left(\mathbf{x}\_{0\prime} p\_{0\prime} \left(2\mathbf{K}(p\_0) - nk\_{\mathbf{B}}T\right)t + \zeta\_{0\prime}v\_0\right),\tag{25c}$$

$$\Phi^{[4]}(t) = \left(\mathbf{x}\_{0\prime}\mathbf{e}^{-\frac{\zeta\_0 t}{Q}}\boldsymbol{p}\_{0\prime}\boldsymbol{\zeta}\_0, \mathbf{e}^{\frac{\kappa\_0 t}{Q}}\boldsymbol{v}\_0\right). \tag{25d}$$

(iii) Turn attention from solutions to phase-space maps.

Consider any decomposed component *j*. Even when we rewrite Φ[*j*] (*t*) as Φ[*j*]*ω*� <sup>0</sup> (*t*) to express the initial value, Φ[*j*]*ω*� <sup>0</sup> (*t*) still denotes the solution as long as we focus on the variation with respect to time *t* by fixing initial value *ω*� <sup>0</sup>. On the other hand, if we focus on the variation with respect to initial value by fixing time, we get a phase-space map Φ[*j*] *<sup>t</sup>* : Ω� → Ω� as,

$$\Phi\_t^{[j]}(\omega') \equiv \Phi^{[j]\omega'}(t),\tag{26}$$

where we have dropped the suffix "0" in the initial value because we consider all the (extended) phase-space point to be impartial. Namely Φ[*j*] *<sup>t</sup>* maps every point *ω*� in Ω� to the point reached by the solution with initial value *ω*� for ODE (23) after a period *t*. Volume-preserving property is ensured by Eq. (20); i.e., for any area *A* ⊂ Ω� , its volume and the volume of the mapped area are equal:

$$\text{Vol}(\Phi\_t^{[j]}(A)) = \text{Vol}(A). \tag{27}$$

For the NH equation, it follows from Eq. (25) that e.g., Φ[1] *<sup>t</sup>* maps *ω*� = (*x*, *p*, *ζ*, *v*) to *tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>*, *<sup>p</sup>*, *<sup>ζ</sup>*, *<sup>v</sup>* ; here, the variable that changes is only *x*, and in this way we can express the essential changes provided by each map Φ[*j*] *<sup>t</sup>* as

*<sup>x</sup>* <sup>→</sup> *tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>* by <sup>Φ</sup>[1] *<sup>t</sup>* , (28a)

$$p \rightarrow -t\nabla \mathcal{U}(\mathbf{x}) + p \text{ by } \Phi\_t^{[2]},\tag{28b}$$

$$\zeta \to (2K(p) - nk\_{\mathbb{B}}T)\mathfrak{t} + \zeta \text{ by } \Phi\_{\mathfrak{t}}^{[3]}\text{.}\tag{28c}$$

$$p \to \mathbf{e}^{-(\mathbb{L}/\mathbb{Q})t} p\_\prime \, v \to \mathbf{e}^{(n\mathbb{L}/\mathbb{Q})t} v \text{ by } \Phi\_t^{[4]}.\tag{28d}$$

(iv) Combine the individual maps as

$$\Phi\_t \equiv \Phi\_t^{[1]} \circ \Phi\_t^{[2]} \circ \cdots \circ \Phi\_t^{[L]} : \Omega' \to \Omega', \tag{29}$$

giving Φ[1] *t* Φ[2] *t* ··· <sup>Φ</sup>[*L*] *<sup>t</sup>* (*ω*� )··· for point *<sup>ω</sup>*�

Through just this procedure we get a first-order integrator for the extended ODE (19). Here, "first-order" means that the map approximates the exact flow of the ODE locally with an order of *t* 1: Φ*t*(*ω*� ) − *Tt*(*ω*� ) = *O*(*t* <sup>2</sup>) as *<sup>t</sup>* <sup>→</sup> 0. In general, map <sup>Φ</sup>*<sup>t</sup>* is said to be order *<sup>p</sup>* if

.

$$
\Phi\_t(\omega') - T\_t(\omega') = O(t^{p+1}) \text{ as } t \to 0. \tag{30}
$$

(v) Combine the maps furthermore to get a higher order integrator.

The simplest manner to do this is, by using the adjoint map

$$\Phi\_t^\* \equiv \left(\Phi\_{-t}\right)^{-1} = \Phi\_t^{[L]} \circ \cdots \circ \Phi\_t^{[2]} \circ \Phi\_t^{[1]} \,. \tag{31}$$

to adopt

$$\Psi\_{\hbar} = \Phi\_{\hbar/2} \circ \Phi\_{\hbar/2'}^\* \tag{32}$$

Therefore, the total straightforward implementation of the algorithms is described as follows:

<sup>51</sup> Numerical Integration Techniques Based on a

⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

. (38)

*x*, *p*, *ζ*, and *v* are the variables at time = *mh*

Short (or very short) time behaviour for the numerical integration is governed by Equation (30). However, long time behaviour is beyond the reach of this equation and not so a trivial matter, especially in non-Hamiltonian systems. Since the results of the fundamental numerical investigations are already described (for the NH equation, see e.g., Queyroy et al., 2009), we here focus on the above issue. To do this, we studied the integrator for the NH equation, Eq. (11), by applying it to a one-dimensional model system, using a double-well potential. In the numerical simulation, we set *k*B*T* = 1 and *Q* = 1, with mass *mi* being unity,

Figure 3 shows the trajectories (time development) of invariant (14) obtained by the integrators with time step *h* = 0.01. It exhibits the behaviours during a relatively short period. The integrators are *P*2*S*1 [equation (32)], its higher-order version, *P*4*S*5, and the local refinement second-order scheme (LR2). Local deviations of each trajectory seem to be similar in the second order schemes, but the critical difference of trajectories between them is in the behaviour of abrupt jumps. The amplitude of such a jump is highest in the *P*2*S*1 scheme, and the amplitude, as well as the frequency, of the jumps is smaller in the other schemes. The local refinement scheme is superior to the *P*2*S*1 scheme. The most accurate integrator is *P*4*S*5, but it needs five times more evaluations of the force function than those by the second order

Long time (10<sup>7</sup> time) trajectories with time step *h* = 0.01 are shown in figure 4. In the conservation behaviour of the invariant, the "jump" appearing in a certain time scale may characterize the essential "deviation" in *larger* time scale. We would like to briefly discuss this view using figures 3 and 4, where the most visible behaviours supporting this view is observed in *P*2*S*1. Each of the jumps in a time scale of about 10<sup>4</sup> observed in figure 3 appears to correspond to an "ordinary" deviation in a time scale of 107 in figure 4. The direction (viz., up or down) of the deviation, in general, seems to be random. However, sometimes, consecutive deviations in a same direction occur, and they eventually may form a jump. Such a relation, viz., the directed consecutive deviations form a jump and the resulting jumps characterize the deviations in a larger time scale, might continue in further larger time scale. The true origin of the jumps in the invariant is not clear, requiring further investigations for longer time simulations. However, the jump would not be directly induced by the "actual" jump of coordinate *x* between the two wells of the one-dimensional double-well potential, nevertheless the jump of *x* induces large motion in the phase space. This is because the jump of *x* with our schemes occurred too frequently (but plausible in the sense of producing the

Boltzmann-Gibbs distribution), compared with the invariant jump.

Do *m* = 1, total time step

*x* = *x*, *p* = *p*, *ζ* = *ζ*, *v* = *v*

Do *k* = 1,*s t* = *βkh* call FAIAD(*t*) *t* = *αkh* call FAI(*t*)

Geometric View and Application to Molecular Dynamics Simulations

enddo

enddo

and initial values were *x*<sup>0</sup> = 0, *p*<sup>0</sup> = 1, *ζ*<sup>0</sup> = 0, *v*<sup>0</sup> = 1.

**2.3 Numerical simulation**

schemes.

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

where *h* stands for a time step. Equation (32) is shown to be a second order integrator, which is a generalized version of the Verlet integrator. In fact, for the NH equation, it can be easily seen from Eq. (28) that Eq. (32) gives the (position) Verlet integrator (Tuckerman et al., 1992) when one ignores the non-canonical variables *ζ* and *v*.

Reaching a construction of first order integrators (29) and (31) enables us to use several general mathematical schemes to obtain higher order integrators (Hairer et al., 2002; McLachlan, 1995). Using a number of parameters, which designate the division of time step *h*, a more general formula for higher order method can be

$$
\Psi\_{\rm li} = \Phi\_{\rm a\_{s}h} \circ \Phi\_{\beta\_{s}h}^{\*} \circ \cdots \circ \Phi\_{\rm a\_{2}h} \circ \Phi\_{\beta\_{2}h}^{\*} \circ \Phi\_{\rm a\_{1}h} \circ \Phi\_{\beta\_{1}h}^{\*}.\tag{33}
$$

We assume the symmetric condition *<sup>α</sup><sup>i</sup>* <sup>=</sup> *<sup>β</sup>s*+1−*<sup>i</sup>* to get the symmetric property:

$$(\Psi\_h)^{-1} = \Psi\_{-h}.\tag{34}$$

Namely, the numerical map is time-reversible as is the exact flow. Since this property is the most fundamental and universal property observed in any smooth ODE, this preservation in a numerical map is highly recommended. Volume-preserving property is achieved by the fact that each map is volume preserving due to Eq. (27):

$$\text{Vol}(\mathbb{Y}\_h(A)) = \text{Vol}(A). \tag{35}$$

Finally, we give the whole explicit algorithm for the NH equation. Equation (31) stands for the mapping from (*x*, *p*, *ζ*, *v*) to (*x*, *p*, *ζ*, *v*), summarized in the following procedure:

$$\text{FALAD}(t): \begin{cases} \mathbf{x} \rightarrow \overline{\mathbf{x}} \equiv tp \cdot \mathbf{M}^{-1} + \mathbf{x} & \text{by } \Phi\_t^{[1]} \\ p \rightarrow \vec{p} \equiv -t\nabla U(\overline{\mathbf{x}}) + p & \text{by } \Phi\_t^{[2]} \\ \updownarrow \rightarrow \overline{\mathfrak{J}} \equiv (2\mathcal{K}(\vec{p}) - nk\_\mathcal{B}\mathbf{T})t + \zeta & \text{by } \Phi\_t^{[3]} \\ \vec{p} \rightarrow \overline{p} \equiv \mathbf{e}^{-(\overline{\mathcal{L}}/Q)t} \not p, \ v \rightarrow \overline{v} \equiv \mathbf{e}^{(n\overline{\mathcal{L}}/Q)t} v & \text{by } \Phi\_t^{[4]} \end{cases} \tag{36}$$

Similarly, Eq. (29) does from (*x*, *p*, *ζ*, *v*) to (*x*, *p*, *ζ*, *v*), summarized in the procedure,

$$\text{FAI}(t): \begin{cases} \overline{p} \to \frac{\widetilde{\mathbf{p}}}{\widetilde{p}} \equiv \text{e}^{-(\widetilde{\zeta}/Q)t} \overline{p}, \ \overline{v} \to \frac{\overline{\mathbf{v}}}{\widetilde{v}} \equiv \text{e}^{(n\widetilde{\zeta}/Q)t} \overline{v} \text{ by } \Phi\_{t}^{[4]}.\\ \overline{\zeta} \to \overline{\overline{\zeta}} \equiv (2K(\overline{\widetilde{p}}) - nk\_{\mathbb{B}}T)t + \overline{\zeta} & \quad \text{by } \Phi\_{t}^{[3]}.\\ \overline{\overline{p}} \to \overline{\overline{p}} \equiv -t\nabla\mathcal{U}(\overline{\boldsymbol{x}}) + \overline{\overline{p}} & \quad \text{by } \Phi\_{t}^{[2]}.\\ \overline{\boldsymbol{x}} \to \overline{\overline{\boldsymbol{x}}} \equiv t\overline{\overline{p}} \cdot \mathbf{M}^{-1} + \overline{\boldsymbol{x}} & \quad \text{by } \Phi\_{t}^{[1]}.\end{cases} \tag{37}$$

⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

Therefore, the total straightforward implementation of the algorithms is described as follows:

$$\begin{array}{l}\text{Do } m = 1, \text{ total time step} \\ \text{Do } k = 1, \text{'s} \\ t = \beta\_k h \\ \text{call FAAD}(t) \\ t = a\_k h \\ \text{call FAI}(t) \\ x = \overline{\overline{\zeta}}, \; p = \overline{\overline{\overline{\nu}}}, \; \zeta = \overline{\overline{\zeta}}, \; v = \overline{\overline{\upsilon}} \\ \text{enddo} \\ x\_\circ p\_\circ \overline{\zeta}\_\circ \text{ and } v \text{ are the variables at time} = mh \\ \text{enddo} \end{array} \tag{38}$$

#### **2.3 Numerical simulation**

8 Will-be-set-by-IN-TECH

<sup>−</sup><sup>1</sup> = Φ[*L*]

Ψ*<sup>h</sup>* = Φ*h*/2 ◦ Φ<sup>∗</sup>

where *h* stands for a time step. Equation (32) is shown to be a second order integrator, which is a generalized version of the Verlet integrator. In fact, for the NH equation, it can be easily seen from Eq. (28) that Eq. (32) gives the (position) Verlet integrator (Tuckerman et al., 1992)

Reaching a construction of first order integrators (29) and (31) enables us to use several general mathematical schemes to obtain higher order integrators (Hairer et al., 2002; McLachlan, 1995). Using a number of parameters, which designate the division of time step *h*, a more

*<sup>β</sup>s<sup>h</sup>* ◦···◦ <sup>Φ</sup>*α*2*<sup>h</sup>* ◦ <sup>Φ</sup><sup>∗</sup>

Namely, the numerical map is time-reversible as is the exact flow. Since this property is the most fundamental and universal property observed in any smooth ODE, this preservation in a numerical map is highly recommended. Volume-preserving property is achieved by the fact

Finally, we give the whole explicit algorithm for the NH equation. Equation (31) stands for

*<sup>p</sup>*˜ <sup>→</sup> *<sup>p</sup>* <sup>≡</sup> <sup>e</sup>−(*ζ*/*Q*)*<sup>t</sup> <sup>p</sup>*˜, *<sup>v</sup>* <sup>→</sup> *<sup>v</sup>* <sup>≡</sup> <sup>e</sup>(*nζ*/*Q*)*<sup>t</sup>*

*<sup>p</sup>* <sup>→</sup> �*<sup>p</sup>* <sup>≡</sup> <sup>e</sup>−(*ζ*/*Q*)*<sup>t</sup> <sup>p</sup>*, *<sup>v</sup>* <sup>→</sup> *<sup>v</sup>* <sup>≡</sup> <sup>e</sup>(*nζ*/*Q*)*<sup>t</sup>*

*<sup>ζ</sup>* <sup>→</sup> *<sup>ζ</sup>* <sup>≡</sup> (2*K*(�*p*) <sup>−</sup> *nk*B*T*)*<sup>t</sup>* <sup>+</sup> *<sup>ζ</sup>* by <sup>Φ</sup>[3]

�*<sup>p</sup>* <sup>→</sup> *<sup>p</sup>* ≡ −*t*∇*U*(*x*) + �*<sup>p</sup>* by <sup>Φ</sup>[2]

*<sup>x</sup>* <sup>→</sup> *<sup>x</sup>* <sup>≡</sup> *tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>* by <sup>Φ</sup>[1]

*<sup>x</sup>* <sup>→</sup> *<sup>x</sup>* <sup>≡</sup> *tp* · **<sup>M</sup>**−<sup>1</sup> <sup>+</sup> *<sup>x</sup>* by <sup>Φ</sup>[1]

*<sup>p</sup>* <sup>→</sup> *<sup>p</sup>*˜ ≡ −*t*∇*U*(*x*) + *<sup>p</sup>* by <sup>Φ</sup>[2]

*<sup>ζ</sup>* <sup>→</sup> *<sup>ζ</sup>* <sup>≡</sup> (2*K*(*p*˜) <sup>−</sup> *nk*B*T*)*<sup>t</sup>* <sup>+</sup> *<sup>ζ</sup>* by <sup>Φ</sup>[3]

We assume the symmetric condition *<sup>α</sup><sup>i</sup>* <sup>=</sup> *<sup>β</sup>s*+1−*<sup>i</sup>* to get the symmetric property:

(Ψ*h*)

the mapping from (*x*, *p*, *ζ*, *v*) to (*x*, *p*, *ζ*, *v*), summarized in the following procedure:

Similarly, Eq. (29) does from (*x*, *p*, *ζ*, *v*) to (*x*, *p*, *ζ*, *v*), summarized in the procedure,

*<sup>t</sup>* ◦···◦ <sup>Φ</sup>[2]

*<sup>t</sup>* ◦ <sup>Φ</sup>[1]

*<sup>β</sup>*2*<sup>h</sup>* ◦ <sup>Φ</sup>*α*1*<sup>h</sup>* ◦ <sup>Φ</sup><sup>∗</sup>

<sup>−</sup><sup>1</sup> <sup>=</sup> <sup>Ψ</sup>−*h*. (34)

*t* ,

⎫ ⎪⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎪⎭

(36)

(37)

*t* ,

*t* ,

*v* by Φ[4] *t* .

*v* by Φ[4] *t* ,

*t* ,

⎫ ⎪⎪⎪⎪⎪⎬

⎪⎪⎪⎪⎪⎭

*t* ,

*t* .

Vol(Ψ*h*(*A*)) = Vol(*A*). (35)

*<sup>t</sup>* , (31)

*<sup>β</sup>*1*h*. (33)

*<sup>h</sup>*/2, (32)

(v) Combine the maps furthermore to get a higher order integrator.

*<sup>t</sup>* ≡ (Φ−*t*)

The simplest manner to do this is, by using the adjoint map

Φ∗

when one ignores the non-canonical variables *ζ* and *v*.

Ψ*<sup>h</sup>* = Φ*αsh* ◦ Φ<sup>∗</sup>

general formula for higher order method can be

that each map is volume preserving due to Eq. (27):

⎧ ⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

⎧ ⎪⎪⎪⎪⎪⎨

⎪⎪⎪⎪⎪⎩

FAIAD(*t*) :

FAI(*t*) :

to adopt

Short (or very short) time behaviour for the numerical integration is governed by Equation (30). However, long time behaviour is beyond the reach of this equation and not so a trivial matter, especially in non-Hamiltonian systems. Since the results of the fundamental numerical investigations are already described (for the NH equation, see e.g., Queyroy et al., 2009), we here focus on the above issue. To do this, we studied the integrator for the NH equation, Eq. (11), by applying it to a one-dimensional model system, using a double-well potential. In the numerical simulation, we set *k*B*T* = 1 and *Q* = 1, with mass *mi* being unity, and initial values were *x*<sup>0</sup> = 0, *p*<sup>0</sup> = 1, *ζ*<sup>0</sup> = 0, *v*<sup>0</sup> = 1.

Figure 3 shows the trajectories (time development) of invariant (14) obtained by the integrators with time step *h* = 0.01. It exhibits the behaviours during a relatively short period. The integrators are *P*2*S*1 [equation (32)], its higher-order version, *P*4*S*5, and the local refinement second-order scheme (LR2). Local deviations of each trajectory seem to be similar in the second order schemes, but the critical difference of trajectories between them is in the behaviour of abrupt jumps. The amplitude of such a jump is highest in the *P*2*S*1 scheme, and the amplitude, as well as the frequency, of the jumps is smaller in the other schemes. The local refinement scheme is superior to the *P*2*S*1 scheme. The most accurate integrator is *P*4*S*5, but it needs five times more evaluations of the force function than those by the second order schemes.

Long time (10<sup>7</sup> time) trajectories with time step *h* = 0.01 are shown in figure 4. In the conservation behaviour of the invariant, the "jump" appearing in a certain time scale may characterize the essential "deviation" in *larger* time scale. We would like to briefly discuss this view using figures 3 and 4, where the most visible behaviours supporting this view is observed in *P*2*S*1. Each of the jumps in a time scale of about 10<sup>4</sup> observed in figure 3 appears to correspond to an "ordinary" deviation in a time scale of 107 in figure 4. The direction (viz., up or down) of the deviation, in general, seems to be random. However, sometimes, consecutive deviations in a same direction occur, and they eventually may form a jump. Such a relation, viz., the directed consecutive deviations form a jump and the resulting jumps characterize the deviations in a larger time scale, might continue in further larger time scale. The true origin of the jumps in the invariant is not clear, requiring further investigations for longer time simulations. However, the jump would not be directly induced by the "actual" jump of coordinate *x* between the two wells of the one-dimensional double-well potential, nevertheless the jump of *x* induces large motion in the phase space. This is because the jump of *x* with our schemes occurred too frequently (but plausible in the sense of producing the Boltzmann-Gibbs distribution), compared with the invariant jump.

Fig. 3. Trajectories of the invariant using the current integration scheme, applied to the NH equation of the one dimensional system.

Fig. 4. Long time behaviours of the invariant shown in Figure 3.

systematic efficient scheme also in such a situation.

necessarily require the potential function.

the assumption itself for the validness of the Liouville equation in our method can actually be removed if we employ the twisting technique (Fukuda & Nakamura, 2006). Using this technique or its variant, we can apply our method to the system that is not Liouvillian. For example, non-equilibrium problems, including the heat flow problem, are important targets. As an alternative method to control the temperature, the equation developed by Berendsen et al. (Berendsen et al., 1984), which is not Liouvillian to the best of our knowledge, is also a target of our method. Berendsen's method is simple, stable, and has been employed by many researchers to perform e.g., biomolecular simulations. Our protocol is expected to introduce a

<sup>53</sup> Numerical Integration Techniques Based on a

Geometric View and Application to Molecular Dynamics Simulations

Note that our protocol is also applicable when the force function defined in the equations of motion does not correspond to an explicit potential function, such as in the case where a tabulated force is used (Queyroy et al., 2004). Then, the numerical integration error check using the ordinary MD protocol (*NVE*, *NTV*, etc.) is not sufficient, since the invariants always use the potential function. However, even in such a situation, our protocol offers a solution for the consistent numerical integration error check. This is because the function *B* (in the definition of the invariant; see Eq. (4)) is an arbitrary function in general and so does not

Regarding the technical issues, we note two points. The first point is the multiple extended-variable formalism developed by Queyroy et al. (Queyroy et al., 2009). Extended

### **3. Discussions**

We have utilized the NH equation to explain the integration techniques. The NH equation is a representative of the Boltzmann-Gibbs dynamics, i.e., an ODE that can directly generate the Boltzmann-Gibbs distribution or *NTV* ensemble (Hoover & Holian, 1996). For such a dynamics, the Nosé-Hoover chain (NHC) equation (Martyna et al. 1992), the Kusnezov, Bulgac, and Bauer equation (Kusnezov et al. 1990), and the generalized Gaussian moment thermostat equation (Liu & Tuckerman, 2000) have been developed, along the line of the generalization of the NH method. Martyna et al. (Martyna et al. 1996) have developed the numerical algorithms to integrate the NHC equations and generalized the scheme to *NTP* ensemble. The current method described so far is also applicable to *NTV* equations other than the NH equation and different ensembles, e.g., *NTP* protocol, including the non-Hamiltonian method derived by Melchionna et al. (Melchionna et al., 1993), and the generalized ensemble (Fukuda & Nakamura, 2002) for the Tsallis statistics (Tsallis, 1988). These targets are well suited to our method, since we have assumed the Liouville equation with *B*≡ − ln *ρ* and since they have densities satisfying the Liouville equation.

The Liouville equation is valid in many MD equations. As well as the Liouvillian case, our scheme is applicable to the locally Liouvillian case such as the Gaussian isokinetic equation (Hoover et al., 1982; Evans, 1983; Evans et al., 1983; Evans & Morris, 1983). However, 10 Will-be-set-by-IN-TECH

Fig. 3. Trajectories of the invariant using the current integration scheme, applied to the NH

We have utilized the NH equation to explain the integration techniques. The NH equation is a representative of the Boltzmann-Gibbs dynamics, i.e., an ODE that can directly generate the Boltzmann-Gibbs distribution or *NTV* ensemble (Hoover & Holian, 1996). For such a dynamics, the Nosé-Hoover chain (NHC) equation (Martyna et al. 1992), the Kusnezov, Bulgac, and Bauer equation (Kusnezov et al. 1990), and the generalized Gaussian moment thermostat equation (Liu & Tuckerman, 2000) have been developed, along the line of the generalization of the NH method. Martyna et al. (Martyna et al. 1996) have developed the numerical algorithms to integrate the NHC equations and generalized the scheme to *NTP* ensemble. The current method described so far is also applicable to *NTV* equations other than the NH equation and different ensembles, e.g., *NTP* protocol, including the non-Hamiltonian method derived by Melchionna et al. (Melchionna et al., 1993), and the generalized ensemble (Fukuda & Nakamura, 2002) for the Tsallis statistics (Tsallis, 1988). These targets are well suited to our method, since we have assumed the Liouville equation with *B*≡ − ln *ρ* and since

The Liouville equation is valid in many MD equations. As well as the Liouvillian case, our scheme is applicable to the locally Liouvillian case such as the Gaussian isokinetic equation (Hoover et al., 1982; Evans, 1983; Evans et al., 1983; Evans & Morris, 1983). However,

equation of the one dimensional system.

they have densities satisfying the Liouville equation.

**3. Discussions**

Fig. 4. Long time behaviours of the invariant shown in Figure 3.

the assumption itself for the validness of the Liouville equation in our method can actually be removed if we employ the twisting technique (Fukuda & Nakamura, 2006). Using this technique or its variant, we can apply our method to the system that is not Liouvillian. For example, non-equilibrium problems, including the heat flow problem, are important targets. As an alternative method to control the temperature, the equation developed by Berendsen et al. (Berendsen et al., 1984), which is not Liouvillian to the best of our knowledge, is also a target of our method. Berendsen's method is simple, stable, and has been employed by many researchers to perform e.g., biomolecular simulations. Our protocol is expected to introduce a systematic efficient scheme also in such a situation.

Note that our protocol is also applicable when the force function defined in the equations of motion does not correspond to an explicit potential function, such as in the case where a tabulated force is used (Queyroy et al., 2004). Then, the numerical integration error check using the ordinary MD protocol (*NVE*, *NTV*, etc.) is not sufficient, since the invariants always use the potential function. However, even in such a situation, our protocol offers a solution for the consistent numerical integration error check. This is because the function *B* (in the definition of the invariant; see Eq. (4)) is an arbitrary function in general and so does not necessarily require the potential function.

Regarding the technical issues, we note two points. The first point is the multiple extended-variable formalism developed by Queyroy et al. (Queyroy et al., 2009). Extended

[6] G. S. Ezra (2006). Reversible measure-preserving integrators for non-Hamiltonian

<sup>55</sup> Numerical Integration Techniques Based on a

[7] I. Fukuda & H. Nakamura (2002). Tsallis dynamics using the Nosé-Hoover approach,

[8] I. Fukuda & H. Nakamura (2004). Efficiency in the generation of the Boltzmann-Gibbs distribution by the Tsallis dynamics reweighting method, *J. Phys. Chem. B* 108: 4162-4170. [9] I. Fukuda & H. Nakamura (2006). Construction of an extended invariant for an arbitrary ordinary differential equation with its development in a numerical integration algorithm,

[10] E. Hairer, C. Lubich, & G. Wanner (2002). *Geometric Numerical Integration*,

[11] W. G. Hoover, A. J. C. Ladd, & B. Moran (1982). High-strain-rate plastic flow studied via

[12] W. G. Hoover (1985). Canonical dynamics: Equilibrium phase-space distributions, *Phys.*

[14] W. G. Hoover & B. L. Holian (1996). Kinetic moments method for the canonical ensemble

[16] F. Kang & S. Zai-jiu (1995). Volume-preserving algorithms for source-free dynamical

[17] D. Kusnezov, A. Bulgac, & W. Bauer (1990). Canonical ensembles from chaos, *Ann. Phys.*

[18] Y. Liu & M. E. Tuckerman (2000). Generalized Gaussian moment thermostatting: A new continuous dynamical approach to the canonical ensemble, *J. Chem. Phys.* 112: 1685-1700. [19] G. J. Martyna, M. L. Klein, & M. Tuckerman (1992). Nosé-Hoover chains: The canonical

[20] G. J. Martyna, M. E. Tuckerman, D. J. Tobias, & M. L. Klein (1996). Explicit reversible

[21] R. I. McLachlan & P. Atela (1992). The accuracy of symplectic integrators, *Nonlinearity* 5:

[22] R. I. McLachlan (1995). On the numerical integration of ordinary differential equations

[23] S. Melchionna, G. Ciccotti, & B. L. Holian (1993). Hoover NPT dynamics for systems

[24] S. Nosé (1984). A unified formulation of the constant temperature molecular dynamics

[25] S. Nosé (1991). Constant temperature molecular dynamics methods, *Prog. Theor. Phys.*

[26] D. I. Okunbor (1995). Energy conserving, Liouville, and symplectic integrators, *J. Comp.*

[27] S. Queyroy, S. Neyertz, D. Brown, & F. Müller-Plathe (2004). Preparing relaxed systems of amorphous polymers by multiscale simulation: Application to cellulose, *Macromolecules*

[28] S. Queyroy, H. Nakamura & I. Fukuda (2009). Numerical examination of the extended phase-space volume-preserving integrator by the Nosé-Hoover molecular dynamics

nonequilibrium molecular dynamics, *Phys. Rev. Lett.* 48: 1818-1820.

[13] Wm. G. Hoover (1991). *Computational Statistical Mechanics,* Elsevier, N.Y.

[15] D. Husemoller (1966). *Fibre Bundles, second edition*, Springer, N.Y.

ensemble via continuous dynamics, *J. Chem. Phys.* 97: 2635.

varying in shape and size, *Mol. Phys.* 78: 533-544.

methods, *J. Chem. Phys.* 81: 511-519.

equations, *J. Comput. Chem*. 30: 1799-1815.

integrators for extended systems dynamics, *Mol. Phys*. 87: 1117-1157.

by symmetric composition methods, *SIAM J. Sci. Comput*. 16: 151-168.

systems, *J. Chem. Phys*. 125: 034104.

Geometric View and Application to Molecular Dynamics Simulations

*Phys. Rev. E* 65: 026105.

*Phys. Rev. E* 73: 026703.

Springer-Verlag, Berlin.

*Rev. A* 31: 1695-1697.

*(N.Y.)* 204: 155-185.

541-562.

*Suppl.* 103: 1-46.

37: 7338-7350.

*Phys.* 120: 375-378.

distribution, *Phys. Lett. A* 211: 253-257.

systems, *Numer. Math.* 71: 451-463.

phase-space volume-preserving integrator should keep the volume of the extended space spanned by the original variables *ω*1, ...., *ω<sup>N</sup>* and the extended variable *v*. Sometimes, the subspace volume change derived by the former variables is very large, which is emphasized as the system size *n* is large. In such a case, to compensate the large change, the single extended variable *v* takes a very large or very small absolute value, as seen in the mapping *<sup>v</sup>* <sup>→</sup> <sup>e</sup>(*nζ*/*Q*)*<sup>t</sup> v* by Φ[4] *<sup>t</sup>* n the NH equation. This situation can be avoided by the multiple extended-variable formalism, where *v*1, ...., *vM* are used instead of *v* (practically *M* = *n*).

The second point is the order of appearance of each map for constructing the basic first order integrator. The first order integrator, Eq. (29), has an arbitrariness for the permutation of Φ[1] *<sup>t</sup>* , <sup>Φ</sup>[2] *<sup>t</sup>* , ··· , <sup>Φ</sup>[*L*] *<sup>t</sup>* , since the decomposition (15) is independent of the order of appearance of each vector field. Regarding the local accuracy, any permutation gives the same order of accuracy in a mathematical sense. However, other properties may be distinguishable in general (e.g., historically, the difference between the velocity Verlet and the position Verlet was discussed). For instance, the long-time behaviour and the robustness would be changed by the permutation. Further, they also depend on the physical system (particle interaction, the number of degrees of freedom, etc.) and the equations of motion (ensembles, system parameters including temperature and pressure, etc.). Thus we require more studies to clarify these issues. The cost of the force evaluation is critical in MD, and the number of force evaluations during a unit period should be small. This number depends on the permutation, in general (but the maps (25) in the NH equation is not the case, and the number is always once in the unit map Φ*h*<sup>1</sup> ◦ Φ<sup>∗</sup> *<sup>h</sup>*<sup>2</sup> in Eq. (33)). Thus the permutation yielding a small number is preferable.

So far, we have considered the volume-preserving property of the integrator in the extended phase space. Considering further the measure-preserving property in the original phase space will be of great value to the improvement of the integration. Moreover, the difficulty we should surmount is the resonance phenomena (see e.g., Schlick, 2006), which disturbs the use of larger unit time step.

#### **4. Acknowledgments**

This research was supported by Research and Development of the Next-Generation Integrated Simulation of Living Matter, a part of the Development and Use of the Next-Generation Supercomputer Project of the Ministry of Education, Culture, Sports, Science and Technology of Japan.

#### **5. References**


12 Will-be-set-by-IN-TECH

phase-space volume-preserving integrator should keep the volume of the extended space spanned by the original variables *ω*1, ...., *ω<sup>N</sup>* and the extended variable *v*. Sometimes, the subspace volume change derived by the former variables is very large, which is emphasized as the system size *n* is large. In such a case, to compensate the large change, the single extended variable *v* takes a very large or very small absolute value, as seen in the mapping

extended-variable formalism, where *v*1, ...., *vM* are used instead of *v* (practically *M* = *n*).

The second point is the order of appearance of each map for constructing the basic first order integrator. The first order integrator, Eq. (29), has an arbitrariness for the permutation of

of each vector field. Regarding the local accuracy, any permutation gives the same order of accuracy in a mathematical sense. However, other properties may be distinguishable in general (e.g., historically, the difference between the velocity Verlet and the position Verlet was discussed). For instance, the long-time behaviour and the robustness would be changed by the permutation. Further, they also depend on the physical system (particle interaction, the number of degrees of freedom, etc.) and the equations of motion (ensembles, system parameters including temperature and pressure, etc.). Thus we require more studies to clarify these issues. The cost of the force evaluation is critical in MD, and the number of force evaluations during a unit period should be small. This number depends on the permutation, in general (but the maps (25) in the NH equation is not the case, and the number is always

So far, we have considered the volume-preserving property of the integrator in the extended phase space. Considering further the measure-preserving property in the original phase space will be of great value to the improvement of the integration. Moreover, the difficulty we should surmount is the resonance phenomena (see e.g., Schlick, 2006), which disturbs the use

This research was supported by Research and Development of the Next-Generation Integrated Simulation of Living Matter, a part of the Development and Use of the Next-Generation Supercomputer Project of the Ministry of Education, Culture, Sports, Science and Technology

[1] H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola, & J. R. Haak (1984). Molecular-dynamics with coupling to an external bath*, J. Chem. Phys.* 81: 3684-3690. [2] Y. Choquet-Bruhat, C. Dewitt-Morette, & M. Dillard-Bleick (1982). *Analysis, Manifolds and*

[3] D. J. Evans (1983). Computer "experiment" for nonlinear thermodynamics of Couette

[4] D. J. Evans, W. G. Hoover, B. H. Failor, B. Moran, & A.J. C. Ladd (1983). Nonequilibrium molecular dynamics via Gauss's principle of least constraint, *Phys. Rev. A* 28: 1016-1021.

*Physics, Part I, Rev. ed.,* North-Holland, Amsterdam.

[5] D. J. Evans & G. P. Morris (1983). *Phys. Lett.* 98A: 433.

flow, *J. Chem. Phys.* 78: 3297-3302.

*<sup>t</sup>* n the NH equation. This situation can be avoided by the multiple

*<sup>t</sup>* , since the decomposition (15) is independent of the order of appearance

*<sup>h</sup>*<sup>2</sup> in Eq. (33)). Thus the permutation yielding a small number is

*<sup>v</sup>* <sup>→</sup> <sup>e</sup>(*nζ*/*Q*)*<sup>t</sup>*

*<sup>t</sup>* , ··· , <sup>Φ</sup>[*L*]

once in the unit map Φ*h*<sup>1</sup> ◦ Φ<sup>∗</sup>

of larger unit time step.

**4. Acknowledgments**

preferable.

of Japan.

**5. References**

Φ[1] *<sup>t</sup>* , <sup>Φ</sup>[2] *v* by Φ[4]


#### 14 Will-be-set-by-IN-TECH 56 Molecular Dynamics – Theoretical Developments and Applications in Nanotechnology and Energy


**0**

**4**

*México*

**Application of Molecular Dynamics**

*Facultad de Ciencias Químicas, Universidad Autónoma de Nuevo León,*

The study of chemical behavior includes answering questions as 'which isomer is the most stable?', 'which relative orientation is the most favorable for such-and-such interaction?', 'which conformer is the global minimum?', 'what are the lowest energy configurations and their relative energies?'. The answer to these questions–and many others–depends on the ability to find and study a variety of configurations of the system of interest. Recently, (Atilgan, 2007) briefly reviewed the use of molecular dynamics simulation for conformational search in the process of drug design, concluding that its use could reduce the errors in estimating binding affinities and finding more viable conformations. In addition, (Corbeil, 2009) considered the need to include ring flexibility in the conformational searches used in flexible docking. Most of the flexible docking algorithms skip searching for conformations in rings, even though a protein may stabilize a conformation other than the most stable one. The need for a tool to examine the diverse configurations of the constituent particles of a system becomes obvious even as we consider relatively small systems (10-100 atoms). Finding by hand all the conformers of cyclohexane is feasible and maybe even instructive; this is somewhat more complex when doing morpholine and even something as small as an eight-membered heterocycle can be prohibitively complex to analyze by manipulating molecular models by hand. One of the methods available to the researcher to tackle this kind of problem is Molecular Dynamics (MD) simulation. MD allows an exploration of the configurational space of a system, respecting chemical constraints. Chemical constraints–such as atomic connectivities–are needed in cases such as conformations of molecular rings and configurations of molecular clusters, e.g., solvation shells. In both cases, atomic connectivities must be kept intact, otherwise we risk breaking up the ring or the molecular constituents of the cluster. In our research group we routinely employ MD simulation, sometimes in its semiempirical variety, to study small systems (systems with fewer than 100 atoms) whether

This review is narrowly focused on current software and methods appropriate for doing MD simulations of small heterocycles and clusters composed of 10-100 molecules. In particular, we try to systematize the tools available to tackle the problem of searching for minima in heterocycles and in molecular clusters. We want to use MD simulations as a tool to explore the energy landscape of a small system, so we can locate the global minimum. We do not include the vast literature simulating water solvation by MD. Even though the aggregates of water

they be solvation shells, inorganic clusters or heterocycles.

**1. Introduction**

**Simulation to Small Systems**

*San Nicolás de los Garza, N. L.*

Víctor M. Rosas-García and Isabel Sáenz-Tavera

