{'fraction_non_alphanumeric': 0.0465839246303899, 'fraction_numerical': 0.02157146586297729, 'mean_word_length': 4.305300264383734, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimized Forest-Ruth- and Suzuki-like algorithms for integration of motion in many-body systems

An approach is proposed to improve the efficiency of fourth-order algorithms for numerical integration 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 compositions 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.


I. INTRODUCTION
Modelling various physical and chemical processes in molecular dynamics (MD) simulations we come to the necessity to integrate the equations of motion for a manybody system of interacting particles. A lot of numerical algorithms have been devised and implemented over the years to perform such an integration. The traditional 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 extra 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 application of RK and PC integrators in long-term MD simulations to very small time steps only, and, thus, reduces significantly the efficiency of the computations. In addition, 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][5][6][7][8]. Within this approach, the time propagation is carried out on the basis of exponential decompositions of evolution propagators. 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 decomposition 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 systems of coupled nonlinear equations must be solved by iteration at each step of the integration process.
Nowadays, the decomposition method should be considered as the main tool for construction of efficient integrators of motion in classical as well as quantum systems [8]. Modified versions of this method have also been introduced. 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 subdynamics 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 interparticle distances, remains to be integrated with the help of original decomposition algorithms. The question of how to derive higher-order integrators by composing lowerorder decomposition schemes has been considered as well [7,13,14]. Moreover, it has been shown recently how to adapt the decomposition approach to integrate not only translational motion in atomic liquids, but also simulate more complicated molecular and spin systems at the presence of orientational degrees of freedom [15][16][17].
The main attention in previous studies has been directed to obtain algorithms which require a minimal number of force evaluations per time step. For instance, the well-known velocity-and position-Verlet integrators [11,18] correspond to a second-order decomposition scheme with one force evaluation per step. The fourth-order algorithm by Forest and Ruth [5] presents a scheme with three such force recalculations. Sixth-order schemes [4,7] can be reproduced beginning from seven evaluations of force for each particle during each time step. It is worth mentioning that sixth-and higher-order schemes involve too large number of force recalculations and generally are not recommended to be used in MD simulations. The reasons are that for a system with a great number of particles, the force evaluations constitute the most time-consuming part of the computations and that the truncations errors decrease with increasing the order of the scheme much slower than the total number of these evaluations. Such high-order algorithms may be efficient only for systems composed of a few bodies when a very high precision in determination of phase trajectories is desirable, for instance, in astronomical applications. In the case of MD simulations, the preference should be given to more simple second-or fourth-order schemes. For many MD applications to achieve a required level of accuracy, the fourth-order integrators appear to be more efficient than second-order schemes.
The minimal numbers of force evaluations per time step do not guarantee, however, the optimization with respect to the overall number of such evaluations which are necessary to perform during a fixed observation time interval. This is so because schemes, in which these numbers exceed the minimum, can be used with larger sizes of the time step in order to obtain the same accuracy in solutions. Only very few papers [7,19,20] were devoted to derivation of such extended schemes. In particular, Suzuki et al. have pointed out that the fourth-order algorithm by Forest and Ruth may not lead to optimal performance since during the evaluation it involves time coefficients which are larger in magnitude than the size of the initial time step [7]. As a result, the original algorithm has been replaced by a new fourth-order integrator composing five, instead of three, second-order Verlet schemes. However, the question of how this replacement influences on the efficiency of the computation has not been studied. Extended Suzuki-like schemes were also the subject of investigations by Kahan and Li [19]. Composing second-order algorithms, they introduced higherorder integrators with time coefficients chosen to provide a minimum for the maximal one among them in magnitude or for the sum of absolute values of these coefficients. Again, it has been unknown whether the increased number of force evaluations in their integrators is compensated by the possibility of using larger step sizes or not.
In Ref. [20] new explicit velocity-and position-Verlet-like algorithms of the second order were proposed to integrate the equations of motion in many-body systems. These algorithms were derived on the basis of an extended decomposition scheme at the presence of a free parameter. The nonzero value for this parameter was obtained by reducing the influence of truncated terms to a minimum. As a result, the new optimized algorithms appear to be more efficient than the original Verlet versions which correspond to a particular case when the introduced parameter is equal to zero.
In the present paper we propose a consequent approach for improving the efficiency of fourth-order decomposition integration. As a result, we derive new Forest-Ruthand Suzuki-like integrators by explicitly reducing the influence of truncation terms to a minimum. Although this requires (as in any extended scheme) extra force evaluations per time step, the overall number of force recalculations keeps unchanged due to the use of larger step sizes. At the same time, the resulting precision increases significantly with respect to that of the original algorithms by Forest and Ruth, and Suzuki. It is demonstrated also that previous criteria used for optimization of decomposition algorithms may have nothing to do with providing the best performance of the computations.

