Modelling and Simulation in Engineering
Volume 2015 (2015), Article ID 904140, 8 pages
http://dx.doi.org/10.1155/2015/904140

Research Article

Continuous Finite Element Methods of Molecular Dynamics Simulations

Qiong Tang,¹ Luohua Liu,¹ and Yujun Zheng²

¹College of Science, Hunan University of Technology, Zhuzhou, Hunan 412007, China
²Department of Mathematics and Computational Science, Hunan University of Science and Engineering, Yongzhou, Hunan 425100, China

Received 13 February 2015; Accepted 26 July 2015

Academic Editor: Ahmed Rachid

Copyright © 2015 Qiong Tang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Molecular dynamics simulations are necessary to perform very long integration times. In this paper, we discuss continuous finite element methods for molecular dynamics simulation problems. Our numerical results about diatomic molecular system and triatomic molecules show that linear finite element and quadratic finite element methods can better preserve the motion characteristics of molecular dynamics, that is, properties of energy conservation and long-term stability. So finite element method is also a reliable method to simulate long-time classical trajectory of molecular systems.

1. Introduction

Molecular dynamics (MD) simulations can provide detailed information on molecular fluctuations and conformational changes and are used to investigate the thermodynamics and structure of chemical and biological molecules. Constructing efficient algorithms suitable for MD applications is an important task. Because of the very rapid oscillatory motions characteristic of MD, very small time steps of about a femtosecond () are typically used in MD. This means that, in order to get estimates for parameters of interest from simulation, it is necessary to perform very long integration times for MD simulations. The total energy , obtained by summing the kinetic and potential contributions, is a conserved quantity for the molecular system. Many numerical methods have been developed but the most effective for use in MD simulations should have superior long-term stability properties and energy conservation and permit a large integration time step [1–6]. It is important to construct discrete algorithms which preserve these basic properties.

Currently symplectic algorithm is the main choice in the long-time calculation to the MD. An important feature is that the molecular systems always remain Hamiltonian dynamical models. The Hamiltonian system has a symplectic structure with strong geometric properties of a dynamical system and conserves the total energy in which the phase-space points are allowed on the constant energy hypersurface satisfying . Ruth, Feng, Sanz-Serna, and Marsden et al. constructed the symplectic algorithm to solve the Hamiltonian system, including the 4-stage, 4th-order explicit symplectic scheme (4SS) to the separable Hamiltonian system and 2-stage, 4th-order implicit symplectic Runge-Kutta scheme (4SRK) [4, 7–13]. It is a difference method that preserves the structure of the system. However, most numerical methods cannot maintain the two properties: Symplectic and energy conservation simultaneously in general according to Ge-Marsde theorem [14]. Symplectic algorithm can preserve symplectic properties, but it can only obtain approximate energy conservation for nonlinear Hamiltonian system. However, the energy conservation is very important to the MD system. If not, do not stay near the energy surface very soon leading to unrealistic averages [4–6, 15].

Applying continuous finite element methods for Hamiltonian systems, there are two important advantages where finite element methods can preserve energy conservation and high accuracy approximation to symplectic structure. Zhong and Chen et al. [16–18] got good results. Based on the importance of energy conservation and long time to the MD system, in this paper, we first use continuous finite element methods to solve MD simulation problems. In numerical experiments of diatomic molecular and triatomic molecules, the linear element and quadratic element methods can better preserve the motion characteristics of molecular system: conservation of energy and long-term stability properties, that reached the microscopic reaction need to be considered time. It is a reliable method for long-time classical trajectory simulations of MD.

2. Continuous Finite Element Methods for Hamiltonian Systems

Consider the first-order ordinary differential equation with initial value in the interval :

Let be a partition of , with , . Consider , , . The constant is independent of and . Define the th degree continuous finite element space:

where is a set of th degree polynomials. In cell , th degree polynomial has parameters. As it is an initial value problem, and we already knew the initial value in , it has only degrees of freedom. Defining th degree continuous finite element satisfies [19, 20]

That is, in cell , it is orthogonal to arbitrary . Taking , then its derivate . In practical computation, take , .

Lemma 1 (see [20]). The th degree continuous finite element of ordinary differential equation has superconvergence in cell node :

We take finite dimension autonomous Hamiltonian canonical systems as

where , , and matrix transpose is defined by “.” In applications to MD, , , and represent the particles positions, the canonical momentum, and the total energy of the systems, respectively.

Let ; (5) can be written as

where , is the order identity matrix, and

Define (6)’s th degree continuous finite element of ; it satisfies orthogonal relation:

Taking , we obtain

By subtracting (9) from (8),

Hence, in each node , we prove that

Then we can prove the following.

Lemma 2 (see [16]). Applying arbitrary degree continuous finite element to solve Hamilton equation, it preserves energy conservation; that is, .
Using continuous finite element methods to solve nonautonomous Hamiltonian systems, from (7), we can obtain

By subtracting (13) from (12),

It is a nonautonomous Hamiltonian system; in general for every cell . Hence, in node ,

Theorem 3. The continuous finite element methods cannot preserve the energy to nonautonomous Hamiltonian systems.

