# Continuous Finite Element Methods of Molecular Dynamics Simulations

Authors: Qiong Tang; Luohua Liu; Yujun Zheng
Journal: Modelling and Simulation in Engineering
 (2015)

---

## 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 aboutAB diatomic molecular system and A2B 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.


---

## Body

## 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 (10-15) 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 H, 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 energyH(q,p) in which the phase-space points (q,p) are allowed on the constant energy hypersurface satisfying H(q,p)=E. 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, (p,q) 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 AB diatomic molecular and A2B 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 intervalJ=[0,T]: (1)dudt=u′=ft,u,u0=u0.Let Jh:t0=0<t1<t2<⋯<tN=T be a partition of J, with Ij=(tj-1,tj), hj=tj-tj-1. Consider h=maxhj, hj≥Ch, 1≤j≤N. The constant C is independent of j and h. Define the mth degree continuous finite element space: (2)Sh=w∣w∈CJ,w∣Ij∈Pm,where Pm is a set of mth degree polynomials. In cell Ij, mth degree polynomial has m+1 parameters. As it is an initial value problem, and we already knew the initial value U(tj-1) in Ij, it has only m degrees of freedom. Defining mth degree continuous finite element U∈Sh satisfies [19, 20](3)∫IjU′-ft,Uvdt=0,v∈Pm-1,U0=u0.That is, in cell Ij, it is orthogonal to arbitrary Pm-1. Taking v∈Sh, then its derivate v′∈Pm-1. In practical computation, take v=(t-tj-1)i, i=0,1,2,…,m-1.Lemma 1 (see [20]).
Themth degree continuous finite element of ordinary differential equation has superconvergence in cell node tj:(4)u-Utj=Oh2mum+1,∞.We take finite dimension autonomous HamiltonianH(p,q) canonical systems as(5)dpdt=-∂H∂q,dqdt=∂H∂p,where p=(p1,p2,…,pn)T, q=(q1,q2,…,qn)T, and matrix transpose is defined by “T.” In applications to MD, q, p, and H represent the n particles positions, the canonical momentum, and the total energy of the systems, respectively.Letz=[p,q]T; (5) can be written as(6)dzdt=z′=J-1Hz=fz,where J=0In-In0, In is the n order identity matrix, and Hz=HpHq.Define (6)’s mth degree continuous finite element Z=[P,Q]T of z; it satisfies orthogonal relation:(7)∫IjZ′+JHzv′dt=0,Z0=z0.Taking v=[Q,P]T, we obtain(8)∫IjP′+HqP,QQ′dt=0,P0=p0,(9)∫IjQ′-HpP,QP′dt=0,Q0=q0.By subtracting (9) from (8),(10)∫IjHpP′+HqQ′dt=∫IjddtHP,Qdt=0.Hence, in each node tj, we prove that(11)HPtj,Qtj=HPtj-1,Qtj-1=⋯=HP0,Q0=Hp0,q0.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,HPtj,Qtj=H(p0,q0).
Using continuous finite element methods to solve nonautonomous HamiltonianH(p,q,t) systems, from (7), we can obtain(12)∫IjP′+HqP,Q,tQ′dt=0,P0=p0,(13)∫IjQ′-HpP,Q,tP′dt=0,Q0=q0.By subtracting (13) from (12),(14)∫IjHpP,Q,tP′+HqP,Q,tQ′dt=∫IjddtHP,Q,t-HtP,Q,tdt=0.It is a nonautonomous Hamiltonian system; in general ∫IjHt(p,q,t)dt≠0 for every cell IJ. Hence, in node tj,(15)HPtj,Qtj,tj≠Hp0,q0,t0.Theorem 3.
The continuous finite element methods cannot preserve the energy to nonautonomous Hamiltonian systems.Theorem3 verifies that the nonautonomous Hamiltonian systems H(p,q,t) do not have the energy conservation, which coincides with the actual case.Definition 4 (see [21]).
A transformationz→z^ of R2n is called canonical, or symplectic, if it is a local diffeomorphism whose Jacobian ∂z^/∂z is everywhere symplectic; that is,(16)∂z^∂zTJ∂z^∂z=J.Utilizing Legendre orthogonal polynomial and Lemmas1 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,(∂Zj+1/∂Zj)TJ(∂Zj+1/∂Zj)=J+O(h3). 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,(∂Zj+1/∂Zj)TJ(∂Zj+1/∂Zj)=J+O(h5). It also possesses fourth-order accuracy and preserves the energy.The specific calculation format of the linear continuous finite element (1FE)Z=1-x/2Zj-1+1+x/2Zj, x∈[-1,1], is as follows:(17)Zj-Zj-1=hj2J∫-11HzZdx,j=0,1,…,N-1.The specific calculation format of the quadratic continuous finite element (2FE)Z=x2-x/2Zj-1+(1-x2)Zj-1/2+x2+x/2Zj, Zj-1/2=Z(tj-1+tj/2), x∈[-1,1], is as follows:(18)Zj-Zj-1=hj2J∫-11HzZdx,56Zj-23Zj-1/2-16Zj-1=hj4J∫-11HzZx+1dx.We utilize high accuracy numerical integration such as at least m+1 points of the Gaussian quadrature formula to the mth degree continuous finite element at the right of (17) or (18). From above equation set, we can solve Zj on cell Ij=(tj-1,tj), and then we can get the value Zj, j=1,…,N.


