Continuous Finite Element Methods of Molecular Dynamics Simulations

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


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 , 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][2][3][4][5][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][8][9][10][11][12][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][5][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][17][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  2  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.

Continuous Finite Element Methods for Hamiltonian Systems
Consider the first-order ordinary differential equation with initial value in the interval  = [0, ]: Let 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  + 1 parameters.As it is an initial value problem, and we already knew the initial value ( −1 ) 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  −1 .Taking V ∈  ℎ , then its derivate V  ∈  −1 .In practical computation, take V = ( −  −1 )  ,  = 0, 1, 2, . . .,  − 1.
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  2 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, ( +1 /  )  ( +1 /  ) =  + (ℎ 3 ).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, ( +1 /  )  ( +1 /  ) =  + (ℎ 5 ).It also possesses fourth-order accuracy and preserves the energy.
The specific calculation format of the linear continuous finite element (1FE The specific calculation format of the quadratic continuous finite element (2FE We utilize high accuracy numerical integration such as at least  + 1 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   = ( −1 ,   ), and then we can get the value   ,  = 1, . . ., .

Numerical Experiments
is the canonical mass.Then  =   and  =  2 /(2) 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 () = ( −2(− 0 ) −  −(− 0 ) ) [6].In each case, the parameters , , and  0 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, 4thorder implicit symplectic Runge-Kutta scheme (4SRK): The 4-stage, 4th-order explicit symplectic difference scheme (4SS) is where ].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).
Figures 1, 3, 5, 7, 8, and 9 show the evolution of the energy and the phase trajectories of CO,  2 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 ( = [0, 3 * 10 8 ]).The energy error computed by 2FE is only 10 −15 (Figures 1 and 7) which indicates conservation of energy over long time intervals  = 3 * 10 8 under a large stepsize   ℎ = 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.

𝐴 2 𝐵 Triatomic Molecules.
Basing Banerjee and Adams constructed a transformation method; we consider the motion of  2  triatomic molecule such as water (H 2 O) within the  2V symmetry; take the Cartesian coordinate system , with origin at the center of mass , and the axis is the  2 axis; the coordinates of the two  atoms and  atom are ( 1 ,  1 ), ( 2 ,  2 ), and ( 3 ,  3 ).The generalized  coordinates  The canonical differential equations are (24)   The second-order symplectic difference scheme (2SS) is

Initial condition: 𝑞
We use the linear element method (1FE), general Leapfrog method [1]   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  = 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  = 1000000, it reached the microscopic reaction needed to be considered time  = ℎ *  = 10 6 , 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.

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.

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

Figure 11 :Figure 12 :
Figure 11: General Leapfrog method, computing  = 10 8 step, the final 1000 nodes classical trajectories in phase space and corresponding energy error curve.
3.1.Diatomic Molecular System.Consider the classical motion of  diatomic molecule in the electronic potential, where the mass of atom  is  1 and that of atom  is  2 .The coordinates of  and  are  1 and  2 , respectively. =  =  1 +  2 is the canonical coordinate; 2FE, computing  = 10 8 step, the final 100000 points' corresponding energy error curve and phase trajectories of CO molecule in the phase space (, ).
, 2SS, and 4SS to compute the classical trajectories of  2  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  2  type molecule vibrate quasi-periodically