Processed Splitting Algorithms for Rigid-Body Molecular Dynamics Simulations
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y Processed Splitting Algorithms for Rigid-Body Molecular Dynamics Simulations
Igor P. Omelyan
Institute for Condensed Matter Physics, 1 Svientsitskii Street, UA-79011 Lviv, Ukraine andInstitute for Theoretical Physics, Linz University, A-4040 Linz, Austria (Dated: November 15, 2018)A new approach for integration of motion in many-body systems of interacting polyatomic molecules isproposed. It is based on splitting time propagation of pseudo-variables in a modified phase space, while the realtranslational and orientational coordinates are decoded by processing transformations. This allows to overcomethe barrier on the order of precision of the integration at a given number of force-torque evaluations per timestep. Testing in dynamics of water versus previous methods shows that the obtained algorithms significantlyimprove the accuracy of the simulations without extra computational costs.
PACS numbers: 02.60.Cb, 02.70.Ns, 05.10.-a, 45.40.-f
I. INTRODUCTION
Systems of rigid bodies are widely used to model variousphenomena on a broad range of length scales: from the mi-croscopic dynamics of molecules in gases and liquids [1, 2],mesoscopic behavior of polymers and other complex collec-tions in chemical and biological physics [3, 4] to macro-scopic movement of astrophysical objects in celestial me-chanics [5, 6]. A lot of approaches, including the traditionalRunge-Kutta and predictor-corrector schemes [1] as well asmore recent splitting techniques [7, 8, 9, 10, 11, 12], havebeen devised over the years to integrate the rigid-body equa-tions of motion.Now it is well established that the most adequate integra-tion can be done by splitting the time propagator into analyt-ically solvable parts [13, 14, 15]. For Hamiltonian systemsthis provides the preservation of such essential properties asconservation of volume in phase space and time reversibility.As a result, the splitting algorithms exhibit remarkable stabil-ity and thus are ideal for long-duration molecular dynamics(MD) simulations. In addition, these algorithms can be sym-plectic, i.e. can exactly conserve the total energy associatedwith a nearby Hamiltonian.The splitting approach however has a limitation on the or-der K of precision at each given number n of force-torqueevaluations per time step. Note that these evaluations presentthe most time-consuming part of the propagation. For thisreason, the rigid-body motion in MD simulations is inte-grated mainly by the simplest ( K = 2 ) Verlet-type algorithms[8, 9, 10, 11, 12] with n = 1 . The optimized algorithms[13, 14, 15] at n = 2 can outperform Verlet schemes. Butsuch an optimization does not rise the order of precision andfor K = 2 only modest accuracy can be reached. Higher-order ( K = 4 ) splitting schemes (note that K should be evento ensure time reversibility) can be derived beginning from n = 3 [14, 15]. The grown computational costs at n = 3 and K = 4 can be compensated by the increased precision whenadding gradient-like terms to the splitting propagator [15].Meanwhile it has been found that the order K of precisioncan be risen by carrying out supplementary (so-called pro-cessed) decompositions apart from the basic (kernel) splitting[16]. For K = 4 , each minimal kernel and processor leadsto one force and one force-gradient evaluations. This yields an effective number n = 2(1 + ν ) , where ν is the relativecost spent on the gradient evaluation with respect to that onthe force calculation. This number can be decreased twiceto n = 1 + ν by constructing cheap approximate processors[17, 18, 19]. Taking into account that the evaluation of oneforce gradient is more expensive at least in a factor of ν = 2 than the calculation of one force [13, 14, 15] gives that n ≥ .However, the gradient evaluation may present a difficulty forsystems with long-range (e.g. Coulomb) interactions, wherethe factor ν can be too large [15] because of the necessity tocalculate cumbersome tail (Ewald-summated) contributions.Note also that the processed algorithms of Refs. [17, 18] wereobtained exclusively for pure translational motion and they arenot suitable for rigid-body dynamics. The processing meth-ods introduced in Refs. [16, 19] for solving ordinary differ-ential equations are more general but need an adaptation tobe exploited in the case of rotational motion. In particular,contrary to free translational dynamics, the propagator of freerotational motion cannot be handled at once and requires ad-ditional splitting into analytically integrable parts [15] or in-volving special functions [20].Up to now, no processing schemes were designed and ap-plied to MD simulations of interacting rigid bodies. Therotational motion is much more complicated than transla-tional displacements and thus demands a separate investiga-tion. Moreover, a fundamental theoretical problem on thepossibility to overcome the barrier n = 3 for the fourth-orderintegration still remains open. This overcoming is importantfrom the practical point of view as well, because smaller val-ues of n could noticeably speed up the calculations in view ofthe restricted capabilities of even supercomputers.In the proposing paper we develop the processing formal-ism in the explicit presence of translational and orientationaldegrees of freedom. We show that using a proper transfor-mation of phase coordinates allows to lower the fourth-orderbarrier to the value n = 2 with no gradient evaluations. It isproven also that in a specific case of quasi-fourth-order inte-gration the number of force-gradient evaluations per step canbe reduced to n = 1 at all.The paper is organized as follows. The new processed algo-rithms are consistently derived in Sec. II. Their applicationsto rigid-body MD simulations and comparison with integra-tors known previously are presented in Sec. III. Concludingremarks are highlighted in Sec. IV. II. THEORY
Let us consider a classical system of N interacting rigidpolyatomic molecules. The dynamical state of such a sys-tem in the laboratory frame is determined by the position r i of the center of mass m of the i th molecule, its atti-tude matrix S i as well as the translational p i and angular q i momenta. The equations of motion can be written inthe following compact form d ρ /dt = L ρ ( t ) . Here ρ = { r , p , S , q ; . . . ; r N , p N , S N , q N } ≡ { r , p , S , q } is theset of phase variables, L = N X i =1 (cid:20) p i m · ∂∂ r i + W ( J − S i q i ) S i · ∂∂ S i (1) + f i ( r , S ) · ∂∂ p i + g i ( r , S ) · ∂∂ q i (cid:21) denotes the Liouville operator, f i and g i are the force andtorque, respectively, acting on the molecule due to atomic in-teractions, W ( Ω ) = Z − Ω Y − Ω Z X Ω Y − Ω X is the skewsymmetric matrix related to the principal compo-nents (Ω X , Ω Y , Ω Z ) of the angular velocity Ω = J − Sq with J = diag( J X , J Y , J Z ) being the matrix of moments of iner-tia. If an initial configuration ρ (0) is specified, the uniquesolution to the equations of motion can formally be cast forany time t as ρ ( t ) = [exp( Lh )] k ρ (0) , where h = t/k is thesize of the time step and k denotes the total number of steps.In the standard splitting approach [13, 14, 15], the Li-ouville operator L = A + B is decomposed into its ki-netic A = m − p · ∂/∂ r + W ( Ω ) S · ∂/∂ S and potential B = f ( r , S ) · ∂/∂ p + g ( r , S ) · ∂/∂ q parts (we will omit the sub-script i for the sake of simplicity). Then the one-step timepropagator e Lh can be factorized as e ( A + B ) h + O ( h K +1 ) = Q n +1 µ =1 e Bb µ h e Aa µ h ≡ Φ K ( h ) , where n ≥ and { a µ , b µ } arechosen in such a way to provide the highest possible order K of precision, and O ( h K +1 ) denotes the local error. Forinstance, the second-order ( K = 2 ) Verlet algorithm is ob-tained at n = 1 by e ( A + B ) h + O ( h ) = e B h e Ah e B h ≡ Φ ( h ) .Note that the decomposition constants a µ and b µ should en-ter symmetrically into the factorization to ensure its time re-versibility. This reduces the total number of independentconstants from n + 1) to n + 1 . In turn the symme-try provides automatic cancellation of all even-order termsin O ( h K +1 ) , leading to evenness of K . For even orders K ≥ , the local error function has the form O ( h K +1 ) = c [ A, [ A, B ]] h + c [ B, [ A, B ]] h + O ( h K +3 ) , where [ , ] designates the commutator operation and the coefficients c and c depend on { a µ , b µ } . At K = 2 , the two conditions P µ a µ = P µ b µ = 1 should be satisfied to exclude thezeroth-order term from O ( h K +1 ) . In order to increase the precision to K = 4 we should satisfy the two additional con-ditions c ( { a µ , b µ } ) = c ( { a µ , b µ } ) = 0 . This can be pro-vided by increasing the number n + 1 of independent con-stants at least to the number of the order conditions, i.e, to .We see thus that fourth-order ( K = 4 ) schemes can be con-structed only beginning from n = 3 and this number cannotbe lowered within the standard splitting method. At n = 3 ,the fourth-order ( K = 4 ) factorization can be presented asthe concatenation Φ ( h ) = Φ ( χh )Φ (1 − χ ) h )Φ ( χh ) + O ( h ) of three Verlet signatures, where χ = 1 / (2 − √ .For arbitrary times t , the solution to the equations of mo-tion can be evaluated by consecutively applying k times theone-step splitting propagation Φ K ( h ) . This yields ρ ( t ) =[Φ K ( h )] k ρ (0) + O ( h K ) , where O ( h K ) ∼ k O ( h K +1 ) is theglobal error due to the accumulation of the local one after k = t/h ≫ steps. The action of the exponential operators e Aτ and e Bτ on a phase space point ρ is given analytically by e Aτ (cid:8) r , p , S , q (cid:9) = (cid:8) r + m − p τ, p , Ξ ( q , τ ) S , q (cid:9) , (2) e Bτ (cid:8) r , p , S , q (cid:9) = (cid:8) r , p + f ( r , S ) τ, S , q + g ( r , S ) τ (cid:9) , were the shift of r corresponds to free translational motion (atconstant p ), while the changes in p and q relate to motion ininstantaneous force-torque fields [15]. The matrix Ξ ( q , τ ) ex-actly propagates S over time τ according to the free rotationaldynamics ( q remains constant) d S /dt = W ( J − Sq ) S . Ex-pressions for Ξ ( q , τ ) in terms of efficient routines for ellipticand theta functions are reported in Ref. [20]. Alternatively, Ξ ( q , τ ) can be replaced by its second- or fourth-order coun-terparts Ξ ( τ ) = Ψ X ( τ ) Ψ Y ( τ ) Ψ Z ( τ ) Ψ Y ( τ ) Ψ X ( τ ) and Ξ ( τ ) = Ξ ( χτ ) Ξ ((1 − χ ) τ ) Ξ ( χτ ) , where Ψ ζ ( τ ) =exp[ W (Ω ζ ) τ ] ≡ Θ (Ω ζ , τ ) is the matrix representing rota-tion on angle Ω ζ τ around axis ζ at constant component Ω ζ of Ω = J − Sq (see Eq. (19) of Ref. [15] for Θ (Ω ζ , τ ) ). Notethat each force-torque recalculation in e Bτ requires ∝ N op-erations that is the most time-taking part of the splitting prop-agation, while the costs for handling e Aτ are negligible (pro-portional to N ). The total number of force-torque recalcula-tions per step in Φ K is equal to n .The commutators [ A, [ A, B ]] and [ B, [ A, B ]] which appearin the local error function O ( h K +1 ) can be calculated explic-itly using the expressions for operators A and B . Then, in thecase of the Verlet algorithm ( K = 2 ) we find c = 1 /
12 =2 c and O ( h ) = − (2 m − ˙f · ∂/∂ r − ¨f · ∂/∂ p ) h /
12 + O ( h ) , where at the moment the orientational degrees of free-dom were frozen to simplify notation. Transferring now thecorresponding parts of O ( h ) from e Lh + O ( h ) to the rightunder the exponentials e Ah and e B h one obtains e Lh =e B h e A h e B h + O ( h ) , where A = A + m − ˙f · ∂/∂ r h / and B = B − ¨f · ∂/∂ p h / are the modified counter-parts of A and B . Thus, the order of the Verlet signaturecan increase from K = 2 to K = 4 when the decomposi-tion is performed for the nearby Liouvillian L = A + B = L (1 + m − f · ∂/∂ r h / − ˙f · ∂/∂ p h / , where the equal-ities ˙f = d f /dt = L f and ¨f = L ˙f for the time derivatives of f have been applied. Note however that the nearby exponentials e A τ and e B τ cannot be handled analytically in ρ -space (unlike e Aτ and e Bτ , see Eq. (2)), because of the existence of compli-cated functions ˙f ≡ ˙f ( ρ ) and ¨f ≡ ¨f ( ρ ) which contrary to theforce field f ( r , S ) depend not only on the positions ( r , S ) buton the momenta ( p , q ) as well.The main idea of our approach consists in finding sucha processing transformation ˜ ρ = T ρ from the phase spacepoint ρ to a new set ˜ ρ of variables to make the action ofthe nearby exponentials analytically calculable. Taking intoaccount the explicit structure for the nearby Liouvillian L ,the general form of the desired transformation reads T =( r + αm − f h ) ∂/∂ r + ( p + β ˙f h ) ∂/∂ p + O ( h ) ≡ T α,β ,where α and β are some coefficients which will be definedbelow. It can be verified readily that in the new variables, theequations of motion become d ˜ ρ /dt = ˜ L ˜ ρ , where ˜ L = ˜ A + ˜ B is the corresponding Liouville operator with ˜ A = m − [ ˜p +( α − β ) ˙f ( ˜r ) h ] · ∂/∂ ˜r and ˜ B = [ f ( ˜r − αm − f ( ˜r ) h ) + β ¨f ( ˜r ) h ] · ∂/∂ ˜p . Then for the nearby counterparts of ˜ A and ˜ B one finds ˜ A = ˜ A + m − ˙f ( ˜r ) · ∂/∂ ˜r h / and ˜ B = ˜ B − ¨f ( ˜r ) · ∂/∂ ˜p h / . We see that the terms with ˙f and ¨f can be killedin ˜ A and ˜ B by putting ( α − β ) = − / and β = 1 / , i.e. α = − / . The orientational degrees of freedom can be includedin a similar manner leading to the total processing transfor-mation T α,β = ( r + αm − f h ) ∂/∂ r + ( p + β ˙f h ) ∂/∂ p + Θ (cid:0) J − Sg ( r , S ) , αh (cid:1) S ∂/∂ S +( q + β ˙g h ) ∂/∂ q + O ( h ) andthe nearby operators ˜ A = m − ˜p · ∂/∂ ˜r + W ( J − ˜S˜q ) ˜S · ∂/∂ ˜S and ˜ B = f ( ˜r γ , ˜S γ ) · ∂/∂ ˜p + g ( ˜r γ , ˜S γ ) · ∂/∂ ˜q ≡ ˜ B γ at α = − / , β = 1 / , and γ = − α =1 / . Here Θ (cid:0) J − Sg ( r , S ) , αh (cid:1) = exp[ W ( Q ) αh ] is the matrix representing three-dimensional rotation (i.e. Θ ( Q , τ ) = I cos( Qτ ) + [1 − cos( Qτ )][ W ( Q ) W ( Q ) /Q + I ] + sin( Qτ ) W ( Q ) /Q with I being the unit matrix) aroundvector Q = J − Sg on angle Qαh , and { ˜r γ , ˜S γ } = T γ, { ˜r , ˜S } are the auxiliary position and attitude matrix.From the aforesaid, we have for the one-step propagationin ˜ ρ -space that ˜ ρ ( t + h ) = e ˜ Lh ˜ ρ ( t ) = e ˜ B γ h e ˜ A h e ˜ B γ h ˜ ρ ( t ) + O ( h ) . In ρ -space the solution can be reproduced by applyingthe inverse transformation T − α,β as ρ ( t + h ) = e Lh ρ ( t ) = T − α,β ˜ ρ ( t + h ) . This leads to the resulting propagation of ρ inthe form e Lh = T − α,β e ˜ B γ h e ˜ A h e ˜ B γ h T α,β + O ( h ) , (3)where α = − / , β = 1 / , and γ = 1 / . The operator T α,β transforms a phase space point ρ to the set ˜ ρ = T α,β ρ ≡{ ˜r , ˜p , ˜S , ˜q } of time-step dependent pseudo-variables, where ˜r = r + αm − f ( r , S ) h , ˜p = p + β ˙f ( ρ ) h , (4) ˜S = Θ (cid:0) J − Sg ( r , S ) , αh (cid:1) S , ˜q = q + β ˙g ( ρ ) h . The action of the exponential operators e ˜ A τ and e ˜ B γ τ can begiven analytically as e ˜ A τ ˜ ρ = (cid:8) ˜r + m − ˜p τ, ˜p , Ξ ( ˜q , τ ) ˜S , ˜q (cid:9) , (5) e ˜ B γ τ ˜ ρ = (cid:8) ˜r , ˜p + f ( ˜r γ , ˜S γ ) τ, ˜S , ˜q + g ( ˜r γ , ˜S γ ) τ (cid:9) . Expressions (5) are similar to Eq. (2), since besides the formalreplacement of ρ by ˜ ρ the only difference between ( A, B ) and ( ˜ A , ˜ B ) lies in the modification of the force f ( ˜r γ , ˜S γ ) and torque g ( ˜r γ , ˜S γ ) . Apart from the calculation of theirbasic values f ( ˜r , ˜S ) and g ( ˜r , ˜S ) , the modification requires(for γ = 0 ) one extra force-torque evaluation at the auxil-iary positional ˜r γ = ˜r + γm − f ( ˜r , ˜S ) h and orientational ˜S γ = Θ (cid:0) J − ˜Sg ( ˜r , ˜S ) , γh (cid:1) ˜S coordinates. This increasesthe number of force-torque calculations in e ˜ B γ τ from n = 1 (at γ = 0 ) to n = 2 (at γ = 0 ), but the order of precisionof the processed splitting propagation grows from K = 2 (at α = β = γ = 0 when it reduces to the genuine Verlet signa-ture) to K = 4 (at − α = β = γ = 1 / ).Because of T − α,β T α,β = 1 , the solution to the equa-tions of motion can now be cast for any t as ρ ( t ) = T − α,β [e ˘ B γ h e Ah e ˘ B γ h ] k T α,β ρ (0) + O ( h ) . Then the process-ing transformation T α,β can be performed only once on thevery beginning, while the inverse transformation T − α,β onlyonce at the end of the considered time interval [0 , t ] . In viewof this, the step by step integration can be interpreted as thetime propagation of pseudo-variables ˜ ρ by the kernel splitting e ˘ B γ h e Ah e ˘ B γ h in the transformed phase space. The real phasecoordinates ρ are not involved explicitly into the consecutiveupdating process. They can be reproduced from ˜ ρ whenever itis necessary (for example, when the measurement is desired)using the inverse transformation ρ = T − α,β ˜ ρ . This transfor-mation reads (cf. to Eq. (4)) r = ˜r − αm − f ( ˜r , ˜S ) h , p = ˜p − β ˙f ( ˜ ρ ) h , (6) S = Θ (cid:0) − J − ˜Sg ( ˜r , ˜S ) , αh (cid:1) ˜S , q = ˜q − β ˙g ( ˜ ρ ) h , where the higher-order terms O ( h ) have been neglectedsince they are not accumulated in ρ ( t ) .The next crucial point concerns the evaluation of timederivatives ˙f ( ˜ ρ ) and ˙g ( ˜ ρ ) which arise in Eq. (6). It is ob-vious that their direct evaluation should be obviated sincethis results in complicated gradient terms. Fortunately, thederivatives can be evaluated at a given t in a quite efficientway by the symmetric interpolation { ˙f , ˙g } ( ˜ ρ ) = [ { ˜f , ˜g } ( t + h ) − { ˜f , ˜g } ( t − h )] / (2 h ) + O ( h ) , where { ˜f , ˜g } ( t ± h ) = { f , g } (cid:0) ˜r ( t ± h ) , ˜S ( t ± h ) (cid:1) . Such an interpolation is indeed re-alizable because the pseudo-variables ˜ ρ ( t ± h ) are determinedstep by step in the course of the kernel propagation indepen-dently of ρ ( t ) . Then the real variables ρ ( t ) can be reproducedfrom ˜ ρ ( t ) with a one-step retardation, when the pseudo-phasecoordinates were already propagated to ˜ ρ ( t + h ) . This avoidsthe calculation of extra forces and torques during the inter-polation and involves only those which already were eval-uated within the kernel propagation. The time derivatives ˙f ( ρ ) and ˙g ( ρ ) in Eq. (4) can be evaluated as { ˙f , ˙g } ( ρ ) =[ { f , g } ( h ) − { f , g } ( − h )] /h + O ( h ) , where { f , g } ( ± h ) = { f , g } (cid:0) r ( ± h ) , S ( ± h ) (cid:1) with r ( ± h ) = r (0) ± m − p (0) h and S ( ± h ) = Θ (cid:0) ± J − S (0) q (0) , h (cid:1) S (0) . This involves two ex-tra forces and torques at ± h/ but exclusively on the first stepof the integration when starting from an initial configuration ρ (0) and performing the direct transformation T α,β .We see therefore that the processed splitting (PS) algorithmderived is truly of the fourth order and requires only n = 2 force-torque evaluations per time step. This overcomes thebarrier n = 3 inherent in standard schemes. Moreover, the al-gorithm is time reversible [because the exponential operatorsenter symmetrically into the propagator (Eq. (3))] and phase-area preserving [since simple shifts and rotations (Eq. (5)) donot change the volume]. In addition, the algorithm is explicit(no iterations) and exactly conserves the rigid molecular struc-ture (because Ξ and Θ are rotational matrices). The kernelsplitting can also be made symplectic, because it is based onthe Verlet-like signature which at γ = 0 conserves a nearbyHamiltonian [15, 20, 21]). For a finite γ = 0 , the potential op-erator can be represented by ˜ B γ = ˜ B + γ [ ˜ B , [ ˜ A , ˜ B ]] h / O ( h ) , where [ ˜ B , [ ˜ A , ˜ B ]] = ( ˜ B ε − ˜ B ) /ε + O ( εh ) with ε ≪ . Then the modified force and torque in ˜ B γ can be evaluated as f ( ˜r γ , ˜S γ ) = f ( ˜r , ˜S ) + γ ∆f ( ˜r , ˜S ) and g ( ˜r γ , ˜S γ ) = g ( ˜r , ˜S ) + γ ∆g ( ˜r , ˜S ) , where the secondaryfields are ∆f ( ˜r , ˜S ) = [ f ( ˜r ε , ˜S ε ) − f ( ˜r , ˜S )] /ε + O ( εh ) and ∆g (cid:0) ˜r , ˜S ) = [ g ( ˜r ε , ˜S ε ) − g ( ˜r , ˜S )] /ε + O ( εh ) . The param-eter ε is typically taken to be of order − for double preci-sion arithmetic to minimize the effect of O ( εh ) -terms whileavoiding round-off truncations. The processing transforma-tions (Eqs. (4) and (6)) need not be necessarily symplectic,since their effects are not propagated ( T − α,β T α,β = 1 ).That is very surprising, within the PS method the num-ber n of force-torque recalculations per time step can bereduced to n = 1 at all when a quasi-fourth order is re-quested. Note that the true fourth order means that the de-viations of the generated trajectories ρ ( t ) from their exactcounterparts are equal to O ( h ) ∼ Ch at t ≫ h . In MDsimulations, this strong requirement may not be so needed,because according to the Lyapunov theorem [3] the coeffi-cient C ∼ e λt grows ( λ > ) exponentially with increas-ing t . Then the concept of the quasi-fourth order can bemore useful. It implies that the deviations apply not to in-dividual variables of each particle but rather to a collectivefunction for which C is independent of t . In microcanonicalsimulations such a function should be the total energy E = P Ni =1 ( p i /m + Ω i JΩ i )+ P N ; Mi = j ; a,b ϕ ab ( | r i − r j | , S i , S j ) of the system, where ϕ ab denotes the intermolecular atom-atom potentials, and M is the number of atoms per molecule.Cumbersome analysis shows that E can be conserved withthe fourth-order accuracy at n = 1 by tuning the parame-ters of the method to α = − / , β = 1 / , γ = 0 , and η = 1 / (then ρ ( t ) and other quantities will not be necessar-ily reproduced up to the fourth order). Here we should add anew η -term when transforming (Eq. (6)) angular momentumas q = ˜q − β ˙g ( ˜ ρ ) h + η ˜S + [ JW ( J − ˜S˜q ) J − − W ( J − ˜S˜q ) − W ( ˜S˜q ) J − ] ˜Sg ( ˜r , ˜S ) h , where ˜S + is the transposed matrix(and correspondingly modify Eq. (4)). This adding presentsno difficulty since g ( ˜r , ˜S ) was already calculated during thekernel splitting. For systems without periodic boundary con-ditions, e.g. in celestial mechanics, the total angular momen-tum is often also conserved. It will be kept with the second-order accuracy by the quasi-fourth integrator ( n = 1 ). Thisis in contrast to the genuine fourth-order algorithm ( n = 2 )which produces all quantities to within the O ( h ) precision.Therefore, the former integrator may be less universally appli- cable than the latter one. The PS algorithms will be referredto as PS1 ( n = 1 ) and PS2 ( n = 2 ), respectively.Further improvements are possible by splitting the atom-atom potentials into short- and long-range parts. Then amultiple-time stepping (MTS) technique [22] can be em-ployed, where the expensive long-range (weak) forces aresampled less frequently using larger time steps, while theshort-range (strong) interactions are integrated more accu-rately inside the kernel propagator using smaller steps. TheMTS implementation within the PS method goes beyond thescope of this paper and will be considered elsewhere. III. NUMERICAL RESULTS
We first present the proposed PS method (see Sec. II)in algorithmic form to simplify its implementation in anumerical code. Thus, starting from an initial configu-ration ρ (0) = { r (0) , p (0) , S (0) , q (0) } at t = 0 andcalculating the three forces f (0) and f ( ± h/ as wellas the three torques g (0) and g ( ± h/ at the positions { r (0) , S (0) } and { r ( ± h/ , S ( h/ } , respectively, where r ( ± h/
2) = r (0) ± m − p (0) h/ and S ( ± h/
2) = Θ (cid:0) ± J − S (0) q (0) , h/ (cid:1) S (0) , we make the direct processingtransformation (Eq. (4)) to ˜ ρ (0) = { ˜r (0) , ˜p (0) , ˜S (0) , ˜q (0) } as ˜r (0) = r (0) + αm − f (0) h , ˜S (0) = Θ (cid:0) J − S (0) g (0) , αh (cid:1) S (0) , (7) ˜p (0) = p (0) + β (cid:0) f ( h/ − f ( h/ (cid:1) h, ˜q (0) = q (0) + β (cid:0) g ( h/ − g ( h/ (cid:1) h. Having ˜ ρ (0) , we calculate the two initial forces ˜f (0) and ˜f ε (0) as well as the two initial torques ˜g (0) and ˜g ε (0) at the positions { ˜r (0) , ˜S (0) } and { ˜r ε (0) , ˜S ε (0) } , respec-tively, where ˜r ε (0) = ˜r (0) + εm − ˜f (0) h and ˜S ε (0) = Θ (cid:0) J − ˜S ( ) ˜g (0) , εh (cid:1) ˜S (0) . Note that the direct transforma-tion (Eq. (7)) as well as the evaluation of the initial forces andtorques should be carried out only once at the very beginning( t = 0 ) of the integration.Now we perform the single-step propagations of ˜ ρ fromtime t to t + h according to the kernel splitting (Eq. (3)) as ˜p t + h = ˜p ( t ) + h ˜f ( t ) + γε (cid:0) ˜f ε ( t ) − ˜f ( t ) (cid:1)i h , ˜q t + h = ˜q ( t ) + h ˜g ( t ) + γε (cid:0) ˜g ε ( t ) − ˜g ( t ) (cid:1)i h , ˜r ( t + h ) = ˜r ( t ) + m − ˜p t + h h, (8) ˜S ( t + h ) = Ξ (cid:0) ˜q t + h , h (cid:1) ˜S ( t ) , ˜p ( t + h )= ˜p t + h + h ˜f ( t + h ) + γε (cid:0) ˜f ε ( t + h ) − ˜f ( t + h ) (cid:1)i h , ˜q ( t + h )= ˜q t + h + h ˜g ( t + h ) + γε (cid:0) ˜g ε ( t + h ) − ˜g ( t + h ) (cid:1)i h , where ˜p t + h and ˜q t + h are the intermediate values, and thetwo new forces ˜f ( t + h ) and ˜f ε ( t + h ) as well as the two newtorques ˜g ( t + h ) and ˜g ε ( t + h ) should be calculated at thenew positions { ˜r ( t + h ) , ˜S ( t + h ) } and { ˜r ε ( t + h ) , ˜S ε ( t + h ) } ,respectively, with ˜r ε ( t + h ) = ˜r ( t + h ) + εm − ˜f ( t + h ) h and ˜S ε ( t + h ) = Θ (cid:0) J − ˜S ( t + h ) ˜g ( t + h ) , εh (cid:1) ˜S ( t + h ) before theevaluation of ˜p ( t + h ) and ˜q ( t + h ) . Saving the forces ˜f ( t + h ) and ˜f ε ( t + h ) as well as the torques ˜g ( t + h ) and ˜g ε ( t + h ) ,we repeat Eq. (8) (with formal replacing t by t + h in it) topropagate ˜ ρ from time t + h to t + 2 h . In such a way, step bystep we can recycle Eq. (8) arbitrarily number k ≥ of timesand obtain the value of ˜ ρ ( t ) for any t = kh . Each recycle willrequire the recalculation of only two ( n = 2 ) new forces andtorques.When at least two recycles of Eq. (8) are done already, wewill have the three consecutive values ˜ ρ ( t − h ) , ˜ ρ ( t ) , and ˜ ρ ( t + h ) for some t = kh . The forces ˜f ( t ) and ˜f ( t ± h ) as well as the torques ˜g ( t ) and ˜g ( t ± h ) will also be alreadyknown because of the kernel propagations. Then we can makethe inverse processing transformation (Eq. (6)) of ˜ ρ ( t ) to thegenuine value ρ ( t ) at a current t according to r ( t ) = ˜r ( t ) − αm − ˜f ( t ) h , S ( t ) = Θ (cid:0) − J − ˜S ( t ) ˜g ( t ) , αh (cid:1) ˜S ( t ) , (9) p ( t ) = ˜p ( t ) − β (cid:0) ˜f ( t + h ) − ˜f ( t − h ) (cid:1) h/ , q ( t ) = ˜q ( t ) − β (cid:0) ˜g ( t + h ) − ˜g ( t − h ) (cid:1) h/ , and calculate at this point all necessary observable quantities(such as the total energy, etc.). This completes the PS2 algo-rithm ( n = 2 ), where α = − / , β = 1 / , and γ = 1 / .The PS1 integrator ( n = 1 ) follows at α = − / , β = 1 / , γ = 0 , and η = 1 / (here the evaluation of the modifiedforce ˜f ε and torque ˜g ε should be omitted in Eq. (8) since γ = 0 , while the inclusion of the η -term in Eqs. (7) and (9) istrivial).For testing of the algorithms we applied the TIP4P model( M = 4 ) of water [23] with N = 512 molecules. The MDsimulations were carried in the microcanonical ( N V E ) en-semble at a density of
N/V = and a temperature of292 K. The Ewald summation [24] was exploited to handlelong-range Coulombic atom interactions. The accuracy of thesimulations was measured by calculating the ratio R of thefluctuations of the total energy E to the fluctuations of its po-tential part [15]. The computational costs Υ were estimatedin terms of the number of force-torque evaluations in a giventime interval, taken to be Λ = 1 ps, so that
Υ = n Λ /h . Theequations of motion were solved at several sizes of the timestep ranging from h = 0 . fs to 5 fs. In total k = t/h = 10 steps were used for each algorithm and each step size.The costs Υ versus precisions R of the integration obtainedwithin the two proposed PS algorithms ( K = 4 ) at the endof the simulations are plotted in Fig. 1 by the curves markedas PS1 ( n = 1 ) and PS2 ( n = 2 ), respectively. The resultscorresponding to the Verlet-type (VT) algorithm ( K = 2 and n = 1 ), its optimized (VO) version ( K = 2 and n = 2 ), theForest-Ruth (FR) scheme ( K = 4 and n = 3 ), as well as thegradient-like (GL) algorithm ( K = 4 and n = 3 ) (these inte-grators are described in Ref. [15, 21]) were also included forthe purpose of comparison. It has been established that other FIG. 1: The cost versus relative error for different algorithms in MDsimulations of water. The circles correspond to the time steps (left toright) h = 1 , , , , and 5 fs. The dashed lines represent the mostcharacteristic levels. known rigid-body integrators [8, 9, 10, 11, 12, 25] ( K = 2 and n = 1 ) behave similarly to the VT algorithm. Higher-orderschemes [15, 26] with K ≥ and n ≥ are less efficient inMD simulations because of the large numbers of costly force-torque recalculations. The processed fourth-order algorithmby Blanes and Casas (BC) et al. [16, 19] with K = 4 and n = 1 + ν = 3 (where the kernel and processor are definedaccording to Eqs. (20) and (21) of Ref. [16]) was adapted torigid-body motion and considered too.As can be seen from Fig. 1, with decreasing Υ (rising h )each curve terminates at some point where the simulationsbegin to exhibit a drift in R . This happens around h ∼ fs (larger h can be used within the MTS). At the minimallypossible costs Υ ∼ , the VT integrator can provide onlya crude energy conservation R ∼ . This level of errors istoo large and generally unacceptable in MD simulations. Itshould be reduced at least to R ∼ , arguably the upperlimit of allowable error for which the dynamics can be sim-ulated adequately. The proposed PS1 algorithm just satisfiesthis criteria even at Υ ∼ . On the other hand, the level R ∼ can be achieved by the VT integrator by increas-ing the load to Υ ∼ , i.e. in a factor of 2.75. Thus thePS1 algorithm may spend considerably smaller CPU time ata given precision. The PS2 algorithm is also superior to theVT scheme. For more accurate ( R < ) simulations, therelative efficient of the PS algorithms ( K = 4 ) with respectto the VT scheme ( K = 2 ) rises further (because R ∼ h K )and reaches a factor of 5 at R ∼ . . At the same time, for Υ ∼ the PS2 and PS1 algorithms are able to lower thenumerical errors from the value R ∼ inherent in the VTintegrator to the levels R ∼ . and 0.02%, respectively, i.e.up in 50 times! The VO integrator is clearly inferior to the PSalgorithms, although it can be better than the VT signature.The BC scheme is superior to the VO integrator but worsethan the PS algorithms. The FR scheme leads to the worst ef-ficiency. The GL algorithm can be used only at Υ > , i.e.when a very high accuracy ( R . . ) is required. Then itappears to be more efficient than the PS2 integrator. However, FIG. 2: The fluctuations (a) and deviations (b) of the total energyversus the length of the MD simulations carried out at h = 4 fs usingdifferent algorithms. the PS1 algorithm is the best in the whole Υ -region.Samples of the relative fluctuations R ( t ) and normalizeddeviations δE ( t ) = ( E ( t ) − E (0)) /E (0) of the instantaneoustotal energy E ( t ) are shown in subsets (a) and (b) of Fig. 2, re-spectively, versus the length t/h of the simulations performedat a typical step h = 4 fs using different integrators. We canobserve in Fig. 2(a) that the functions R ( t ) are flat with nodrift on the entire time domain. The PS algorithms thus apartfrom their high efficiency, exhibit also excellent stability prop-erties. As is illustrated in Fig. 2(b) for the PS1 method, thetotal energy E ( t ) continues to keep near its initial value E (0) even after an extremely long period of time with k = 10 steps. The magnitude of the deviations δE ( t ) is quite smalland does not exceed a level of 0.01%, making the energy con-servation almost exact. IV. CONCLUSION
In this paper we have proposed a novel method for theintegration of motion in rigid-body MD simulations. Itcombines standard splitting techniques with special phase-space processing transformations. Comparison with the well-recognized previous schemes has demonstrated that the newmethod allows to significantly improve the efficiency of theintegration with no extra computational costs. The algorithmsobtained are easy in implementation and can readily be incor-porated into existing MD codes. They can also be applied tohybrid Monte-Carlo, MD simulations of simple fluids and toother fields mentioned in the introduction as well as be ex-tended to more complicated systems with flexible molecules.
ACKNOWLEDGMENT
The author acknowledges support by the Fonds zurF¨orderung der Wissenschaftlichen Forschung under theProject No. P18592-TPH. [1] M. P. Allen and D. J. Tildesley,
Computer Simulation of Liquids (Clarendon, Oxford, 1987).[2] D. C. Rapaport,
The Art of Molecular Dynamics Simulation (Cambridge University Press, Cambridge, 1995).[3] D. Frenkel and B. Smit,
Understanding Molecular Simulation:from Algorithms to Applications (Academic Press, New York,1996).[4] S. Essiz and R. D. Coalson, J. Chem. Phys. , 144116 (2006).[5] E. Celledoni and N. S¨afstr¨om, J. Phys. A, Math. Gen. , 5463(2006).[6] S. A. Chin, Phys. Rev. E , 036701 (2007).[7] S. Reich, Fields Inst. Commun. , 181 (1996).[8] A. Kol, B. B. Laird, and B. J. Leimkuhler, J. Chem. Phys. ,2580 (1997).[9] A. Dullweber, B. Leimkuhler, and R. McLachlan, J. Chem.Phys. , 5840 (1997).[10] N. Matubayasi and M. Nakahara, J. Chem. Phys. , 3291(1999).[11] T. F. Miller III, M. Eleftheriou, P. Pattnaik, A. Ndirango, D.Newns, and G. J. Martyna, J. Chem. Phys. , 8649 (2002).[12] H. Kamberaj, R. J. Low, M. P. Neal, J. Chem. Phys. ,224114 (2005).[13] I. P. Omelyan, I. M. Mryglod, and R. Folk, Comput. Phys. Com- mun. , 272 (2003).[14] I. P. Omelyan, Phys. Rev. E , 036703 (2006).[15] I. P. Omelyan, J. Chem. Phys. , 044102 (2007).[16] S. Blanes, F. Casas, and J. Ros, SIAM (Soc. Ind. Appl. Math.)J. Sci. Stat. Comput. , 711 (1999).[17] R. D. Skeel, G. Zhang, and T. Schlick, SIAM (Soc. Ind. Appl.Math.) J. Sci. Stat. Comput. , 203 (1997).[18] M. A. L´opez-Marcos, J. M. Sanz-Serna, and R. D. Skeel, SIAM(Soc. Ind. Appl. Math.) J. Sci. Stat. Comput. , 223 (1997).[19] S. Blanes, F. Casas, and A. Murua, SIAM (Soc. Ind. Appl.Math.) J. Sci. Stat. Comput. , 531 (2004); , 1817 (2006).[20] R. van Zon and J. Schofield, Phys. Rev. E , 056701 (2007).[21] R. van Zon, I. P. Omelyan, and J. Schofield, J. Chem. Phys. ,136102 (2008).[22] M. E. Tuckerman, B. J. Berne, and G. J. Martyna, J. Chem.Phys. , 1990 (1992).[23] W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey,and M. L. Klein, J. Chem. Phys. , 926 (1983).[24] I. P. Omelyan, Comput. Phys. Commun. , 113 (1997).[25] J. P. Ryckaert, G. Ciccotti, and H. J. C. Berendsen, J. Comput.Phys. , 327 (1977).[26] S. Blanes, and F. Casas, J. Phys. A, Math. Gen.39