Higher-order symplectic integration techniques for molecular dynamics problems

Abstract

The explicit symplectic difference schemes with a number of stages from 1 to 5 are considered for the numerical solution of molecular dynamics problems described by systems with separable Hamiltonians. A general method for finding symplectic schemes of high order of accuracy using Gröbner bases is proposed. It is shown that it is possible to significantly increase the accuracy of symplectic schemes without increasing the number of stages by obtaining all possible real schemes with the aid of Gröbner bases. For the numbers of stages 2, 3, and 4, the parameters of Runge–Kutta–Nyström (RKN) schemes and partitioned Runge–Kutta (PRK) schemes are obtained using the Gröbner basis technique. It is shown that there exists only one explicit invertible (symmetric) two-stage RKN scheme. In the cases of the vanishing Vandermonde determinant, 20 new real four-stage RKN schemes and 18 new PRK schemes of Forest–Ruth and Yoshida were found in the analytic form with the aid of Gröbner bases. In the cases of the nonzero Vandermonde determinant, four symmetric four-stage RKN schemes were derived in the analytic form with the aid of the Gröbner basis technique. In addition, for the number of stages 4 and 5, new nonsymmetric schemes were found using the Nelder–Mead numerical optimization method. In particular, four new schemes were obtained for the number of stages 4. For the number of stages 5, three new schemes were obtained in addition to four schemes known earlier. For each specific number of stages, a scheme has been found that is the best in terms of the minimum of the leading term of the approximation error. The relations have been established between the polynomials entering the expansions of the errors of each scheme for the momentum and for the particle coordinate. Verification of the schemes was carried out on a problem that has an exact solution. It is shown that the symplectic five-stage RKN scheme provides a more accurate conservation of the total energy balance of the particle system than schemes of lower orders of accuracy. The stability studies of the schemes were performed using the Mathematica software package.