Theorem 3 verifies that the nonautonomous Hamiltonian systems do not have the energy conservation, which coincides with the actual case.

Definition 4 (see [21]). A transformation of is called canonical, or symplectic, if it is a local diffeomorphism whose Jacobian is everywhere symplectic; that is,

Utilizing Legendre orthogonal polynomial and Lemmas 1 and 2, we can prove the following.

Lemma 5 (see [17]). The linear finite element for nonlinear Hamiltonian systems is an approximately symplectic algorithm which has the accuracy of third order to Hamiltonian systems symplectic structure; that is, It also possesses second-order accuracy and preserves the energy.

Lemma 6 (see [17]). The quadratic finite element for nonlinear Hamiltonian systems is an approximately symplectic algorithm which has the accuracy of fifth order to Hamiltonian systems symplectic structure; that is, It also possesses fourth-order accuracy and preserves the energy.

The specific calculation format of the linear continuous finite element (1FE) , , is as follows:

The specific calculation format of the quadratic continuous finite element (2FE) , , , is as follows:

We utilize high accuracy numerical integration such as at least points of the Gaussian quadrature formula to the th degree continuous finite element at the right of (17) or (18). From above equation set, we can solve on cell , and then we can get the value , .

3. Numerical Experiments

3.1. Diatomic Molecular System

Consider the classical motion of diatomic molecule in the electronic potential, where the mass of atom is and that of atom is . The coordinates of and are and , respectively. is the canonical coordinate; is the canonical mass. Then and are the momentum and canonical kinetics, respectively. Morse potentials have been used to describe vibrational levels of diatomic molecules. We take the Morse potential as [6]. In each case, the parameters , , and of the corresponding Morse potential are specified. The total energy of the diatomic molecule system is

And the Hamiltonian equation is

Since it is a separated Hamiltonian system, symplectic scheme can be used to solve this system [6]. We consider 2-stage, 4th-order implicit symplectic Runge-Kutta scheme (4SRK):

The 4-stage, 4th-order explicit symplectic difference scheme (4SS) is

where   ,  ,  ,  ,  and

Take the potential parameters , , , , , and to the CO molecule [5]. Initial condition , , stepsize , computing step, and integral interval . Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error; two atoms oscillate versus time and classical trajectories of the CO diatomic systems as follows (Figures 1–6).

Figure 1: 2FE, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space .

Figure 2: 2FE, computing step, the evolution of the final 100000 points .

Figure 3: 4SRK, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space .

Figure 4: 4SRK, computing step, the evolution of the final 100000 points .

Figure 5: 4SS, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space .

Figure 6: 4SS, computing step, the evolution of the final 100000 points .

Take the potential parameters , , , , , and to the molecule [5]. Initial condition , , stepsize , computing step, and integral interval . Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error and classical trajectories of the diatomic systems as follows (Figures 7–9).

Figure 7: 2FE, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of molecule in the phase space .

Figure 8: 4SRK, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of molecule in the phase space .

Figure 9: 4SS, computing step, the final 100000 points’ corresponding energy error curve and phase trajectories of molecule in the phase space .

Figures 1, 3, 5, 7, 8, and 9 show the evolution of the energy and the phase trajectories of CO, molecule versus time. Figures 2, 4, and 6 indicate that two atoms oscillate in the evolution process, respectively, utilizing 2FE, 4SRK, and 4SS. It is observed from Figures 1–9 that the numerical results computed by 2FE are in agreement with the theoretical analysis, the phase trajectories in phase space are steady, and two atoms oscillate periodically in the evolution process (Figures 1, 2, and 7) which indicates that 2FE can keep approximate high accuracy of symplectic structure just as the symplectic difference method (Figures 3, 4, 5, 6, 8, and 9) which preserves the structure of phase space in a long time . The energy error computed by 2FE is only (Figures 1 and 7) which indicates conservation of energy over long time intervals under a large stepsize , but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to (4SRK, Figures 3 and 8) and (4SS, Figures 5 and 9), which is only approximately conservative.

3.2. Triatomic Molecules

Basing Banerjee and Adams constructed a transformation method; we consider the motion of triatomic molecule such as water () within the symmetry; take the Cartesian coordinate system , with origin at the center of mass , and the -axis is the axis; the coordinates of the two atoms and atom are , , and . The generalized coordinates , , and the generalized momenta , Hamiltonian function of type molecule is , where is the kinetic energy and is the potential energy:

The canonical differential equations are

Initial condition: , , , , and Taking stepsize , computing step, integral interval .

The second-order symplectic difference scheme (2SS) is

We use the linear element method (1FE), general Leapfrog method [1], 2SS, and 4SS to compute the classical trajectories of type molecule in phase spaces and energy error as follows (Figures 10–13).

Figure 10: 1FE, computing step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.

Figure 11: General Leapfrog method, computing step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.

Figure 12: 2SS, computing step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.

Figure 13: 4SS, computing step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.