## 3. Numerical Experiments

### 3.1. Diatomic Molecular System

Consider the classical motion ofAB diatomic molecule in the electronic potential, where the mass of atom A is m1 and that of atom B is m2. The coordinates of A and B are x1 and x2, respectively. q=R=x1+x2 is the canonical coordinate; m=(m1+m2)/(m1m2) is the canonical mass. Then p=mqt and U=p2/(2m) 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 Vq=D(e-2α(q-q0)-e-α(q-q0)) [6]. In each case, the parameters D, α, and q0 of the corresponding Morse potential are specified. The total energy of the AB diatomic molecule system is(19)Hp,q=Up+Vq=12mp2+Vq.And the Hamiltonian equation is(20)dpdt=-∂H∂q=-∂Vq∂q,dqdt=∂H∂p=∂Up∂p.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):(21)Y1=zn+h14fY1+14-36fY2,Y2=zn+h14+36fY1+14fY2,zn+1=zn+h2fY1+fY2.The 4-stage, 4th-order explicit symplectic difference scheme (4SS) is(22)Hp,q=Tp+Vq,α=2-21/3-1,β=1-2α,p1=pn-c1h∂Vq∂qqn,q1=qn+d1h∂Tp∂pp1,p2=pn-c2h∂Vq∂qq1,q2=qn+d2h∂Tp∂pp2,p3=pn-c3h∂Vq∂qq2,q3=qn+d3h∂Tp∂pp3,pn+1=pn-c4h∂Vq∂qq3,qn+1=qn+d4h∂Tp∂ppn+1,where   c1=c4=α/2,  c2=c3=α+β/2,  d1=d2=α,  d3=β,  and d4=0.Take the potential parametersD=88575.905cm-1, α=2.3425473Å-1, m1=22036a.u., m2=29376a.u., m=m1∗m2/m1+m2, and q0=1.128323Å to the CO molecule [5]. Initial condition q(0)=q0, p(0)=2MD-0.01, stepsize h=3, computing K=108 step, and integral interval T=[0,3∗108]. Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error; two atoms oscillate versus time and (q,p) classical trajectories of the CO diatomic systems as follows (Figures 1–6).Figure 1
2FE, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 2
2FE, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Figure 3
4SRK, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 4
4SRK, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Figure 5
4SS, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 6
4SS, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Take the potential parametersD=51443.117cm-1, α=24.251821nm-1, m1=25704a.u., m2=m1, m=m1/2, and q0=0.12203nm to the N2 molecule [5]. Initial condition q(0)=q0, p(0)=2MD-0.01, stepsize h=3, computing K=108 step, and integral interval T=[0,3∗108]. Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error and (q,p) classical trajectories of the N2 diatomic systems as follows (Figures 7–9).Figure 7
2FE, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figure 8
4SRK, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figure 9
4SS, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figures1, 3, 5, 7, 8, and 9 show the evolution of the energy and the phase trajectories of CO, N2 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 (T=[0,3∗108]). The energy error computed by 2FE is only 10-15 (Figures 1 and 7) which indicates conservation of energy over long time intervals T=3∗108 under a large stepsize h=3, but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to 10-9 (4SRK, Figures 3 and 8) and 10-7 (4SS, Figures 5 and 9), which is only approximately conservative.