II. DECOMPOSITION INTEGRATION
We will deal with a classical N -body system described by the Hamiltonian where r i and v i denote the position and velocity, respectively, of particle i with mass m, and ϕ(r ij ) is the interparticle potential of interaction with r ij = |r i − r j |. The equations of motion for such a system can be cast in the following compact form Here ρ ≡ {r i , v i } is the full set (i = 1, 2, . . . , N ) of phase variables, [ • ] represents the Poisson bracket and is the Liouville operator with f i = − N j(j =i) ϕ ′ (r ij )r ij /r ij being the force acting on particles due to the interaction.
If an initial configuration ρ(0) is provided, the unique solution to Eq. (2) can be presented as where h denotes the time step, and the Liouville operator L = A + B has been split into the free-motion A = v·∂/∂r and potential B = f /m·∂/∂v parts with v ≡ {v i }, r ≡ {r i }, and f ≡ {f i }. The significance of such a splitting will be understood below. Then the evolution of the system can be investigated during arbitrary times t by repeating the single-step propagation, Of course, solution (4) is quite formal because the exponential propagator e Lh does not allow to be evaluated exactly at any h (solutions in quadratures are possible only for N = 2 that is not relevant for our case of manybody systems when N ≫ 1). However, at small enough values of h, the total propagator can be decomposed [4][5][6][7][8] using the formula The coefficients a p and b p in this formula should be chosen in such a way to provide the highest possible value for K ≥ 1 at a given integer number P ≥ 1. Then integration (4) can be performed approximately with the help of Eq. (5) by neglecting truncation terms O(h K+1 ). The main advantage of the above decomposition is that the exponential subpropagators e Aτ and e Bτ , appearing in the right-hand-side of Eq. (5), are analytically integrable. Indeed, as can be shown readily, e Bτ ρ ≡ e f /m·∂/∂vτ {r, v} = {r, v + τ f /m} with τ being equal to a p h or b p h. That is why we explicitly presented L as the sum of free-motion A and potential B terms. Another important feature of the decomposition integration is that it leads to symplectic trajectories. This is so since Eq. (6) represents, in fact, simple shifts of position and velocity, and these shifts do not change the volume in phase space. The time reversibility S(−t)ρ(t) = ρ(0) of solutions (following from the property S −1 (t) = S(−t) of evolution operator S(t) = e Lt ) can also be reproduced by imposing additional constraints on the coefficients a p and b p , namely, a 1 = 0, a p+1 = a P −p+1 , and b p = b P −p+1 , or a p = a P −p+1 , and b p = a P −p at b P = 0. Then the subpropagators e Aτ and e Bτ will enter symmetrically in the decompositions and, thus, provide automatically the required reversibility. Note also that such constraints lead to automatic disappearing even-order terms ∝ h 2k in the function O(h K+1 ) for any k ≥ 0. For this reason, the order K of time-reversible algorithms may accept only even numbers, (K = 2, 4, 6, . . .). The cancellation of odd-order terms ∝ h 2k+1 in O(h K+1 ) up to a required finite number k will be provided by fulfilling a set of basic conditions for a k and b k . For example, the condition P p=1 a p = P p=1 b p = 1 is necessary to cancel the firstorder truncation uncertainties.
The decomposition method is quite general to build numerical integrators of arbitrary orders. In particular, the second-order (K = 2) velocity-Verlet (VV) algorithm [11,18] is obtained from Eq. (5) at P = 2 and a 1 = 0, b 1 = b 2 = 1/2, a 2 = 1, i.e., e (A+B)h+O(h 3 ) = e Bh/2 e Ah e Bh/2 . The case when the operators A and B are replaced by each other (A ↔ B) is also possible, and we come [11] to the position-Verlet (PV) integrator corresponding to the choice a 1 = a 2 = 1/2, b 1 = 1, and b 2 = 0. The fourth-order (K = 4) algorithm by Forest and Ruth (FR) [5] is immediately reproduced from Eq. (4) at P = 4 (8) can be related to the position version of the FR integrator (PFR), because putting formally θ = 0 it transforms to the second-order PV algorithm (7). For A ↔ B, Eq. (8) will represent the velocity FR counterpart (VFR), which can be derived directly from Eq.

