arXiv:cond-mat/0110585v1 [cond-mat.stat-mech] 29 Oct 2001 Optimized Forest-Ruth- and Suzuki-like algorithms for integration of motion in many-body systems I. P. Omelyan, 1,2 I. M. Mryglod, 1,2 and R. Folk 2 1 Institute for Condensed Matter Physics, 1 Svientsitskii Street, UA-79011 Lviv, Ukraine 2 Institute for Theoretical Physics, Linz University, A-4040 Linz, Austria (February 1, 2008) An approach is proposed to improve the efficiency of fourth-order algorithms for numerical in- tegration of the equations of motion in molecular dynamics simulations. The approach is based on an extension of the decomposition scheme by introducing extra evolution subpropagators. The extended set of parameters of the integration is then determined by reducing the norm of truncation terms to a minimum. In such a way, we derive new explicit symplectic Forest-Ruth- and Suzuki-like integrators and present them in time-reversible velocity and position forms. It is proven that these optimized integrators lead to the best accuracy in the calculations at the same computational cost among all possible algorithms of the fourth order from a given decomposition class. It is shown also that the Forest-Ruth-like algorithms, which are based on direct decomposition of exponential propagators, provide better optimization than their Suzuki-like counterparts which represent com- positions of second-order schemes. In particular, using our optimized Forest-Ruth-like algorithms allows us to increase the efficiency of the computations more than in ten times with respect to that of the original integrator by Forest and Ruth, and approximately in five times with respect to Suzuki’s approach. The theoretical predictions are confirmed in molecular dynamics simulations of a Lennard-Jones fluid. A special case of the optimization of the proposed Forest-Ruth-like algorithms to celestial mechanics simulations is considered as well. Pacs numbers: 02.60.Cb; 02.70.Ns; 05.10.-a; 45.50.Pk; 95.75.Pq Keywords: molecular dynamics, fourth-order algorithms, decomposition scheme, optimized algo- rithms, celestial mechanics simulations. I. INTRODUCTION Modelling various physical and chemical processes in molecular dynamics (MD) simulations we come to the ne- cessity to integrate the equations of motion for a many- body system of interacting particles. A lot of numeri- cal algorithms have been devised and implemented over the years to perform such an integration. The tradi- tional high-order explicit Runge-Kutta (RK) and implicit predictor-corrector (PC) schemes [1,2] were applied in early investigations. Is was soon realized that the ex- tra orders obtained in these schemes are not relevant since the truncation errors accumulate drastically on MD scales of time [3]. This high instability restricts the appli- cation of RK and PC integrators in long-term MD simu- lations to very small time steps only, and, thus, reduces significantly the efficiency of the computations. In addi- tion, the RK and PC algorithms produce solutions which, unlike exact phase trajectories, are neither symplectic nor time reversible. In 1990, a new approach to the integration of motion in many-body systems has been proposed [4–8]. Within this approach, the time propagation is carried out on the basis of exponential decompositions of evolution propaga- tors. The main advantage of the decomposition method is that for an arbitrary order in the time step it allows to construct algorithms which are exactly symplectic and time-reversible. The preservation of symplecticity and reversibility appears to be very important because, as is now well established, this closely relates to the stability of an algorithm [9]. Another nice property of the decom- position integration is its explicitness and simplicity in implementation. This is in a sharp contrast to implicit time-reversible symplectic algorithms obtained recently [10] within the RK approach, where cumbersome sys- tems of coupled nonlinear equations must be solved by iteration at each step of the integration process. Nowadays, the decomposition method should be con- sidered as the main tool for construction of efficient inte- grators of motion in classical as well as quantum systems [8]. Modified versions of this method have also been in- troduced. In particular, it was shown that for atomic systems with long-range interactions the efficiency of the integration can be improved by additionally splitting the Liouville operator into slow and fast components [11,12]. In such a multiple scale propagation, the slow subdynam- ics is treated in a specific way using larger step sizes in view of the weakness of long-range forces. The fasted motion, caused by the interactions at short interparti- cle distances, remains to be integrated with the help of original decomposition algorithms. The question of how to derive higher-order integrators by composing lower- order decomposition schemes has been considered as well [7,13,14]. Moreover, it has been shown recently how to 1