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

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.


Introduction
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.
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.

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 ), viz., dω i /dt = X i (ω 1 , ...., ω N ) for i = 1, ..., N, consider an extended ODĖ which is defined, on an extended phase space where a point in Ω ′ is expressed as ω ′ =(ω, v) [R × denotes non-zero numbers; D i the partial differentiation with respect to ω i (all quantities are supposed to be sufficiently smooth)]. Here, B is any phase-space function. Then we can show that L(ω, v)=B(ω)+ln |v| (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, holds for a certain density function ρ, and put B≡−ln ρ as the choice of B. Thus, it follows from Eq. (3b) that X e (ω ′ )=−divX(ω)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, ω ′ t ≡ φ ′ (t) is a point of the solution at time t with initial value ω ′ , and the similar holds for ω t ≡ φ(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. 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 ω 0 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 45 Numerical Integration Techniques Based on a Geometric View and Application to Molecular Dynamics Simulations www.intechopen.com 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 Ω × R × . A level set (i.e., a generalization of energy surface) is represented by 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 This indicates that L ± c is equivalent, in a topological sense, to the original space Ω. This situation is different in the Hamiltonian system in which a level set {(x, p) ∈ R 2n | H(x, p)= c} often becomes compact and is not necessarily equivalent to R 2n−1 nor R 2n . 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.

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 Ψ h 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 Ψ h constitute the integration [T 0 ≡identity, T −m ≡ (T m ) −1 can be defined]. This iteration corresponds to an approximation to Exact flow ≡{ T t : Ω ′ → Ω ′ } t: real numbers ≡ Continuous dynamical system (generated by ODE).
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 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 bẏ where ω ≡ (x, p, ζ) and M ≡diag(m 1 , ..., m n ); x ≡ (x 1 ,...,x n ), p ≡ (p 1 ,...,p n ), U(x), and K(p) ≡ ∑ n i=1 p 2 i /2m i represent the coordinates, momenta, potential energy, and kinetic energy, respectively, of a physical system; k B 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 is a density for the Liouville equation.
and invariant (4) becomes Now, the five procedures are described as follows: (i) For a given ODE (1), consider a solvable decomposition

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

www.intechopen.com
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 j=1 X [j] defined as: This is indeed the solvable decomposition as is easily seen; e.g.,ω = X [1] (ω) iṡ so we have its solution φ [1] with an initial value ω 0 =( x 0 , p 0 , ζ 0 ) in an explicit form for all time t,as φ [1] (ii) Consider a solvable decomposition for, in turn, the extended ODE (2), which is now represented byω 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: Although, in a strict sense, regarding a solution Φ [j] oḟ 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).
Consider any decomposed component j. Even when we rewrite denotes the solution as long as we focus on the variation with respect to time t by fixing initial value ω ′ 0 . On the other hand, if we focus on the variation with respect to initial value by fixing time, we get a phase-space map Φ 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]
t 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: For the NH equation, it follows from Eq. (25) that e.g., Φ [1] t maps ω ′ =( x, p, ζ, v) to tp · M −1 + x, p, ζ, v ; here, the variable that changes is only x, and in this way we can express the essential changes provided by each map Φ

[j]
t as (iv) Combine the individual maps as giving Φ 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 (ω ′ ) − T t (ω ′ )=O(t 2 ) as t → 0. In general, map Φ t is said to be order p if (v) Combine the maps furthermore to get a higher order integrator.
The simplest manner to do this is, by using the adjoint map to adopt 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  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 We assume the symmetric condition α i = β s+1−i to get the symmetric property: 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): Vo l (Ψ h (A)) = Vo l (A).
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: Similarly, Eq. (29) does from (x, p, ζ, v) to (x, p, ζ, v), summarized in the procedure, Therefore, the total straightforward implementation of the algorithms is described as follows: ⎡ (38)

Numerical simulation
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 m i being unity, and initial values were  (32)], its higher-order version, P4S5, 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 P2S1 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 P2S1 scheme. The most accurate integrator is P4S5, but it needs five times more evaluations of the force function than those by the second order schemes.
Long time (10 7 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 P2S1. Each of the jumps in a time scale of about 10 4 observed in figure 3 appears to correspond to an "ordinary" deviation in a time scale of 10 7 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.

51
Numerical Integration Techniques Based on a Geometric View and Application to Molecular Dynamics Simulations www.intechopen.com

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 ), 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 & Morris, 1983). However, 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. phase-space volume-preserving integrator should keep the volume of the extended space spanned by the original variables ω 1 , ...., ω N 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 v → e (nζ/Q)t v by Φ [4] t n the NH equation. This situation can be avoided by the multiple extended-variable formalism, where v 1 , ...., v M 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 Φ t , 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 1 • Φ * h 2 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. Molecular Dynamics is a two-volume compendium of the ever-growing applications of molecular dynamics simulations to solve a wider range of scientific and engineering challenges. The contents illustrate the rapid progress on molecular dynamics simulations in many fields of science and technology, such as nanotechnology, energy research, and biology, due to the advances of new dynamics theories and the extraordinary power of today's computers. This first book begins with a general description of underlying theories of molecular dynamics simulations and provides extensive coverage of molecular dynamics simulations in nanotechnology and energy. Coverage of this book includes: Recent advances of molecular dynamics theory Formation and evolution of nanoparticles of up to 106 atoms Diffusion and dissociation of gas and liquid molecules on silicon, metal, or metal organic frameworks Conductivity of ionic species in solid oxides Ion solvation in liquid mixtures Nuclear structures

How to reference
In order to correctly reference this scholarly work, feel free to copy and paste the following: Ikuo Fukuda and Séverine Queyroy (2012). Numerical Integration Techniques Based on a Geometric View and Application to Molecular Dynamics Simulations, Molecular Dynamics -Theoretical Developments and Applications in Nanotechnology and Energy, Prof. Lichang Wang (Ed.), ISBN: 978-953-51-0443-8, InTech, Available from: http://www.intechopen.com/books/molecular-dynamics-theoretical-developments-andapplications-in-nanotechnology-and-energy/numerical-integration-techniques-based-on-a-geometric-view-inmolecular-dynamics-simulations