### 3.2.A2B Triatomic Molecules

Basing Banerjee and Adams constructed a transformation method; we consider the motion ofA2B triatomic molecule such as water (H2O) within the C2v symmetry; take the Cartesian coordinate system yOz, with origin at the center of mass O, and the z-axis is the C2 axis; the coordinates of the two A atoms and B atom are A(y1,z1), A(y2,z2), and B(y3,z3). The generalized coordinates q1=y1-y2, q2=z1-2z3+z2, and the generalized momenta p1=0.25dq1/dt, p2=0.5dq2/dt. Hamiltonian function of A2B type molecule is H(q,p)=K(p)+V(q), where K(p)=2p12+p22 is the kinetic energy and V(q) is the potential energy: (23)Vq=5π2D2-5D+6.5+4D-1+0.5π2q2-1.52+q2-1,D=q12+q221/2.The canonical differential equations are(24)dq1dt=∂H∂p1,dq2dt=∂H∂p2,dp1dt=-∂H∂q1,dp2dt=-∂H∂q2.Initial condition: q1(0)=3, q2(0)=3/2, p1(0)=0, p2(0)=0, and H(0)=50.1951. Taking stepsize h=0.01, computing K=108 step, integral interval T=[0,106].The second-order symplectic difference scheme (2SS) is(25)zk+1=zk+hJ-1Hzzk+1+zk2.We use the linear element method (1FE), general Leapfrog method [1], 2SS, and 4SS to compute the classical trajectories of A2B type molecule in phase spaces and energy error Hh-H as follows (Figures 10–13).Figure 10
1FE, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 11
General Leapfrog method, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 12
2SS, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 13
4SS, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.It is observed from Figures10–13 that the numerical results computed by the linear element methods are in agreement with the theoretical analysis; atom B and two A atoms in A2B 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 T=1000000, just as the symplectic difference method (Figures 11 and 12). The energy error computed by the linear element methods is only 10-10 when T=1000000, it reached the microscopic reaction needed to be considered time T=h∗K=106, but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to 10-1 (Leapfrog method, Figure 11), 10-2 (2SS, Figure 12), and 10-3 (4SS, Figure 13), which is only approximately conservative.


## 3.1. Diatomic Molecular System

