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 aswhere , , 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 aswhere , is the order identity matrix, and
Define (6)’s th degree continuous finite element of ; it satisfies orthogonal relation:Taking , we obtainBy subtracting (9) from (8),Hence, in each node , we prove thatThen 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 obtainBy 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 isAnd the Hamiltonian equation isSince 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) iswhere , , , , 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).
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).
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 areInitial condition: , , , , and Taking stepsize , computing step, integral interval .
The second-order symplectic difference scheme (2SS) isWe 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).
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.