III. OPTIMIZATION OF FOURTH-ORDER ALGORITHMS WITHIN DIRECT DECOMPOSITION
Let us consider now an extended decomposition scheme of the fourth order by allowing to accept a value for P which exceeds the necessary minimum (P = 4) on unity, i.e., letting P = 5. Remember that we cannot choose the number P to be too big, because this results in too larger number, namely P − 1, of expensive force recalculations. Chossing P = 5 we hope simply to reduce the truncation errors O(h 5 ) significantly in a little additional computation cost, rather than to increase the order of the decomposition scheme (note that sixth-order integrators are derivable [4] beginning up from P = 8).
For P = 5, the extended decomposition can be presented in the form Here the symmetry of time coefficients and the condition 5 p=1 a p = 5 p=1 b p = 1 have already been taken into account. Again, the propagation with A ↔ B is also acceptable, and then a 1 = a 5 = ξ, , and b 5 = 0. The operator C 3 , appearing in the left-hand-side of Eq. (9), is responsible for the cancellation of third-order truncation uncertainties. The explicit expression for it is where [ , ] denotes the commutator of two operators, and So that at C 3 = 0, i.e. when α(ξ, λ, χ) = 0 and β(ξ, λ, χ) = 0, formula (9) represents a whole family of symplectic time-reversible integrators of the fourth-order. A particular member of the above family can be obtained by choosing corresponding values for ξ, λ, and χ. As far as there are three parameters and only two constraints, α = 0 and β = 0, one from these parameters, ξ say, can be treated to be free. Then, for example, putting ξ = 0, Eq. (9) reduces to the original FR algorithm (8) in position or (when A ↔ B) velocity forms. The extended (when ξ = 0) propagation will require already four, instead of three, force recalculation per time step. However, having a room in varying ξ, we can overcompensate the increased computational efforts by minimizing the fifth-order truncation uncertainties C 5 h 5 .
In order to show this, let us analyze in detail the influence of these uncertainties on the result. Expanding both the sides of Eq. (9) into Taylor's series with respect to h, one finds where explicit expressions for γ-multipliers are: Assuming that all the fifth-order commutators in (12) are nonzero valued, the norm of C 5 with respect to fifthorder commutators arising in Eq. (12) can be written as follows Then the norm of local uncertainties C 5 ρh 5 appearing in phase trajectory ρ during a single-step propagation given by Eqs. (4) and (9) can be expressed in terms of γ and h as g = γh 5 . During a whole integration over a fixed time interval t, the total number l of such single steps is proportional to h −1 . As a result, the local fifth-order uncertainties will accumulate step by step leading at t ≫ h to the fourth-order global errors Γ = gh −1 , i.e., Extended propagation (9) can now be optimized with respect to time coefficients ξ, λ, and χ by finding the global minimum for the function γ(ξ, λ, χ), provided α = 0 and β = 0, i.e., solving the system of equations A way to simplify the problem is to solve analytically the first two equations of (15) with respect to χ and ξ, then substitute expressions (16) into the third equation, and find numerically the minimum considering already γ(ξ(λ), λ, χ(λ)) as a function of only one variable λ. The result within sixteenth significant digits is The global minimum of γ(ξ, λ, χ), corresponding to solutions (17) found (they satisfy Eq. (16) at sign "+" in the right hand side of expression for χ(λ)), consists γ EFRL min ≈ 0.00092, where the superscript EFRL refers to the extended FR-like integration (9). On the other hand, the value of γ corresponding to usual FR scheme (8), i.e. when ξ = 0, λ = (1 − θ)/2 and χ = θ (then Eq. (9) reduces to Eq. (8)), is equal to γ FR ≈ 0.039. We see, therefore, that applying the extended integration allows one to decrease the truncation errors approximately in γ FR /γ EFRL min ≈ 42 times. Taking into account that such a decreasing has been achieved increasing the number of force evaluations per step from three to four, the extended propagation must be performed with step sizes which are in factor 4/3 higher than those of the FR algorithm, in order to provide the same number of total force recalculations during the fixed overall interval of integration. Thus, we will come to more efficient calculations if the following inequality Γ EFRL min (4h/3) < Γ FR (h) takes place. In view of Eq. (14), such an inequality can be rewritten as γ FR /γ EFRL min > (4/3) 4 ≈ 3.16, so that it is fulfilled completely in the optimization regime. In particular, indicating that the global errors can be reduced more than in 10 times with respect to the FR integration without spending any additional overall computational costs. The proposed procedure can be used, in principle, for the decompositions of arbitrary two noncommutative operators A and B to achieve the best performance. Note that the necessity in these decompositions may arise not only when considering the integration of motion in classical systems, but also in mathematical and quantum mechanical calculations. In some cases, further improvement of the efficiency of the computations may also be possible. For instance, an extra optimization of the decomposition scheme can be achieved in celestial mechanics applications due to a specific character of motion in the solar system (see Appendix).
In the case of MD simulations with velocity-independent forces, an additional optimization can be carried out as well using specific properties of operators A and B. Taking into account explicit expressions for these operators, it can be verified readily that two of the sixth five-order commutators vanish in Eq. (12) (9). This scheme, by analogy to Verlet integrators, we will refer to the velocity version of the EFRL algorithm and abbreviate in our notations as VEFRL. Then, letting formally γ 4 = 0 and γ 6 = 0 to exclude the above zeroth commutators in Eq. (13), and resolving problem (15) on condition minimum yields This leads to the global minimum γ VEFRL min ≈ 0.00065 (which is achieved at sign "−" in Eq. (16)), while the value of γ corresponding to usual FR integration (8) reduces (at γ 4 = 0 and γ 6 = 0) to γ FR ≈ 0.028. Again, the truncation uncertainties decrease approximately in the same factor, γ FR /γ VEFRL where the values for ξ, λ, and χ should be taken from Eq. (19). The optimized PEFRL algorithm (when A ↔ B in Eq. (9)) reads: and the parameters ξ, λ, and χ should accept their values from Eq. (20). The algorithms are explicit, simple in implementation, and require only slight modification with respect to the original FR scheme.