It is observed from Figures 10–13 that the numerical results computed by the linear element methods are in agreement with the theoretical analysis; atom and two atoms in type molecule vibrate quasi-periodically and the phase spaces are not squeezed together (Figure 10) which indicates that 1FE can for a long time preserve the characteristic of the oscillatory motions when , just as the symplectic difference method (Figures 11 and 12). The energy error computed by the linear element methods is only when , it reached the microscopic reaction needed to be considered time , but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to (Leapfrog method, Figure 11), (2SS, Figure 12), and (4SS, Figure 13), which is only approximately conservative.

4. Conclusion

The above computations indicate that the finite element methods can preserve less energy drift and are extremely stable over long runs to MD, which leads to an excellent long-time behaviour. These just meet the requirements of MD simulations which have superior long-term stability properties, conservation of energy, and a large integration time step.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China (no. 11101136) and aid program for Science and Technology Innovative Research Team in Higher Educational Institutions of Hunan. The authors thank the anonymous reviewer for their useful comments which were helpful to the presentation of this paper.

References

S. D. Bond, B. J. Leimkuhler, and B. J. Laird, “The nosé-poincaré method for constant temperature molecular dynamics,” Journal of Computational Physics, vol. 151, no. 1, pp. 114–134, 1999.

S. D. Bond and B. J. Leimkuhler, “Molecular dynamics and the accuracy of numerically computed averages,” Acta Numerica, vol. 16, pp. 1–65, 2007.

J. B. Sturgeon and B. B. Laird, “Symplectic algorithm for constant-pressure molecular dynamics using a Nosé-Poincaré thermostat,” Journal of Chemical Physics, vol. 112, no. 8, pp. 3474–3482, 2000.

B. Leimkuhler and S. Reich, Simulating Hamiltonian Dynamics, Cambridge University Press, Cambridge, UK, 2004.

E. Ley-Koo, S. Mateos-Cortés, and G. Villa-Torres, “Vibrational-rotational levels and Franck-Condon factors of diatomic molecules via Morse potentials in a box,” International Journal of Quantum Chemistry, vol. 58, no. 1, pp. 23–28, 1996.

X.-S. Liu, Y.-Y. Qi, J.-F. He, and P.-Z. Ding, “Recent progress in symplectic algorithms for use in quantum systems,” Communications in Computational Physics, vol. 2, no. 1, pp. 1–53, 2007.

K. Feng and M. Z. Qin, Symplectic Geometry Algorithm for Hamiltonian systems, Zhejiang Press of Science and Technology, Hangzhou, China, 2004, (Chinese).

O. Gonzalez and J. C. Simo, “On the stability of symplectic and energy-momentum algorithms for non-linear Hamiltonian systems with symmetry,” Computer Methods in Applied Mechanics and Engineering, vol. 134, no. 3-4, pp. 197–222, 1996.

C. Kane, J. E. Marsden, and M. Ortiz, “Symplectic-energy-momentum preserving variational integrators,” Journal of Mathematical Physics, vol. 40, no. 7, pp. 3353–3371, 1999.

M. Z. Qin, D. L. Wang, and M. Q. Zhang, “Explicit symplectic difference schemes for separable Hamiltonian system,” Journal of Computational Mathematics, vol. 9, no. 3, pp. 211–221, 1991.

S. Reich, “Multi-symplectic Runge-Kutta collocation methods for Hamiltonian wave equations,” Journal of Computational Physics, vol. 157, no. 2, pp. 473–499, 2000.

J. M. Sanz-Serna and M. P. Calvo, Numerical Hamiltonian Problems, St Edmundsbury Press, Bury St Edmunds, UK, 1994.

A. M. Stuart and A. R. Humphries, Dynamical Systems and Numerical Analysis, Cambridge University Press, 1998.

G. Zhong and J. E. Marsden, “Lie-poisson hamilton-jacobi theory and lie-poisson integrators,” Physics Letters A, vol. 133, no. 3, pp. 134–139, 1988.

E. Hairer, C. Lubich, and G. Wanner, Geometric Numerical Integration: Structure—Preserving Algorithms for Ordinary Differential Equations, vol. 31 of Springer Series in Computational Mathematics, Springer, Berlin, Germany, 2nd edition, 2006.

Q. Tang, C.-m. Chen, and L.-h. Liu, “Energy conservation and symplectic properties of continuous finite element methods for Hamiltonian systems,” Applied Mathematics and Computation, vol. 181, no. 2, pp. 1357–1368, 2006.

Q. Tang, C. M. Chen, and L. H. Liu, “Finite element methods for Hamiltonian systems,” Mathematica Numerica Sinica, vol. 31, no. 4, pp. 393–406, 2009.

W. X. Zhong and Z. Yao, “Time domain FEM and symplectic conservation,” Journal of Mechanical Strength, vol. 27, no. 2, pp. 178–183, 2005.

C. M. Chen and Y. Q. Huang, High Accuracy Theory of Finite Element, Hunan Press of Science and Technology, Changsha, China, 1995 (Chinese).

C. M. Chen, Finite Element Superconvergence Construction Theory, Hunan Press of Science and Technology, Changsha, China, 2001 (Chinese).

K. Feng, Collected Works of Feng Kang, National Defence Industry Press, Beijing, China, 1995.