Using finite element codes as a numerical platform to run molecular dynamics simulations

F. Niederhöfer, J. Wackerfuß

A mathematically rigorous methodology for embedding the governing equations of molecular dynamics in the formalism of the finite element method is presented. Only one generalized finite element type is needed to cover all different types of existing interatomic potentials. The finite element type is simply specified by two parameters characterizing the type of the interatomic potential to be considered. Built on this formulation a partitioned Runge-Kutta method-summarizing a wide range of explicit and implicit, single- and multi-stage, lower and higher order time integration schemes-is embedded in a unified manner. The required finite element residual vector and the related Jacobian matrix are stated explicitly. The related FE-mesh coincides with the neighborhood lists used in standard molecular dynamics enabling the use of common tools. The range, versatility and performance of the proposed finite element formulation have been demonstrated by means of several numerical examples.

Today, molecular dynamics simulations are an important tool for characterizing known and developing new molecular structures in many areas of materials science, biotechnology and nanotechnology. The efficiency of such numerical simulations depends, among other things, on the choice of the underlying time integration method. Various factors must be taken into account when making this choice (total number of atoms, number of time steps, required accuracy, memory resources, etc.). On the other hand, influences must be taken into account that are independent of the choice of time integration method (e.g. calculation inaccuracies).

One aspect of this is the sensitivity of the numerical results to minimal changes in the initial values. The following animation illustrates this issue.

It shows the results of 3 different molecular dynamics simulations for a buckyball "pumped up" by 10 % compared to the equilibrium position. The Störmer-Verlet method with a time step size of Δt = 0.025 femtoseconds (deliberately small!) and a total of 1000 time steps was selected for time integration. The only difference between the 3 simulations is the accuracy of the initial data: In the left-hand simulation, these are specified to 4 decimal places, in the middle simulation to 8 decimal places and in the right-hand simulation to 12 decimal places. The visual comparison shows that the time of transition from a purely radial to a "chaotic" movement depends on the accuracy of the initial data. This behavior is independent of the selected time step method.

In another simulation example, a (8.0) carbon nanotube (CNT) with a length of approx. 4 nanometers is set into axial oscillation. The CNT, which is supported differently at the right edge, is excited at the left edge by a specified (axial) deflection. A total of 10,000 time steps were calculated with a time step width of Δt = 0.25 femtoseconds. The video sequence clearly shows a longitudinal wave propagating from left to right through the CNT, which is reflected at the right edge. This reflection initiates a successive transition to an indefinite grating oscillation.

Publication

Niederhöfer, F & Wackerfuß, J 2012, 'High-order time integration methods in molecular dynamics', PAMM Proceedings in Applied Mathematics and Mechanics, Volume 12, pp. 47-48.

Wackerfuß, J., Niederhöfer, F., 2018, Using finite element codes as a numerical platform to run molecular dynamics simulations. Computational Mechanics 18, 1-30, https://doi.org/10.1007/s00466-018-1594-5.