IV. OPTIMIZATION BY COMPOSING SECOND-ORDER SCHEMES
Another way to construct high-order algorithms consists in composing lower-order schemes. In particular, employing second-order Verlet integrator (7), the composition can be performed [7] for an arbitrary higher order K > 2 as Here the coefficients d q are chosen by providing a maximum for K at a given number Q > 1. This leads to the necessity of fulfilling a set of order constraints, and, for instance, the condition Q q=1 d q = 1 must be satisfied to avoid the first-order truncation terms in O(h K+1 ). For time-reversible compositions, the coefficients d q should appear symmetrically, i.e., satisfy the property d q = d Q−q+1 . Then, as in the case of direct decomposition (5), even-order terms will be absent in the function O(h K+1 ), and the order of reversible composition schemes will also accept only even numbers.
Fourth-order (K = 4) composition integrators are derivable from Eq. (23) at Q ≥ 3. For example, putting Q = 3 as well as d 1 = d 3 = θ and d 2 = 1 − 2θ one comes to the scheme which coincides entirely with the FR algorithm (as can be seen by comparing with Eq. (8)). Thus, in this particular situation, the direct decomposition and secondorder-based composition approaches exhibit to be identical (but with increasing K or Q both the approaches will lead to different results). We will now consider the construction of extended (when Q > 3) composition algorithms of the fourth order. It can be shown that no real solutions exist for d q at Q = 4 and K = 4. So putting Q = 5, one obtains the simplest extended fourth-order composition formula with α(ξ, λ) = 2ξ 3 + 2λ 3 + 1 − 2(ξ + λ) 3 .
Formula (25) constitutes a family of composition timereversible integrators of the fourth-order, provided C 3 = 0, i.e. α(ξ, λ) = 0. Any one from the two parameters ξ and λ can be chosen, in principle, arbitrarily since we have only one constraint. For ξ = 0 and λ = θ, we come to usual FR integration (8). Suzuki [7], considering an extended scheme like (25), has imposed the additional constraint ξ = λ ≡ ϑ, and obtained ϑ = 1/(4 − 3 √ 4) ≈ 0.41449. We will show below that although Suzuki's approach leads to an increased efficiency with respect to the original FR scheme, it is not the best choice for fourth-order integration.
The operator C 5 , which forms the fifth-order term of truncation uncertainties in composition (25), can again be presented in form of Eq. (12), where now It is worth remarking that the operators C 3 and C 5 cannot be reduced to zero simultaneously at Q = 5, since we cannot satisfy all the three equations α = 0, s = 0, and w = 0 having only two parameters ξ and λ (such a reduction, resulting in sixth-order composition integrators, is possible within real coefficients d q beginning up from Q = 7 [7]). Our task is to minimize the norm γ(ξ, λ) = γ 2 1 + γ 2 2 + γ 2 3 + γ 2 4 + γ 2 5 + γ 2 6 (29) of C 5 in space of fifth-order commutators, provided C 3 = 0, i.e., to solve the problem The simplest way to do this is to transform the first of two equalities of system (30) from the cubic (see Eq. (27)) to square equation replacing the variable ξ by the new independent quantity χ = 1 − 2(ξ + λ), and find solutions for λ considering χ as a parameter. As a result, we obtain Then, substituting these expressions into the second equation, we come to the function γ ξ(χ), λ(χ) ≡ γ(χ) which depends already on only one variable χ. The global minimum of this function is achieved at χ = −0.7269082885036828E + 00 (with sign "+" for λ(χ) in Eq. (31)) and consists γ ESL min ≈ 0.0011, where the superscript ESL refers to the extended Suzuki-like integration (25). So that according to Eq. (31), the corresponding solutions are ξ = 0.3221375960817984E+00 (32) λ = 0.5413165481700430E+00 .
As was mentioned above, Suzuki has used the parameters ξ = λ = ϑ ≈ 0.41449 in his integration. This corresponds to the value γ S ≈ 0.0015 of function (29) that is approximately in factor 1.5 higher than its optimized counterpart γ ESL min ≈ 0.0011. Therefore, his choice was not optimal. But the main conclusion we may make in this context is that the original Suzuki approach as well its optimized ESL version are less efficient than the EFRL integration described in the preceding section. Indeed, taking into account the value γ EFRL min ≈ 0.00092 and the fact that the Suzuki approach requires up five (see Eq. (25)), instead of four, force evaluations per time step, we obtain in view of Eq. (14) that efficiency of the EFRL integration is in factor better with respect to the ESL scheme, and, thus, approximately in 5 times higher with respect to the original approach by Suzuki. Note that results (32) and (33) correspond to a general case of integration of motion when forces may explicitly depend, rigorously speaking, on velocities. for which γ PS ≈ 0.0014, the efficiency increases more than in 5 times. It is necessary to point out that solutions (34) correspond to a local minimum of γ (when γ 4 = 0 and γ 6 = 0). The global minimum is achieved at ξ = 1.453335388755048, λ = −2.154220603244978, with χ = 2.401770428979862 and consists about 0.00096 that is only slightly smaller than the presented above value 0.00109. But now all the time coefficients appear to be significantly larger than their counterparts (34). In such a situation, the preference should be given to the local minimum, because it leads to a significantly smaller value for the norm of higher-order truncation uncertainties O(h 7 ), and thus the precision of composition (25) will be much less sensitive to increasing the size of the time step h (see the next section).

