Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

A molecular dynamics model for symplectic intergrators

A molecular dynamics model for symplectic intergrators Mathematical Modelling of Systems 138 1-2424197IO304-282$12.00 @wets & Zeitlinger 1997, Vol. 3, NO. 4, pp. 282-296 A Molecular Dynamics Model for Symplectic Integrators M. HANKEL*, B. KARASOZEN~, P. RENTROP* AND U. SCHMI?T$ ABSTRACT Mechanical systems in the very large scale like in celestial mechanics or in the very small scale like in the molecular dynamics can be modelled without dissipation. The resulting Hamiltonian systems possess conservation properties, which are characterized with the term symplecticness. Numerical integration schemes should preserve the symplecticness. Different methods are introduced and their performance is studied for constant and variable step size. As test examples two systems from molecular dynamics are introduced. Keywords: Hamiltonian systems, molecular models, symplectic integrators. 1 INTRODUCTION TO HAMILTONIAN SYSTEMS Hamiltonian systems play a key role in applications arising from the classical, statistical and quantum dynamics, from optics and plasma physics. In principle all phenomena where dissipation is absent or can be ignored may be modelled by Hamiltonian systems. The symplectic methods developed for the Hamiltonian systems in the recent years have become already a well established part of numerical analysis (see the monographs [3], [I 11 ). But they are far from the notion of a code (package) in the sense of standard numerical integrators with automatic error control and step size change. The numerical examples given in the literature for symplectic integrators consist mostly of systems with few degrees of freedom, like the harmonic oscillator, Kepler's problem, etc. The use of symplectic integrators for large systems of equations in molecular dynamics and for systems arising by space discretization of partial differential equations is in the beginning stage now. 'Fachbereich Mathematik, Technische Hochschule Darmstadt, Germany. l~e~axtment of Mathematics, Middle East Technical University, Ankara, Turkey. *~h~sikalische Chemie I, Technische Hochschule Darmstadt, Germany. MOLECULAR DYNAMICS MODEL 283 The Hamiltonian equations are a special kind of ordinary differential equations. The initial value problem is in Hamiltonian form if its dimension D is even, D = 2d, where d denotes the number of degree of freedom, and the vector field f (x) is of the form with the 2d x 2d skew symmetric matrix J VH(x) is the gradient of the real valued Hamilton function H(x). With x = [PI, . . . , pd, 41, . . . , qdlT, the equations (11, (2) are The flow 4, of a Hamiltonian system is a symplectic transformation for each t. In general, a transformation Q in lR2d is symplectic if its Jacobian Q1(x) satisfies elT (x) JQ1(x) = J. An equivalent definition in terms of differential forms exists too. The flow of the Hamiltonian system (3) is symplectic, if the differential two-form is preserved. Differential forms provide a geometric interpretation of the symplectic- ness in terms of conservation of areas. Another preserved quantity is the Hamiltonian, i.e if p(t) and q (t) are solutions of (3) then H (p(t), q (t)) is constant for all t. In the following we will consider special types of Hamiltonian systems occurring in molecular dynamics. Molecular dynamics (MD) simulations have been proven to be a powerful tool to obtain detailed insight into the behavior of dynamical processes in isolated molecules, clusters, biomolecules, liquids and surfaces [I], [17]. The dynamics of a molecular system can be either formulated in the classical, quantum or mixed quantum-classical framework. In a classical molecular dynamics or mechanics formulation the Hamiltonian is given by 284 M. HANKEL ET AL. where M contains the masses mi of the molecules and is an invertible symmetric matrix. The first part of the Hamiltonian denotes the kinetic energy T(p), the second part the potential energy V(q). Usually the variables q; denote the generalized coordinates and pi the generalized momenta. The resulting Hamiltonian system has the form The quantum dynamics problems can be formulated as the discretized Schrodinger equation (h = 1) in matrix representation with the Hamilton operator Replacing the complex vector c(t) by a linear combination of formal momentum and position coordinates and separating real and imaginary parts, the classical Hamiltonian formulation with HQ is then obtained where Hij are the elements of the Hamilton operator % in (7) and N stands for the number of pseudo particles. By adding the classical Hamilton function Hc(p, q) governing the dynamics of the classical part of the system to the Hamilton function HQ(p, q) describing the quantum degrees of freedom the mixed quantum-classical system can be interpreted as a set of purely classical Hamilton equations of motion. The matrix elements Hij contain now the coupling between the classical and quantum modes. An important property of the Hamiltonian systems (6) is the reversibility. The system (1) is called p-reversible, if there exists an involution p, that is a linear mapping in lRW such that p2 = I with the property for all x. f (px) = -p f (x), This implies that the flow $,(x) of the system (1) is a p-reversible mapping for each t MOLECULAR DYNAMICS MODEL 285 Since 4;' = 4-, the system is also time-reversible, shortly reversible. The Hamilto- nian equations in (6) are reversible with respect to In MD simulations the number of particles necessary to obtain statistically reliable results can easily exceed N lo5, the 6N coupled first order differential equations (3) have to be integrated numerically with an efficient and robust algorithm. Therefore the need for improved algorithms has received greater attention. Recent studies have shown that symplectic techniques can be used to derive stable and efficient integrators useful in the classical and mixed quantum-classical MD simulations, and furthermore a time-scale separation can be readily overcome by introducing different time step lengths [9], [12]. Figure 1 shows a snapshot from a mixed quantum-classical MD simulation of a charged Ar,19] sir n luster performed with a multiple time step sym- plectic integrator. The charge distribution is visualized as gray scales on the atomic sites ranging from black (low charge) to white (high charge). Fig. 1. Aqls1 luster. 2 NUMERICAL HAMILTONIAN SOLVERS There exists mainly three classes of symplectic integrators; methods based on gen- erating functions, methods using compositions and symplectic Runge-Kutta methods (see the monographs [l 11, [3], Section 11.16). The generating function methods seem to be too costly. The most effective symplectic integrators are belonging to the other 286 M. HANKEL ET AL. two classes; some of them specially designed for the integration of the Hamiltonian systems (6) will be given here. Symplectic integrators using composition The separable Hamiltonian (5) provides the most common application of the com- position approach. For many systems different parts of the vector field f (x) often correspond to physically different contributions, so that a natural splitting in form of H = H1 + H2, with HI = T(p) and H2 = V(q) is available. Many Hamilto- nian systems (6) are integrable in closed form after splitting. The most widespread second order composition method is the so called VerletJStormerlleap-frog method. Application of the Verlet method to (6) gives The Verlet method is symplectic and reversible. Higher order composition methods are obtained by using the basic method qh, here the Verlet method. For example the fourth order symplectic, time-reversible composition method is given by with the coefficients cl = &. e2 = -A 2-fi. Several types of composlhon methods were derived and their properties are inves- tigated using the Baker-Campbell-Hausdorff (BCH) formula, (see [8], [14], [15]). Symplectic Runge-Kutta methods For separable Hamiltonian problems (6) there exist explicit partitioned Runge-Kutta (PRK) methods, characterized by two Butcher tableaus MOLECULAR DYNAMICS MODEL 287 where the first one is used for the p variable and the other for the q variables. The PRK method applied to the Hamiltonian system (6) results in The symplecticness condition for the irreducible PRK methods is given by ( see [l 11, Section 6.3) Remark: Hamiltonian flows exactly preserve the symplectic structure and the Hamiltonian. In general for a symplectic integrator it is impossible to preserve the Hamiltonian or the energy exactly. In practice the energy error remains constant in long terms for symplectic integrators, whereas it increases or decreases for non-symplectic integrators. This can be explained by the backward error analysis, which states that the numerical solution obtained by the symplectic integrator can be formally inter- preted as the exact solution of a perturbed nearby Hamiltonian system [4]. Symplectic integrators also have more favourable properties concerning the global error growth of the solutions when the integration is performed over very long time intervals . For integrable systems and for problems with periodic solutions, the error growth is only linear in case of symplectic integrators, compared with the quadratic error growth in the general case [2]. 3 VARIABLE STEP SIZE STRATEGIES The favourable properties of the symplectic methods for Hamiltonian systems are lost in a standard variable step size implementation based on local error estimators; they are no longer superior to the non-symplectic algorithms. Therefore the development of variable time symplectic integrators was indispensable to be competitive with the standard methods. In the following we will describe three different non-standard variable step size strategies, which preserve the geometric structure of the Hamiltonian system (6) in long time integration. There exist several variable step size strategies for reversible systems using symmetric Runge-Kutta methods. It was shown in [2], [5], [13] that these methods retain the good long time behavior of the symplectic methods. We consider here an explicit adaptive step size strategy based on the time repara- metrization of the original Hamiltonian (6). Huang and Leimkuhler considered in [7] a modification of the Hamilton equations (6) of the following form with a re-parametrization s(p, q. E), where $ = &. In general, reparametri- zation will destroy the Hamiltonian structure, but the time-reversibility of the flow is preserved. The modified Hamiltonian system (13) can be solved by the second order Lobatto IIIa-b pair partitioned Runge-Kutta formula, which can be seen as a generalization of the Verlet method. Because of the coupling structure of the modified equation, the Verlet method becomes implicit. It was shown in [7] that for the case where the reparametrization is either based on p or q variables, the Verlet method becomes explicit. The generalized Verlet scheme for q-dependent reparametrizations s (q, E) is. given by (see [6]): where on+1/2 is defined recursively by: This scheme is explicit and symmetric, but it is not symplectic. As time reparame- trization function we have used the following q-dependent arclength parametrizations: In [6] Hairer introduced a new variable step size method, which preserves the symplectic structure. The new strategy is very cheap, since the solution of one scalar nonlinear equation is demanded in each step. The variable step size strategies discussed here are in development. They are tested mostly for Hamiltonian systems with few degrees of freedom, like Kepler's problem ([5], [6], [7], [13]) or the double pendulum and the three body problem [7]. Their effectiveness for large Hamiltonian system occurring in molecular dynamics have to be investigated in future. MOLECULAR DYNAMICS MODEL 4 BENCHMARK MODEL FROM MOLECULAR DYNAMICS As a benchmark we introduce an ideal and linear molecule of three atoms, see Fig. 2. Fig. 2. Molecule model. Linear model means that only motions in x-direction relatively to an equilibrium point are considered. The yi , z; coordinates i = 1,2, 3 are fixed. The difference xi - xi+], i = 1,2 stand for the deflections, the abbreviations AX: denote the equilibrium distances. The classical Hamilton function is represented by an usual kinetic energy term and by a potential enery term: The potential Vi,i+l represents the energy of oscillations between the neighbouring atoms. For further studies it is helpful to transform to new coordinates. n1, I Abbreviations: M = E:=, mi ; pi = m,+;+, , i= 1,2 1 " x3 -+ ql = - mixi "coordinate for the center of gravity" M. [=I For the transformed momenta pi it is assumed, that the motion of the center of gravity p3 is neglected, i.e. p3 = 0. The remaining momenta pl, pz are defined via the formula 290 M. HANKEL ET AL. This results in the transformed Hamiltonian For the potential V(q) the Morse potential is chosen. The constants Di stand for the dissociation energy and bi are experimentally obtained parameters of the Morse potential. The Hamiltonian equations of motion are given by Another way of dealing with the model described above allows the normal mode approach, where the Morse potentials are expanded in Taylor series and coordinate depending coupling between both modes is added. This model represents a strongly nonharmonic vibrational system with nonlinear cou- pling, and was studied in [16] in both the classical and the quantum mechanical senses. Remark: The choice of the Morse potential offers strong connections to the theory of Toda flows and to classical mechanics. The presented simplified model from MD can be reinterpreted as a mass chain. The mass points are coupled by springs with a nonlinear exponential characteristic. Moreover there are bridges to the standard eigenvalue problem in Numerical Linear Algebra. Using the Flaschka transformation (which destroys symplecticness) the resulting differential equations stand for a continuous version of the QR-algorithm. For further details see [lo]. In less simplified MD models these nice properties and connections to different areas are lost. MOLECULAR DYNAMICS MODEL 5 TEST RESULTS In a first test series the constant step size behavior of four integrators were studied for the two systems (20) and (21). The chosen integrators are (i) Verlet method (9) of order two, which is marked by + in the Figs. 3 to 6. (ii) Ruths method (1 1) of order three, which is marked by in the Figs. 3 to 6. (iii) Composite Verlet method (10) of order four, which is marked by 0 in the Figs. 3 to 6. (iv) For comparison the nonsymplectic classical Runge-Kutta method of order four is given, which is marked by x in the Figs. 3 to 6. The initial values and parameters of the decoupled system (20) and the coupled system (21) are listed in Table 1. Table 1. Parameters and initial values for (20) and (21). Decoupled (20) Coupled (2 1 ) pi(0) = -5.81602698201787 PI@) = 4 p2 (0) = - 1.460798931644 p2(0) = 2.593642227 ql(0) = -0.586227950237599 ql(0) = 0 q2(0) = 2.20649342825819 qz(0) = 1 bl = 0.531 dl = 5.93 b2 = 0.462 d2 = 10.01 same ml = 5.0079 m2 = 12.01 1 mj = 12.011 The relative energy error FE in the figures is measured at discrete points Eo stands for the initial energy and Ei denote the computed energy at time ti, i = 1, . . . , N. N was chosen as 1000. Figures 3 and 4 are related to the decoupled problem (20). In Fig. 3 the energy error is plotted against the constant step sizes The time integration starts with to = 0 and stops at tend = 5000. The straight lines correspond to our expectations. The composite Verlet performs best and the slopes of 292 M. HANKEL ET AL. Verlet and Ruth are two respectively three for large errors > lo-'. For smaller errors (20). The classical there seem to be some round-off in the exponential function in Runge-Kutta method works irregularly as expected for tolerances > In general the symplectic methods are superior. Figure 4 shows the similar picture for the energy error vs. CPU-time. Fixing computation time, the symplectic methods are superior. Fig. 3. Energy error vs. step size for (20). Fig. 3. Energy error vs. CPU time for (20). Figure 5 and 6 are related to the coupled problem (21). The step sizes are The slopes of the symplectic integrators are very regular and outperform the classi- cal Runge-Kutta. There are no round-off problems, since (2 1) only includes polynomial terms. In a second test series the adaptive step size strategies are studied. In general one would expect the adaptive strategies only superior in the presence of strong local forces, i.e. if collision models for molecules or nearly singular potential functions V(q) are incorporated. Since the Morse potential (19) describes bounded states, we do not expect significant differences. This can be seen in Table 2. The arclength parametrization (16) is comparable with constant step sizes. In this example (16) is more efficient than (1 5). To have some impressions on the solution structure, we have added the phase plots for ql - pl and 92 - p2 in Figs. 7 to 10. In case of one degree of freedom, which holds for the harmonic oscillator, the symplecticness can be easily visualized. A symplectic integrator will compute the classical circle, whereas a nonsymplectic method produces a cloud of points, see the very nice pictures in [3] and [Ill. For two or more degrees of freedom the figures are less clear. Therefore we choose the phase plots for the generalized momenta p, and generalized coordinates q,, i = 1,2. Symplecticness requires that the sum of the volumes, which is defined by the points, is constant for all MOLECULAR DYNAMICS MODEL 293 CPU 1ms Fig. 5. Energy emor vs. step size for (21). Fig. 6. Energy error vs. CPU time for (21). times t. Figures 7 and 8, respectively 9 and 10 show the time development. The initial values define closed orbits, see Figs. 7 and 9 with volumes Vl + V2 = const. where VI := VI (PI, ql; t = 0) , Fig. 7, V2 := V~(p2,q~; t = 0) , Fig. 9. The time history shows at t = T a snapshot with V3 := V3(pl, ql; t = T) , Fig. 8, V4 := V4(p2, 92; t = T) , Fig. 10. Due to the symplectic behavior it must hold: which can be easily seen in the figures. A visualization of the time history would show that Vl is gettting smaller, whereas V2 is still growing. This property is only reproduced by symplectic integration methods. Table 2. Adaptive step size strategies. h Number h,,, hmax max CPU t of steps IH - Hol constant step size .O1 500.000 ,316-1 11.30 arclength ,009 2.102.984 ,103-2 ,181-1 ,304-1 58.27 (15) .01 1.924.252 .I 16-2 ,219-1 ,381-1 53.30 arclength .05 697.372 ,571-3 ,720-2 ,675-2 25.66 .I441 (16) I 348.667 ,114-2 .271-1 12.83 294 M. HANKEL ET AL. Fig. 7. VI-phase plot (41, PI) at r = 0. Fig. 8. V3-phase plot (41, PI) at time T. Fig. 9. Vz-phase plot (42, p2) at t = 0. Fig. 10. V4-phase plot (92, p2) at time T. 6 CONCLUSIONS Development of symplectic integration methods for molecular dynamics and rigid body simulations is an actual research field. Although explicit schemes, such as the Verlet method and its higher order compositions are simple to formulate and propagate fast, they impose a severe constraint on the maximum step size. We expect especially for the design of variable step size strategies, which preserve the symplectic and reversible MOLECULAR DYNAMICS MODEL 295 structure of the system further progress and more numerical experimentation is needed in this area. In the low dimensional examples (20), (21) composite Verlet performs best. More complex molecules use the main computational time for the computation of force terms. This requires special solution techniques like multiple time stepping. ACKNOWLEDGEMENTS B. Karasozen is grateful to the DAAD (Deutscher Akademischer Austauschdienst) and SFB 298 for the support during the time this paper was written. REFERENCES 1. Brickmann, J., Kast, S.M., Vollhardt, H., Reiling S., Trends in Molecular Dynamics Simulation Technique, in Frontiers in Chemical Dynamics, E. Yurtsever (ed.), Kluwer Academic Press, 1995, pp. 217-253. 2. Calvo, M.P., Hairer, E.: Accurate Long-Term Integration of Dynamical Systems, Applied Numerical Analysis, 18(1995), pp. 95-105. 3. Hairer, E., Nprrsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems, Springer, 1993. 4. Hairer, E.: Backward analysis of numerical integrators and symplectic methods, Annals of Numerical Mathematics, 1 (1994), pp. 107-1 32. 5. Hairer, E., Stoffer, D.: Reversible long-term integration with variable step sizes, to appear in SIAM J. Sci. Comput. 1997. 6. Hairer, E.: Variable time step integration with symplectic methods, Universite' de Gdneve, Dkpt. de Mathhatiques, Technical Report, 1996. 7. Huang, W., Leimkuhler, B.: The adaptive Verlet method, to appear in SIAM J. Sci. Comput. -. 1997. . . . 8. McLachlan, R.1: On the numerical integration of ordinary differential equations by sym- metric composition methods, SIAM J. Sci. Comput., 16(1995), pp. 151-168. 9. Nettesheim, P., Bomemann F.A., Schmidt, B., Schiitte, C.: An explicit and symplectic inte- grator for quantum-classical molecular dynamics, Chem. Phys. Lett., 256(1996), pp. 581-588. 10. Nanda, T.: Dtfferential equations and the QR algorithm, SIAM J. Numer. Anal., 12(1985), pp. 310-321. 11. Sanz-Sema, J.M., Calvo, M.P.: Numerical Hamiltonian Problems, Chapman & Hall, 1994. Discrete time-reversible propagation scheme for mixed 12. Schmitt, U., Brickmann, J.: quantum-classical dynamics, Chemical Physics, 208(1996), pp. 45-56. 13. Stoffer, D.: Variable Steps for Reversible Integration Methods, Computing, 55(1995), pp. 1-22. 14. Suzuki, M.: General theory of higher-order decompositions of exponential operators and symplectic integrators, Physics Letters A, 165(1992), pp. 387-395. 15. Yoshida, H.: Construction of higher order symplectic integrators, Physics Letters, 159(1990), pp. 262-268. 296 M. HANKEL ET AL. 16. Yurtsever, E., Brickmann J.: Does Quantum Mechanics Select Out Regularity and Lo- cal Mode Behaviour in Nonlinearly Coupled Vibrational Systems?, Berichte der Bunsen-Gesellschaft fiir physikalische Chemie, 94(1990), pp. 804-827. 17. Warshel, A.: Computer Modeling of Chemical Reactions in Proteins and Solitons, Wiley, New York, 199 1. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical Modelling of Systems Taylor & Francis

A molecular dynamics model for symplectic intergrators

Download PDF
15 pages

Loading next page...
 
/lp/taylor-francis/a-molecular-dynamics-model-for-symplectic-intergrators-Oioi8FmPYh

References (16)

Publisher
Taylor & Francis
Copyright
Copyright Taylor & Francis Group, LLC
ISSN
1381-2424
DOI
10.1080/13873959708837062
Publisher site
See Article on Publisher Site

Abstract

Mathematical Modelling of Systems 138 1-2424197IO304-282$12.00 @wets & Zeitlinger 1997, Vol. 3, NO. 4, pp. 282-296 A Molecular Dynamics Model for Symplectic Integrators M. HANKEL*, B. KARASOZEN~, P. RENTROP* AND U. SCHMI?T$ ABSTRACT Mechanical systems in the very large scale like in celestial mechanics or in the very small scale like in the molecular dynamics can be modelled without dissipation. The resulting Hamiltonian systems possess conservation properties, which are characterized with the term symplecticness. Numerical integration schemes should preserve the symplecticness. Different methods are introduced and their performance is studied for constant and variable step size. As test examples two systems from molecular dynamics are introduced. Keywords: Hamiltonian systems, molecular models, symplectic integrators. 1 INTRODUCTION TO HAMILTONIAN SYSTEMS Hamiltonian systems play a key role in applications arising from the classical, statistical and quantum dynamics, from optics and plasma physics. In principle all phenomena where dissipation is absent or can be ignored may be modelled by Hamiltonian systems. The symplectic methods developed for the Hamiltonian systems in the recent years have become already a well established part of numerical analysis (see the monographs [3], [I 11 ). But they are far from the notion of a code (package) in the sense of standard numerical integrators with automatic error control and step size change. The numerical examples given in the literature for symplectic integrators consist mostly of systems with few degrees of freedom, like the harmonic oscillator, Kepler's problem, etc. The use of symplectic integrators for large systems of equations in molecular dynamics and for systems arising by space discretization of partial differential equations is in the beginning stage now. 'Fachbereich Mathematik, Technische Hochschule Darmstadt, Germany. l~e~axtment of Mathematics, Middle East Technical University, Ankara, Turkey. *~h~sikalische Chemie I, Technische Hochschule Darmstadt, Germany. MOLECULAR DYNAMICS MODEL 283 The Hamiltonian equations are a special kind of ordinary differential equations. The initial value problem is in Hamiltonian form if its dimension D is even, D = 2d, where d denotes the number of degree of freedom, and the vector field f (x) is of the form with the 2d x 2d skew symmetric matrix J VH(x) is the gradient of the real valued Hamilton function H(x). With x = [PI, . . . , pd, 41, . . . , qdlT, the equations (11, (2) are The flow 4, of a Hamiltonian system is a symplectic transformation for each t. In general, a transformation Q in lR2d is symplectic if its Jacobian Q1(x) satisfies elT (x) JQ1(x) = J. An equivalent definition in terms of differential forms exists too. The flow of the Hamiltonian system (3) is symplectic, if the differential two-form is preserved. Differential forms provide a geometric interpretation of the symplectic- ness in terms of conservation of areas. Another preserved quantity is the Hamiltonian, i.e if p(t) and q (t) are solutions of (3) then H (p(t), q (t)) is constant for all t. In the following we will consider special types of Hamiltonian systems occurring in molecular dynamics. Molecular dynamics (MD) simulations have been proven to be a powerful tool to obtain detailed insight into the behavior of dynamical processes in isolated molecules, clusters, biomolecules, liquids and surfaces [I], [17]. The dynamics of a molecular system can be either formulated in the classical, quantum or mixed quantum-classical framework. In a classical molecular dynamics or mechanics formulation the Hamiltonian is given by 284 M. HANKEL ET AL. where M contains the masses mi of the molecules and is an invertible symmetric matrix. The first part of the Hamiltonian denotes the kinetic energy T(p), the second part the potential energy V(q). Usually the variables q; denote the generalized coordinates and pi the generalized momenta. The resulting Hamiltonian system has the form The quantum dynamics problems can be formulated as the discretized Schrodinger equation (h = 1) in matrix representation with the Hamilton operator Replacing the complex vector c(t) by a linear combination of formal momentum and position coordinates and separating real and imaginary parts, the classical Hamiltonian formulation with HQ is then obtained where Hij are the elements of the Hamilton operator % in (7) and N stands for the number of pseudo particles. By adding the classical Hamilton function Hc(p, q) governing the dynamics of the classical part of the system to the Hamilton function HQ(p, q) describing the quantum degrees of freedom the mixed quantum-classical system can be interpreted as a set of purely classical Hamilton equations of motion. The matrix elements Hij contain now the coupling between the classical and quantum modes. An important property of the Hamiltonian systems (6) is the reversibility. The system (1) is called p-reversible, if there exists an involution p, that is a linear mapping in lRW such that p2 = I with the property for all x. f (px) = -p f (x), This implies that the flow $,(x) of the system (1) is a p-reversible mapping for each t MOLECULAR DYNAMICS MODEL 285 Since 4;' = 4-, the system is also time-reversible, shortly reversible. The Hamilto- nian equations in (6) are reversible with respect to In MD simulations the number of particles necessary to obtain statistically reliable results can easily exceed N lo5, the 6N coupled first order differential equations (3) have to be integrated numerically with an efficient and robust algorithm. Therefore the need for improved algorithms has received greater attention. Recent studies have shown that symplectic techniques can be used to derive stable and efficient integrators useful in the classical and mixed quantum-classical MD simulations, and furthermore a time-scale separation can be readily overcome by introducing different time step lengths [9], [12]. Figure 1 shows a snapshot from a mixed quantum-classical MD simulation of a charged Ar,19] sir n luster performed with a multiple time step sym- plectic integrator. The charge distribution is visualized as gray scales on the atomic sites ranging from black (low charge) to white (high charge). Fig. 1. Aqls1 luster. 2 NUMERICAL HAMILTONIAN SOLVERS There exists mainly three classes of symplectic integrators; methods based on gen- erating functions, methods using compositions and symplectic Runge-Kutta methods (see the monographs [l 11, [3], Section 11.16). The generating function methods seem to be too costly. The most effective symplectic integrators are belonging to the other 286 M. HANKEL ET AL. two classes; some of them specially designed for the integration of the Hamiltonian systems (6) will be given here. Symplectic integrators using composition The separable Hamiltonian (5) provides the most common application of the com- position approach. For many systems different parts of the vector field f (x) often correspond to physically different contributions, so that a natural splitting in form of H = H1 + H2, with HI = T(p) and H2 = V(q) is available. Many Hamilto- nian systems (6) are integrable in closed form after splitting. The most widespread second order composition method is the so called VerletJStormerlleap-frog method. Application of the Verlet method to (6) gives The Verlet method is symplectic and reversible. Higher order composition methods are obtained by using the basic method qh, here the Verlet method. For example the fourth order symplectic, time-reversible composition method is given by with the coefficients cl = &. e2 = -A 2-fi. Several types of composlhon methods were derived and their properties are inves- tigated using the Baker-Campbell-Hausdorff (BCH) formula, (see [8], [14], [15]). Symplectic Runge-Kutta methods For separable Hamiltonian problems (6) there exist explicit partitioned Runge-Kutta (PRK) methods, characterized by two Butcher tableaus MOLECULAR DYNAMICS MODEL 287 where the first one is used for the p variable and the other for the q variables. The PRK method applied to the Hamiltonian system (6) results in The symplecticness condition for the irreducible PRK methods is given by ( see [l 11, Section 6.3) Remark: Hamiltonian flows exactly preserve the symplectic structure and the Hamiltonian. In general for a symplectic integrator it is impossible to preserve the Hamiltonian or the energy exactly. In practice the energy error remains constant in long terms for symplectic integrators, whereas it increases or decreases for non-symplectic integrators. This can be explained by the backward error analysis, which states that the numerical solution obtained by the symplectic integrator can be formally inter- preted as the exact solution of a perturbed nearby Hamiltonian system [4]. Symplectic integrators also have more favourable properties concerning the global error growth of the solutions when the integration is performed over very long time intervals . For integrable systems and for problems with periodic solutions, the error growth is only linear in case of symplectic integrators, compared with the quadratic error growth in the general case [2]. 3 VARIABLE STEP SIZE STRATEGIES The favourable properties of the symplectic methods for Hamiltonian systems are lost in a standard variable step size implementation based on local error estimators; they are no longer superior to the non-symplectic algorithms. Therefore the development of variable time symplectic integrators was indispensable to be competitive with the standard methods. In the following we will describe three different non-standard variable step size strategies, which preserve the geometric structure of the Hamiltonian system (6) in long time integration. There exist several variable step size strategies for reversible systems using symmetric Runge-Kutta methods. It was shown in [2], [5], [13] that these methods retain the good long time behavior of the symplectic methods. We consider here an explicit adaptive step size strategy based on the time repara- metrization of the original Hamiltonian (6). Huang and Leimkuhler considered in [7] a modification of the Hamilton equations (6) of the following form with a re-parametrization s(p, q. E), where $ = &. In general, reparametri- zation will destroy the Hamiltonian structure, but the time-reversibility of the flow is preserved. The modified Hamiltonian system (13) can be solved by the second order Lobatto IIIa-b pair partitioned Runge-Kutta formula, which can be seen as a generalization of the Verlet method. Because of the coupling structure of the modified equation, the Verlet method becomes implicit. It was shown in [7] that for the case where the reparametrization is either based on p or q variables, the Verlet method becomes explicit. The generalized Verlet scheme for q-dependent reparametrizations s (q, E) is. given by (see [6]): where on+1/2 is defined recursively by: This scheme is explicit and symmetric, but it is not symplectic. As time reparame- trization function we have used the following q-dependent arclength parametrizations: In [6] Hairer introduced a new variable step size method, which preserves the symplectic structure. The new strategy is very cheap, since the solution of one scalar nonlinear equation is demanded in each step. The variable step size strategies discussed here are in development. They are tested mostly for Hamiltonian systems with few degrees of freedom, like Kepler's problem ([5], [6], [7], [13]) or the double pendulum and the three body problem [7]. Their effectiveness for large Hamiltonian system occurring in molecular dynamics have to be investigated in future. MOLECULAR DYNAMICS MODEL 4 BENCHMARK MODEL FROM MOLECULAR DYNAMICS As a benchmark we introduce an ideal and linear molecule of three atoms, see Fig. 2. Fig. 2. Molecule model. Linear model means that only motions in x-direction relatively to an equilibrium point are considered. The yi , z; coordinates i = 1,2, 3 are fixed. The difference xi - xi+], i = 1,2 stand for the deflections, the abbreviations AX: denote the equilibrium distances. The classical Hamilton function is represented by an usual kinetic energy term and by a potential enery term: The potential Vi,i+l represents the energy of oscillations between the neighbouring atoms. For further studies it is helpful to transform to new coordinates. n1, I Abbreviations: M = E:=, mi ; pi = m,+;+, , i= 1,2 1 " x3 -+ ql = - mixi "coordinate for the center of gravity" M. [=I For the transformed momenta pi it is assumed, that the motion of the center of gravity p3 is neglected, i.e. p3 = 0. The remaining momenta pl, pz are defined via the formula 290 M. HANKEL ET AL. This results in the transformed Hamiltonian For the potential V(q) the Morse potential is chosen. The constants Di stand for the dissociation energy and bi are experimentally obtained parameters of the Morse potential. The Hamiltonian equations of motion are given by Another way of dealing with the model described above allows the normal mode approach, where the Morse potentials are expanded in Taylor series and coordinate depending coupling between both modes is added. This model represents a strongly nonharmonic vibrational system with nonlinear cou- pling, and was studied in [16] in both the classical and the quantum mechanical senses. Remark: The choice of the Morse potential offers strong connections to the theory of Toda flows and to classical mechanics. The presented simplified model from MD can be reinterpreted as a mass chain. The mass points are coupled by springs with a nonlinear exponential characteristic. Moreover there are bridges to the standard eigenvalue problem in Numerical Linear Algebra. Using the Flaschka transformation (which destroys symplecticness) the resulting differential equations stand for a continuous version of the QR-algorithm. For further details see [lo]. In less simplified MD models these nice properties and connections to different areas are lost. MOLECULAR DYNAMICS MODEL 5 TEST RESULTS In a first test series the constant step size behavior of four integrators were studied for the two systems (20) and (21). The chosen integrators are (i) Verlet method (9) of order two, which is marked by + in the Figs. 3 to 6. (ii) Ruths method (1 1) of order three, which is marked by in the Figs. 3 to 6. (iii) Composite Verlet method (10) of order four, which is marked by 0 in the Figs. 3 to 6. (iv) For comparison the nonsymplectic classical Runge-Kutta method of order four is given, which is marked by x in the Figs. 3 to 6. The initial values and parameters of the decoupled system (20) and the coupled system (21) are listed in Table 1. Table 1. Parameters and initial values for (20) and (21). Decoupled (20) Coupled (2 1 ) pi(0) = -5.81602698201787 PI@) = 4 p2 (0) = - 1.460798931644 p2(0) = 2.593642227 ql(0) = -0.586227950237599 ql(0) = 0 q2(0) = 2.20649342825819 qz(0) = 1 bl = 0.531 dl = 5.93 b2 = 0.462 d2 = 10.01 same ml = 5.0079 m2 = 12.01 1 mj = 12.011 The relative energy error FE in the figures is measured at discrete points Eo stands for the initial energy and Ei denote the computed energy at time ti, i = 1, . . . , N. N was chosen as 1000. Figures 3 and 4 are related to the decoupled problem (20). In Fig. 3 the energy error is plotted against the constant step sizes The time integration starts with to = 0 and stops at tend = 5000. The straight lines correspond to our expectations. The composite Verlet performs best and the slopes of 292 M. HANKEL ET AL. Verlet and Ruth are two respectively three for large errors > lo-'. For smaller errors (20). The classical there seem to be some round-off in the exponential function in Runge-Kutta method works irregularly as expected for tolerances > In general the symplectic methods are superior. Figure 4 shows the similar picture for the energy error vs. CPU-time. Fixing computation time, the symplectic methods are superior. Fig. 3. Energy error vs. step size for (20). Fig. 3. Energy error vs. CPU time for (20). Figure 5 and 6 are related to the coupled problem (21). The step sizes are The slopes of the symplectic integrators are very regular and outperform the classi- cal Runge-Kutta. There are no round-off problems, since (2 1) only includes polynomial terms. In a second test series the adaptive step size strategies are studied. In general one would expect the adaptive strategies only superior in the presence of strong local forces, i.e. if collision models for molecules or nearly singular potential functions V(q) are incorporated. Since the Morse potential (19) describes bounded states, we do not expect significant differences. This can be seen in Table 2. The arclength parametrization (16) is comparable with constant step sizes. In this example (16) is more efficient than (1 5). To have some impressions on the solution structure, we have added the phase plots for ql - pl and 92 - p2 in Figs. 7 to 10. In case of one degree of freedom, which holds for the harmonic oscillator, the symplecticness can be easily visualized. A symplectic integrator will compute the classical circle, whereas a nonsymplectic method produces a cloud of points, see the very nice pictures in [3] and [Ill. For two or more degrees of freedom the figures are less clear. Therefore we choose the phase plots for the generalized momenta p, and generalized coordinates q,, i = 1,2. Symplecticness requires that the sum of the volumes, which is defined by the points, is constant for all MOLECULAR DYNAMICS MODEL 293 CPU 1ms Fig. 5. Energy emor vs. step size for (21). Fig. 6. Energy error vs. CPU time for (21). times t. Figures 7 and 8, respectively 9 and 10 show the time development. The initial values define closed orbits, see Figs. 7 and 9 with volumes Vl + V2 = const. where VI := VI (PI, ql; t = 0) , Fig. 7, V2 := V~(p2,q~; t = 0) , Fig. 9. The time history shows at t = T a snapshot with V3 := V3(pl, ql; t = T) , Fig. 8, V4 := V4(p2, 92; t = T) , Fig. 10. Due to the symplectic behavior it must hold: which can be easily seen in the figures. A visualization of the time history would show that Vl is gettting smaller, whereas V2 is still growing. This property is only reproduced by symplectic integration methods. Table 2. Adaptive step size strategies. h Number h,,, hmax max CPU t of steps IH - Hol constant step size .O1 500.000 ,316-1 11.30 arclength ,009 2.102.984 ,103-2 ,181-1 ,304-1 58.27 (15) .01 1.924.252 .I 16-2 ,219-1 ,381-1 53.30 arclength .05 697.372 ,571-3 ,720-2 ,675-2 25.66 .I441 (16) I 348.667 ,114-2 .271-1 12.83 294 M. HANKEL ET AL. Fig. 7. VI-phase plot (41, PI) at r = 0. Fig. 8. V3-phase plot (41, PI) at time T. Fig. 9. Vz-phase plot (42, p2) at t = 0. Fig. 10. V4-phase plot (92, p2) at time T. 6 CONCLUSIONS Development of symplectic integration methods for molecular dynamics and rigid body simulations is an actual research field. Although explicit schemes, such as the Verlet method and its higher order compositions are simple to formulate and propagate fast, they impose a severe constraint on the maximum step size. We expect especially for the design of variable step size strategies, which preserve the symplectic and reversible MOLECULAR DYNAMICS MODEL 295 structure of the system further progress and more numerical experimentation is needed in this area. In the low dimensional examples (20), (21) composite Verlet performs best. More complex molecules use the main computational time for the computation of force terms. This requires special solution techniques like multiple time stepping. ACKNOWLEDGEMENTS B. Karasozen is grateful to the DAAD (Deutscher Akademischer Austauschdienst) and SFB 298 for the support during the time this paper was written. REFERENCES 1. Brickmann, J., Kast, S.M., Vollhardt, H., Reiling S., Trends in Molecular Dynamics Simulation Technique, in Frontiers in Chemical Dynamics, E. Yurtsever (ed.), Kluwer Academic Press, 1995, pp. 217-253. 2. Calvo, M.P., Hairer, E.: Accurate Long-Term Integration of Dynamical Systems, Applied Numerical Analysis, 18(1995), pp. 95-105. 3. Hairer, E., Nprrsett, S.P., Wanner, G., Solving Ordinary Differential Equations I, Nonstiff Problems, Springer, 1993. 4. Hairer, E.: Backward analysis of numerical integrators and symplectic methods, Annals of Numerical Mathematics, 1 (1994), pp. 107-1 32. 5. Hairer, E., Stoffer, D.: Reversible long-term integration with variable step sizes, to appear in SIAM J. Sci. Comput. 1997. 6. Hairer, E.: Variable time step integration with symplectic methods, Universite' de Gdneve, Dkpt. de Mathhatiques, Technical Report, 1996. 7. Huang, W., Leimkuhler, B.: The adaptive Verlet method, to appear in SIAM J. Sci. Comput. -. 1997. . . . 8. McLachlan, R.1: On the numerical integration of ordinary differential equations by sym- metric composition methods, SIAM J. Sci. Comput., 16(1995), pp. 151-168. 9. Nettesheim, P., Bomemann F.A., Schmidt, B., Schiitte, C.: An explicit and symplectic inte- grator for quantum-classical molecular dynamics, Chem. Phys. Lett., 256(1996), pp. 581-588. 10. Nanda, T.: Dtfferential equations and the QR algorithm, SIAM J. Numer. Anal., 12(1985), pp. 310-321. 11. Sanz-Sema, J.M., Calvo, M.P.: Numerical Hamiltonian Problems, Chapman & Hall, 1994. Discrete time-reversible propagation scheme for mixed 12. Schmitt, U., Brickmann, J.: quantum-classical dynamics, Chemical Physics, 208(1996), pp. 45-56. 13. Stoffer, D.: Variable Steps for Reversible Integration Methods, Computing, 55(1995), pp. 1-22. 14. Suzuki, M.: General theory of higher-order decompositions of exponential operators and symplectic integrators, Physics Letters A, 165(1992), pp. 387-395. 15. Yoshida, H.: Construction of higher order symplectic integrators, Physics Letters, 159(1990), pp. 262-268. 296 M. HANKEL ET AL. 16. Yurtsever, E., Brickmann J.: Does Quantum Mechanics Select Out Regularity and Lo- cal Mode Behaviour in Nonlinearly Coupled Vibrational Systems?, Berichte der Bunsen-Gesellschaft fiir physikalische Chemie, 94(1990), pp. 804-827. 17. Warshel, A.: Computer Modeling of Chemical Reactions in Proteins and Solitons, Wiley, New York, 199 1.

Journal

Mathematical Modelling of Systems – Taylor & Francis

Published: Jan 1, 1997

Keywords: Hamiltonian systems; molecular models; symplectic integrators

There are no references for this article.