Consider the classical motion ofAB diatomic molecule in the electronic potential, where the mass of atom A is m1 and that of atom B is m2. The coordinates of A and B are x1 and x2, respectively. q=R=x1+x2 is the canonical coordinate; m=(m1+m2)/(m1m2) is the canonical mass. Then p=mqt and U=p2/(2m) 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 Vq=D(e-2α(q-q0)-e-α(q-q0)) [6]. In each case, the parameters D, α, and q0 of the corresponding Morse potential are specified. The total energy of the AB diatomic molecule system is(19)Hp,q=Up+Vq=12mp2+Vq.And the Hamiltonian equation is(20)dpdt=-∂H∂q=-∂Vq∂q,dqdt=∂H∂p=∂Up∂p.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):(21)Y1=zn+h14fY1+14-36fY2,Y2=zn+h14+36fY1+14fY2,zn+1=zn+h2fY1+fY2.The 4-stage, 4th-order explicit symplectic difference scheme (4SS) is(22)Hp,q=Tp+Vq,α=2-21/3-1,β=1-2α,p1=pn-c1h∂Vq∂qqn,q1=qn+d1h∂Tp∂pp1,p2=pn-c2h∂Vq∂qq1,q2=qn+d2h∂Tp∂pp2,p3=pn-c3h∂Vq∂qq2,q3=qn+d3h∂Tp∂pp3,pn+1=pn-c4h∂Vq∂qq3,qn+1=qn+d4h∂Tp∂ppn+1,where   c1=c4=α/2,  c2=c3=α+β/2,  d1=d2=α,  d3=β,  and d4=0.Take the potential parametersD=88575.905cm-1, α=2.3425473Å-1, m1=22036a.u., m2=29376a.u., m=m1∗m2/m1+m2, and q0=1.128323Å to the CO molecule [5]. Initial condition q(0)=q0, p(0)=2MD-0.01, stepsize h=3, computing K=108 step, and integral interval T=[0,3∗108]. Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error; two atoms oscillate versus time and (q,p) classical trajectories of the CO diatomic systems as follows (Figures 1–6).Figure 1
2FE, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 2
2FE, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Figure 3
4SRK, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 4
4SRK, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Figure 5
4SS, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of CO molecule in the phase space (q,p).Figure 6
4SS, computingK=108 step, the evolution of the final 100000 points p(t),q(t).Take the potential parametersD=51443.117cm-1, α=24.251821nm-1, m1=25704a.u., m2=m1, m=m1/2, and q0=0.12203nm to the N2 molecule [5]. Initial condition q(0)=q0, p(0)=2MD-0.01, stepsize h=3, computing K=108 step, and integral interval T=[0,3∗108]. Respectively, we utilize 2FE, 4SRK, and 4SS to compute the energy error and (q,p) classical trajectories of the N2 diatomic systems as follows (Figures 7–9).Figure 7
2FE, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figure 8
4SRK, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figure 9
4SS, computingK=108 step, the final 100000 points’ corresponding energy error curve and phase trajectories of N2 molecule in the phase space (q,p).Figures1, 3, 5, 7, 8, and 9 show the evolution of the energy and the phase trajectories of CO, N2 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 (T=[0,3∗108]). The energy error computed by 2FE is only 10-15 (Figures 1 and 7) which indicates conservation of energy over long time intervals T=3∗108 under a large stepsize h=3, but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to 10-9 (4SRK, Figures 3 and 8) and 10-7 (4SS, Figures 5 and 9), which is only approximately conservative.


## 3.2.A2B Triatomic Molecules

Basing Banerjee and Adams constructed a transformation method; we consider the motion ofA2B triatomic molecule such as water (H2O) within the C2v symmetry; take the Cartesian coordinate system yOz, with origin at the center of mass O, and the z-axis is the C2 axis; the coordinates of the two A atoms and B atom are A(y1,z1), A(y2,z2), and B(y3,z3). The generalized coordinates q1=y1-y2, q2=z1-2z3+z2, and the generalized momenta p1=0.25dq1/dt, p2=0.5dq2/dt. Hamiltonian function of A2B type molecule is H(q,p)=K(p)+V(q), where K(p)=2p12+p22 is the kinetic energy and V(q) is the potential energy: (23)Vq=5π2D2-5D+6.5+4D-1+0.5π2q2-1.52+q2-1,D=q12+q221/2.The canonical differential equations are(24)dq1dt=∂H∂p1,dq2dt=∂H∂p2,dp1dt=-∂H∂q1,dp2dt=-∂H∂q2.Initial condition: q1(0)=3, q2(0)=3/2, p1(0)=0, p2(0)=0, and H(0)=50.1951. Taking stepsize h=0.01, computing K=108 step, integral interval T=[0,106].The second-order symplectic difference scheme (2SS) is(25)zk+1=zk+hJ-1Hzzk+1+zk2.We use the linear element method (1FE), general Leapfrog method [1], 2SS, and 4SS to compute the classical trajectories of A2B type molecule in phase spaces and energy error Hh-H as follows (Figures 10–13).Figure 10
1FE, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 11
General Leapfrog method, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 12
2SS, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.Figure 13
4SS, computingK=108 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.It is observed from Figures10–13 that the numerical results computed by the linear element methods are in agreement with the theoretical analysis; atom B and two A atoms in A2B 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 T=1000000, just as the symplectic difference method (Figures 11 and 12). The energy error computed by the linear element methods is only 10-10 when T=1000000, it reached the microscopic reaction needed to be considered time T=h∗K=106, but the energy deviation is comparatively larger by symplectic difference scheme, energy error up to 10-1 (Leapfrog method, Figure 11), 10-2 (2SS, Figure 12), and 10-3 (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.



---

*Source: 904140-2015-08-13.xml*