V. APPLICATION TO MOLECULAR DYNAMICS SIMULATIONS
In order to verify theoretical predictions presented in sections III and IV, we have applied the velocity and position versions of our EFRL and ESL algorithms to MD simulations of a Lennard-Jones fluid (LJ), and compared the results with those of the original FR and Suzuki's approaches. The system under consideration was a collection of N = 256 particles interacting through a shifted LJ-like potential, ϕ(r) = Φ(r) − Φ(r c ) at r < r c and ϕ(r) = 0 otherwise, where Φ(r) = 4u (σ/r) 12 − (σ/r) 6 is the genuine LJ potential. The particles were placed in a basic cubic box of volume V = L 3 , and the modification of Φ(r) with r c = L/2 ≈ 3.36σ as well as the periodic boundary conditions have been used to exclude the finite-size effects. The simulations were performed in a microcanonical ensemble at a reduced density of n * = N V σ 3 = 0.845 and a reduced temperature of T * = k B T /u = 1.7. All runs of the length in l = 10 000 time steps each were started from an identical well equilibrated initial configuration ρ(0). The precision of the algorithms was measured in terms of the relative total energy fluctuations E = ( N j(j =i) ϕ(r ij ) and denotes the microcanonical averaging. Note that if the equations of motion could be solved exactly, the above fluctuations should vanish, because in microcanonical ensembles the total energy is an integral of motion, E(t) = E(0). So that during approximate MD integrations, smaller values of E will correspond to a better precision in evaluation of phase trajectories. The equations of motion were first integrated with the help of extended decomposition scheme (9) in which the parameters ξ, λ, and χ, being constant within each run, varied from one run to another. For convenience of presentation, we have chosen ξ (instead of λ as earlier) to be a free parameter, so that the two other quantities λ(ξ) and χ(ξ) should be treated as depending on ξ according to constraints (16). The total energy fluctuations obtained in such simulations at the end of the runs for three (fixed within each run) undimensional time steps, namely h * = h(u/mσ 2 ) 1/2 = 0.00125, 0.0025, and 0.005, are shown in Fig. 1 as functions of parameter ξ. The subsets (a) and (b) of this figure correspond to VEFRL (Eq. (21)) and PEFRL (Eq. (22)) versions, respectively. As can be seen, all the dependencies E(ξ, h) have the global minimum which locates at the same point, ξ ≈ 0.164 for VE-FRL or ξ ≈ 0.179 for PEFRL, independently on the size h of the time step. This point coincides completely with the minimum given by Eq. (19) or (20) for the function γ(ξ) ≡ γ(ξ, λ(ξ), χ(ξ)) (see Eq. (13), in which we should put γ 4 = 0 and γ 6 = 0 for VEFRL or γ 1 = 0 and γ 3 = 0 for PEFRL) which is also included in Fig. 1 (dashed  curves in the subsets). Moreover, the fluctuations E(ξ, h) appear to be proportional to the norm Γ = γh 4 (Eq. (14)) of global errors, and the coefficient of this proportionality almost does not depend on ξ and h. In addition, at each time step the energy fluctuations decrease at the minimum in 40-60 times (in dependence on the versions considered) with respect to those at ξ = 0, that is in agreement with our predicted values γ FR /γ VEFRL min ≈ 43 and γ FR /γ PEFRL min ≈ 62. The result for the total energy fluctuations as functions of the length l = t/h of the simulations corresponding to the optimized VEFRL (Eqs. (19) and (21))) and PEFRL (Eqs. (20) and (22)) algorithms is presented in subsets (a), (b), (c), and (d) of Fig. 2 for the time steps h * = 0.00125, 0.0025, 0.005, and 0.01, respectively. For the purpose of comparison, the same dependencies corresponding to original VFR and PFR integrators (see Eq. (8)) are also shown there. Note that for the original integrators, the time step for each subset was chosen to be always in factor 4/3 smaller than that of the optimized versions, namely, h * = 0.0009375, 0.001875, 0.00375, and 0.0075. The last condition is necessary to provide the same number of total force recalculations during the same observation time t ≫ h within both the approaches. Note also that within these approaches, the fifth-and higherorders truncation uncertainties become too big at step sizes h * > 0.01. In particular, then the ratio of the total energy fluctuations to the fluctuations in potential energy (the standard quantity to estimate the precision of the calculations) exceeds a few per cent. Such large step sizes cannot be used in precise MD simulations and, for this reason, are not considered here. As can be seen from Fig. 2, both the original (VFR and PFR) and optimized (VEFRL and PEFRL) algorithms exhibit very good stability properties, that is explained by the symplecticity and time reversibility of the produced solutions. No systematic deviations in the total energy fluctuations can be observed for all the integrators. Instead, in each of the cases the amplitude of the fluctuations quickly tends to its own value which does not increase with further increasing the length of the simulations. However, this value is significantly larger within the original VFR and PFR integration. At the same time, using the optimized VEFRL and PEFRL algorithms allows us to decrease the unphysical energy fluctuations approximately in factor 10-15 (for VEFRL) or 20-25 (for PEFRL), despite the larger time steps. This is in excellent accord with our theoretical predictions Γ VEFRL min (4h/3)/Γ FR (h) ≈ 0.073 and Γ PEFRL min (4h/3)/Γ FR (h) ≈ 0.051, obtained in section III. It is worth pointing out also that the PEFRL algorithm is slightly better in energy conservation than its VEFRL counterpart (whereas the VFR integrator is better with respect to its PFR version).
Consider now the results for composition integration. The total energy fluctuations obtained in this case at the end of the simulation runs for the same set (as in Fig. 1) of undimensional time steps, h * = h(u/mσ 2 ) 1/2 = 0.00125, 0.0025, and 0.005, are shown in Fig. 3 as depending on parameters ξ (subsets (a) and (c)) and χ ((b) and (d)). The subsets (a)-(b) and (c)-(d) relate to the VESL and PESL versions, respectively, of composition scheme (25). Mention that the time coefficients ξ(χ) and λ(χ) for this scheme should satisfy equality (31) in which the quantity χ = 1 − 2(ξ + λ) is treated as a free parameter. So that only one χ-dependency is enough, in principle, to present the results. But, in view of somewhat complicated structure of function E(χ), we have included the ξ-dependency as well for more clarity. Note also that according to constraint (31), the parameter χ cannot accept values from the interval ]1 − 4ϑ ≈ −0.658, 1], because otherwise this will lead to imaginary solutions for ξ and λ. As can easily be seen from Fig. (3), all the dependencies have up three minima and their locations on ξ-and χ-axes practically do not depend on the size h of the time step. Moreover, the fluctuations E(ξ, h) again appear to be proportional to the norm Γ = γh 4 of global errors, where the function γ is now defined by Eq. (29) and plotted in the figure as well. Among the three minima, two of them have nearly the same depth, and for VESL integration the global minimum of E is achieved at ξ ≈ 0.323 or χ ≈ −.726, i.e., from the left of the interval ]1 − 4ϑ, 1] for χ. This coincides with the predicted values given by Eq. (35) corresponding to the global minimum of γ. In the case of PESL integration the pattern is somewhat different. At small enough values of the time step, namely at h ≤ 0.0025, the global minimum for the functions E and γ is identified from the right of the above χ-internal, namely, at ξ ≈ 1.45 or χ ≈ 2.4. But with further increasing h, the position of this minimum for E begin to shift with respect to its initial value. Such a behavior can be explained by the influence of higher-order O(h 7 ) truncation errors, which appears to be significant in this case due to the largeness of absolute values of time coefficients ξ, λ, and χ (in particular, then all these coefficients are larger in amplitude than unity). This is contrary to a left-lying minimum of E which is achieved at the same point ξ ≈ 0.316 or χ ≈ −0.737 as that of function γ (see Eq. ((34) independently on the size of the time steps considered. Then taking into account that both the minima have almost the same depth, and the fact that the left-lying minimum of E, being local at small h, transforms into the global minimum at moderate and larger time steps, solutions given by Eq. (34) should be considered as optimal.
As expected, at the minima found the energy fluctuations decrease approximately in 1.5 times with respect to those at ξ = λ ≡ ϑ ≈ 0.41 (and χ ≈ −0.66) corresponding to the original Suzuki approach. It is interesting to point out that the constraint ξ = λ minimizes the function 5 q=1 |d q | ≡ 2|ξ| + 2|λ| + |χ| and provides a minimum for the maximal value among the quantities |ξ|, |λ|, and |χ|. Using these last two criteria and the additional requirement that |d q | < 1, Kahan and Li [19] have derived extended composition schemes like (23) up to the tenth order in the time step. In view of our more precise analysis performed in this paper we see, therefore, that the previous criteria for optimization of algorithms may have nothing to do with providing the best performance of the computations. Finally, comparing the results presented in Figs. 1 and 3 between themselves, it can be pointed out that the values for minima of E corresponding to the VESL and PESL integration are approximately in 1.5 times higher than those of the VEFRL and PEFRL algorithms. In addition, VESL and PESL integrators require more (five, instead of four) force evaluations per time step, so that their efficiency is in factor 1.5 (5/4) 4 ≈ 3 worse with respect to the VEFRL and PEFRL schemes. Therefore, as was predicted (see. Eq. (33)), the preference should be given to the decomposition VEFRL and PEFRL algorithms rather than to their composition counterparts VESL and PESL.

VI. CONCLUSION
We have formulated a new approach for improving the efficiency of decomposition integration of the equations of motion in many-body systems. As a result, the Forest-Ruth-like algorithms have been derived and optimized for MD as well as celestial mechanics simulations. The algorithms are explicit, simple in implementation, and produce phase trajectories which are exactly symplectic and time reversible. Their main advantage over the original integrators by Forest and Ruth is the possibility to generate much more precise solutions at the same overall computational efforts. It has been predicted theoretically and confirmed in actual MD simulations of a LJ fluid that the optimized Forest-Ruth-like algorithms allow to decrease significantly the unphysical fluctuations of the total energy using even greater time steps. The question of how to optimize the integration by composing second-order schemes has also been studied. As has been shown, the Suzuki-like algorithms, obtained in this case, although improve the efficiency of the computations with respect to the original integrators by Suzuki as well as Forest and Ruth, appear to be less efficient than the optimized Forest-Ruth-like algorithms. The last algorithms should be considered as optimal among all possible decomposition integrators of the fourth-order with single splitting of the Liouville operator.
The approach presented can also be adapted for optimization of the integration of motion in more complicated systems for which splitting of the Liouville operator into more than two parts may pay dividends. The examples are multi-component fluids and systems with longrange interactions, where characteristic time intervals of the dynamical processes can differ by many orders from each other. In some cases, for instance, for systems with orientational degrees of freedom, the additional splitting may be necessary to obtain analytically integrable parts. These and related problems will be considered in a separate publication.
Part of this work was supported by the Fonds zur Förderung der wissenschaftlichen Forschung under Project No. P15247.

Appendix. Optimization of Forest-Ruth-like algorithms in specific cases
In section III, we have derived the velocity-and position-Forest-Ruth-like algorithms and optimized them for integration of motion in MD simulations. Now, one considers two specific cases when the efficiency of decomposition scheme (9) can be improved additionally. One example is a collection of weakly interacting particles, in which the potential part B ≡ εB of the Liouville operator can be treated as a small perturbation with respect to its kinetic counterpart A, where ε ≪ 1. Another case is the integration of motion in the solar system, where the interactions between planets are negligible small in comparison with the force that acts on them with respect to the sun. Then, the Liouville operator can be presented as L = A + εB, where now A is the sum of two-bodies Liouville's commutative operators, each of which describes the motion of an isolated body in an external field of the central mass, whereas εB corresponds to weak interactions of the bodies between themselves. Here, as in the usual case, the separate exponential operators e Aτ and e Bτ appear to be exactly integrable, so that the decomposition scheme will also lead to explicit integration. Therefore, assuming that the parameter ε is sufficiently small and considering the velocity version of extended Forest-Ruth-like decomposition (9), we can neglect by fifth-order commutators arising in the truncation term C 5 at the multipliers γ 2 , γ 3 , γ 4 , γ 5 , and γ 6 (see Eq. (12)). The reason is that such commutators are of order ε 2 or higher, whereas the commutator at γ 1 is proportional only to the first power of ε. In such a case, it is quite natural to reduce the multiplier γ 1 to zero exactly, rather than to find a minimum for the norm (13) of γ with respect to all the commutators. Then the problem (15) transforms into the following system of three equations      α(ξ, λ, χ) = 0 , β(ξ, λ, χ) = 0 , γ 1 (ξ, λ, χ) = 0 , where α and β are defined according to Eq. (11), and the expression for γ 1 follows immediately after Eq. (12). Since there exist up four sets of real solutions to Eq. (A1), an additional requirement should be imposed to choose the optimal values among them. Obviously, such values must minimize the norm of the main term of truncated uncertainties, which now is equal to γ 2 2 + γ 2 3 (note that commutators at γ 2 and γ 3 are proportional to ε 2 , while to the third power of ε or higher at γ 4 , γ 5 , and γ 6 ). As a result, one finds the optimal solutions ξ = 0.5316386245813512 λ = 0.5437514219173741 (A2) χ = −0.3086019704406066 . A ↔ B), the commutators at γ 1 and γ 4 are interchanged, so that it is necessary to solve the system α = 0, β = 0, and γ 4 = 0. No real solutions to this system have been found. In such a situation, we could try to find the global minimum for |γ 4 |.

For the position version (when
But since it appears to be greater than zero, the position version is not recommended to be used for this case and the preference should be given to the extended integration in velocity form (with time coefficients defined by Eq. (A2)). Note that the main terms of uncertainties for this integration are O(ε 2 h 5 ) and O(εh 7 ), so that in view of the smallness of ε, the scheme derived can be related to a quasi sixth-order algorithm. Such an algorithm, contrary to usual sixth-order schemes, will require four, instead of seven [4,7], force evaluations per time step.
In the above derivation, we tentatively assumed that the order of smallness of ε is nearly the same as that of the size for the time step h. In systems for which the parameter ε is extremely small, the strategy of optimization should be reconstructed. Namely, we can now ignore all the terms of the second order with respect to ε, even that which appears at β (see Eq. (10)) in thirdorder uncertainties O(h 3 ) with respect to h. Then within the velocity version, one comes to two equations, α = 0 and γ 1 = 0, for three variables, ξ, λ, and χ. Thus, the order of the decomposition scheme can be increased additionally by canceling an extra-order term O(εh 7 ) (which is assumed to be much more significant than the terms O(ε 2 h 3 ) and O(ε 2 h 5 )). An explicit expression for O(εh 7 ) (it is denoted simply as O(h 7 ) in decomposition formula (9)) has the form