A Symplectic Integrator for Molecular Dynamics on a Hypersphere

We present a reversible and symplectic algorithm called ROLL, for integrating the equations of motion in molecular dynamics simulations of simple fluids on a hypersphere $\mathcal{S}^d$ of arbitrary dimension $d$. It is derived in the framework of Geometric Algebra and shown to be mathematically equivalent to algorithm RATTLE. An application to molecular dynamics simulation of the one component plasma is briefly discussed

So, for a system of N point particles (positions q α , linear momenta p α = m αqα (as usual dots denote time derivative), α = 1 . . . N ) of E d , we are led to define the angular momentum as [6,10] K is a bivector of the GA algebra G d . In dimension d = 3 the dual of K is the usual 1-vector angular momentum.
It follows fromq α ∧ p α = 0 that a well known theorem of classical mechanics [6,10,12] is generalized to arbitrary dimension : where f α is the total force acting on particle α. Following Goldstein [12] we decompose f α as the sum of an external force f e α and the sum of the pair-wise additive forces due to the other particles β =α f βα . If Newton's weak law of action and reaction holds, i.e. if f βα = −f αβ then we have In the case where there are no external forces, the angular momentum bivector K is a conserved dynamical variable of the motion. This means the conservation of its d 2 scalar components. As well known, this is a consequence of the rotational invariance of the system, as expressed by Noether's theorem [10,12]. As a warm up, let us derive this theorem in the framework of the GA formalism.

B. Rotational invariance and the conservation of the angular momentum bivector
Henceforth in this section, we restrict ourselves to systems made of N classical point particles with identical masses m living in Euclidean space E d . The Lagrangian reads : where the pair potential depends only on the distance |q α − q β | between the particles. Extending Euler-Lagrange equations to an arbitrary dimension d causes no trouble, with the well known result [12] d dt We apply now Noether theorem to the Lagrangian function (5). Clearly L (q α ,q α ) is invariant under any rotation of O(d) of center O. An elemental rotation of space E d is described firstly by its invariant plane, containing the origin, and represented by some unit bivector i (with i 2 = −1), and secondly, by a (scalar) magnitude θ. Bivector Ω = iθ is called the angle of the rotation. The action of this rotation is described by a member of the geometric algebra G d , R, called a rotor (or a generalized complex number) defined as [8][9][10]13] Under the rotation every element A of the algebra is transformed according to the law where the inverse of R, R −1 = R † = exp (+Ω/2). Therefore, for an infinitesimal rotation δΩ, the variation of q α is given by where the medium-dot "·" denotes the inner product of the algebra, cf equation (A3a). idem, the variations of the velocities are given by The infinitesimal rotation yields the following variation for the lagrangian where we took care of the non-associativity of the inner product. However, by taking into account Euler-Lagrange equations as well as the identity (A8) of appendix A, equation (11) can be recast as were we introduced momentum p α = ∂L/∂q α = mq α . Requiring that δL = 0 for any arbitrary infinitesimal bivector δΩ yields, as expected, the conservation of the total angular momentum K :

III. CLASSICAL MECHANICS ON A HYPERSPHERE
A. The bivector angular velocity The hypersphere S d (O, R) of center O and radius R, i.e. the set of points M of E d+1 such that − − → OM 2 = q 2 = R 2 , may be seen as a d-dimensional vector manifold of the euclidean space E d+1 [7]. It is convenient to introduce the unit hypersphere S d = S d (O, R = 1) of unit vectors ξ = q/R. The space E d ξ tangent to S d (O, R) at point q = Rξ plays a central role, here live the particle velocities. It is an euclidean hyperplane of dimension d, spanned by the d orthonormal vectors {e i (ξ)}, i = 1, . . . d. When complemented with vector ξ, this system of vectors forms a local orthonormal basis of full space E d+1 . E d ξ is represented by its unit pseudo-scalar I d (ξ) = e 1 (ξ)e 2 (ξ) . . . e n (ξ) [7,10]. Vector ξ and the d-blade I d (ξ) are related by a duality relation (A6), hence ξ ⋆ = ±I d (ξ), the sign depending of the orientation of the local basis of E d+1 . Therefore, the projection of any blade A of G d+1 onto E d ξ is given by (A7) which gives, in the case of a vector a ∈ E d+1 Let us consider the dynamics of a single point particle q on S d (O, R). Its dynamics is that of a rotator of E d+1 : its position at time t can be obtained from that at time t = 0 by a rotation of O(d + 1) of center O described by some rotor R(t). We thus have [6,10] where R † (t) = R −1 (t) is the reverse of R(t) [6,9,10]. Taking the time derivative of equation (16) one finds thaṫ where the angular velocity Ω(t) is a bivector defined by At this point it must be stressed that equation (16) does not fully determine the rotor R(t) (and thus Ω(t)). Recall that an elemental rotation of O(d + 1) is defined by a 2-plane and an amplitude. If this plane is a subspace of the tangent hyperplane ξ ⋆ then it lets vector ξ unchanged. Therefore there is some freedom and, in order to fix it, we decompose Ω upon the standard basis of the bivectors of the algebra One readily obtains from expansion (19b) that ξ · Ω ⊥ = ξΩ ⊥ 1 = 0 and thus ξ · Ω = ξ · Ω . Similarly, from (19c) one deduces that that ξ ∧ Ω = 0 and thus ξ ∧ Ω = ξ ∧ Ω ⊥ . It is thus possible to chose Ω ⊥ = 0 which fully determines the angular velocity. It amounts to impose the condition Indeed, from Ω ⊥ = 0 andξ = ξ · Ω one obtains all the coefficients of the expansions (19b) and (19c) : Ω ⊥ αβ = 0 and Ω α = −ξ · e α (ξ), which means that Ω = ξ ∧ξ = ξξ. This expression for Ω shows that the two additional conditionṡ ξ ∧ Ω = 0 and ξ ∧Ω = 0 also hold.
It is interesting to consider more closely the case d = 2 for a comparison with the work of Lee et al. [3]. In this case the dual of bivector Ω is an axial vector ω. And one recovers the usual expression for the vector angular velocity ω = Ω ⋆ = ξ ×ξ and the usual relationξ = ω × ξ [3]. The condition ξ ∧ Ω = 0 becomes ξ · ω = 0.

B. Equation of motion on a hypersphere
To simplify the notations we consider the case of a single point particle. For the matters discussed here, the extension to the case of N particles is trivial. This particle, of mass m, position q, is subjected to some external, unspecified force f and constrained to move on the surface of the sphere S d (O, R). Again ξ = q/R denotes its position on the unit hypersphere S d .
Two points of view can be adopted : the motion can be either described as a motion in E d+1 with holonomic constraints or by the theorem (3) derived in section II A and rederived in section III B 2 using a lagrangian formalism in the framework of GA.

Constraints
The mechanical situation is described by the holonomic constraint If f = −∇ q V (q) the Newton's equation of motion can be obtained as the Euler-Lagrange equation associated to the constrained Lagrangian The force due to the constraint is thus −λ∇ q σ = −2λq, it is of course perpendicular to the tangent plane ξ ⋆ . Equation (22) can be simplified as where f ⊥ is the projection of the force f in the tangent hyperplane E d ξ (cf equation (15)) : By passing we note that, taking the scalar product of both sides of (23a) with ξ yields λ = (m/2)ξ 2 where we used the "hidden constraint" ξ ·ξ +ξ 2 = 0 obtained by taking twice the time derivative of the constraint equation (21) valid for all times. λ has been eliminated and the equation of motion reads as

Euler-Lagrange equations
We consider a particle confined on is an arbitrary potential. Euler-Lagrange equations are obtained by minimizing the action with respect to infinitesimal variations δq(t) such that δq(0) = δq(T ) = 0. As noted in reference [3] the expression of δq(t) must agree with the geometric structure of the configuration manifold, which is not a linear vector space but a group. Referring to our discussion of section III A we are led to consider δq = q · δΘ where δΘ(t) is, as the angular velocity Ω(t), an arbitrary infinitesimal time-dependent bivector subject to the condition q ∧ δΘ = 0, i.e. δΘ ⊥ = 0 with the notations of section III A. Then the variation of velocities should read δq =q · δΘ(t) + q · δΘ. At the first order (q + δq) 2 = R 2 and (q + δq) · (q + δq) = 0. The first variation of the lagrangian is therefore where we made use of equation (A8) and of the antisymmetry of the outer product. It follows that the variation of the action is where the usual integration by part on time t and the use of the boundary conditions at t = 0 and t = T were made. Equation (28) can be rewritten in a more transparent way where K = mq ∧q = mR 2 Ω is the angular momentum of the particle. Demanding that δS = 0 for an arbitrary variation δΘ(t) does not automatically ensures that the bivector in brackets A = −q ∧ f ⊥ +K in equation (29) vanishes. To see that, we decompose A = A ⊥ + A as we did in section III A, idem for δΘ. We then note that Then the stationary condition says nothing about A ⊥ but implies that all the A i are zero. Therefore we have where the coefficients c ij are unknown a priori. However since bothK = q ∧q and q ∧ f ⊥ have a null projection onto the plane E d (ξ), it implies that all the c ij of equation (30) vanish. Euler-Lagrange equation is thusK or, in term of the angular velocityΩ We recover the angular momentum theorem of equation (4). Note that taking the inner product of both sides of the above equation with vector q yields back the Newton equation (25) derived in section III B 1.

IV. ROLL : AN INTEGRATOR OF THE EQUATIONS OF MOTION ON A HYPERSPHERE
Let us begin with a short historical digression. The Verlet algorithm was originally devised to integrate numerically the Newton equations of motion of a system of N interacting point particles in the first MD simulations of simple classical fluids or plasmas [15]. The geometry considered in these pioneer works was the d = 3 Euclidean space E 3 , within the usual and convenient periodic boundary conditions [15]. The so called velocity-Verlet algorithm derived later by Swipe et al. is mathematically equivalent to that of Verlet but, numerically is more stable and precise [16]. SHAKE generalizes Verlet integration to systems of point particles with mechanical constraints, for instance molecules and polymers [14], while RATTLE generalizes the velocity-Verlet algorithm for the same class of systems [14]. The analysis of Andersen shows that RATTLE and SHAKE are not numerically equivalent [11]. It appears that SHAKE fails to consider hidden constraints on velocities which are correctly taken into account by RATTLE.
SHAKE or RATTLE can both be used to study the dynamic of particles confined on a hypersphere, but we consider only the superior RATTLE in appendix B. In next section IV A we derive ROLL from Taylor expansions of the equation of motion (25) (c.f. section IV A), and, in section IV B from the exact minimization of a discretized action.
A. Taylor expansions

ROLL-0
We first discretize equation of motion (25) and project it on hyperplane ξ ⋆ n P ξ ⋆ n (ξ n ) = f ⊥ n /(mR) , which gives the projection of the acceleration at time t = nh, where h denotes the time increment. This allows to approximate the positions at discrete times t n±1 with Taylor expansions : Following Verlet and adding both equations gives us [15] It is important to reckon that the knowledge of the projection P ξ ⋆ (ξ n+1 ) implies the full knowledge of the new position ξ n+1 of the particle at time t n+1 . From equation (15) we have In order to satisfy the position constraint ξ 2 n+1 = 1 at time t n+1 , the unknown parameter α must be : For a sufficient small time increment h, ξ n · ξ n+1 > 0 and the sign + should be retained. Therefore Just as the Verlet algorithm, the integrator (34) says nothing about velocities at time t n+1 , but generates a sequence of positions ξ n on the hypersphere. As proposed by Verlet, velocities at time t n can be obtained by subtracting equations (33a) and (33a) yieldingξ Note that these expressions automatically satisfy the hidden constraint on velocities ξ n ·ξ n = 0. To summarize, ROLL-0 is a discrete mapping (ξ n−1 , ξ n ) → ξ n+1 which can be written, for a sufficiently small time increment h : it satisfies both constraints ξ 2 n+1 = 1 and ξ n ·ξ n = 0.

ROLL-2
We skip ROLL-1 which will be derived in next section IV B and we present now ROLL-2. We start from a Taylor expansion about time t n rewritten with the help of bivector angular velocity Ω n We project the equation on the tangent plane E n ξ obtaining As explained in section IV A 1 P ξ ⋆ n (ξ n+1 ) determines ξ n+1 . We need an expression for Ω n+1 . In this aim we discretize equation (31) and adopt the velocity-Verlet recipe : To summarize, ROLL-2 is a discrete mapping (ξ n , Ω n ) → (ξ n+1 , Ω n+1 ) which can be written, for a sufficiently small time increment h : B. The discrete action: ROLL as a variational integrator

A variational integrator
We extend the analyses of Lee et al. concerning variational integrators of the discrete action on two-spheres [3] to the case of hyperspheres of arbitrary dimension. Viewing Verlet integrator as an exact variational principle of a discrete action rather than the approximate Taylor expansion of exact equations of motion is an idea due to Marsden and collaborators, see e.g. reference [17]. Here it yields an unification of ROLL-0 and ROLL-2, lets us discover ROLL-1, and shows the mathematical equivalence of all these versions of ROLL.
As in section IV A we consider only the case of one particle in the potential V (ξ) since the extension to N particles is trivial. The discrete lagrangian of the particle can be defined as [17] where V n ≡ V (ξ n ) is the total potential energy of the system in the configuration ξ n and at discrete time t n = nh. The discrete action between times t 0 = 0 and T = n max h is defined as We now compute the variation δS induced by some infinitesimal variation of the positions δξ n . As in the continuous case discussed section III B 2, we must be careful in the definition of δξ n and reckon that the configurational space is not a vector space but a group. Therefore, at each time t n , we consider an infinitesimal variation in the tangent n-plane ξ ⋆ n : δξ n = ξ n · δΘ n , where δΘ n is some bivector of the form δΘ n ≡ δΘ n = 1≤i≤d ǫ n i e i (ξ n ) ∧ ξ n where the scalars ǫ n i are d arbitrary infinitesimal variations. As discussed in section III B 2 this expansion of δΘ n ensures that the projection of bivector δΘ n in the tangent plane E n ξn is zero, in other words that ξ n ∧ δΘ n = 0. Discrete Hamilton's principle states that δS = 0 for arbitrary ǫ n i submitted to boundary conditions ǫ 0 i = ǫ nnmax−1 i = 0. The variation of the lagrangian is where we have used Marsden's notations D ξn = ∂/∂ξ n and D ξn+1 = ∂/∂ξ n+1 . It follows from (46) and a discrete integration by parts, (i.e. a relabeling of the second sum), that the variation of the action is where we made use of equation (A8). The condition δS = 0 for an arbitrary δΘ with the condition ξ n ∧ δΘ n = 0 yields, as discussed at length in section III B 2 to the equations After taking the inner product of both sides of (49) we obtain which is nothing , when complemented with equation (37), but the integrator ROLL-0 (cf equations (39)).

linear and angular momenta
The linear momentum p n is defined in the tangent plane E d (ξ n ) as follows [3,17] : This is satisfied for any δξ n perpendicular to ξ n , i.e. of the form δξ n = ξ n · δΘ n with ξ n ∧ δΘ n = 0. Therefore we have This equation must be satisfied for any δΘ n ≡ δΘ n and, following the same technique than in section IV B 1 we find after some algebra Similarly one defines the momentum p n+1 as and one obtains, with computations in the same vein as those used to compute p n , the results From equations (53) and (53) one obtains easily two versions of ROLL integrator.
• ROLL-1 From equations (53a) and (55a) one obtains a discrete mapping (ξ n , p n ) → (ξ n+1 , p n+1 ) which can be written, for a sufficiently small time increment h : Equations 56 constitute ROLL-1 integrator. After definingξ n by the identity p n ≡ mRξ n , it turns out that there are identical to the equations (B10), constitutive of RATTLE algorithm, which are derived in appendix B.

Properties of ROLL
It follows from the analysis of section IV B 1 that ROLL-1 and ROLL-2 are two mathematically equivalent integrators. ROLL-1 is also mathematically equivalent to RATTLE, therefore all these algorithms share the same properties and are thus reversible and symplectic, as shown for RATTLE in reference [11,18]. In the case d = 2 it is a simple exercice to show that ROLL-2 coincides with the integrators of references [3,4].

V. NUMERICAL APPLICATION
The main application of ROLL should be for MD simulations of simple fluids. Actually, the results of preceding sections can easily be extended to a system of N particles. In case of pairwise additive potentials, for instance, the perpendicular force f ⊥ α acting on particle "α" can be written where it had been recognized that the pair potential v(ψ αβ ) depends only on the geodesic distance ψ αβ = cos −1 (ξ α · ξ β ) on the unit hypersphere. In equation (57) t αβ (ξ α ) denotes the tangent to the geodesic (ξ α , ξ β ) at point ξ α and orientated from ξ α to ξ β [19,20].
We applied the ROLL algorithm to MD simulations of the three dimensional version of the OCP. The OCP is a model which consists of identical point charges q of mass m, immersed in a uniform neutralizing background of charge density ρ = −nq where n = N/V is the numerical density of particles (V = 2π 2 R 3 volume of the system) [21,22]. In the thermodynamic limit, its properties depend on the sole dimensionless coupling parameter Γ = βq 2 /a (β = 1/kT , T temperature, k Boltzman constant, and a ion-sphere radius, defined in E 3 , by a = (4πn/3) −1/3 ).
The expression for the configurational energy the OCP of in S 3 (O, R) was derived in [23] and will not be repeated here. We only mention that the pair potential v(ψ) entering the configurational energy is essentially the electrostatic potential created by a point charge embedded in its neutralizing background. It is obtained by solving Poisson's Equation analytically in S 3 (O, R) [19,20].
In reference [23] MC simulations in the canonical ensemble, performed on the hypersphere S 3 (O, R), have produced very accurate results for the thermodynamic limit of the internal energy of the model with a relative precision p ∼ 10 −5 in the range 1 ≤ Γ ≤ 190. In these MC simulations the parameter Γ is fixed while, in the MD simulations, the total energy is fixed and the potential and kinetic energy (and therefore parameter Γ) must both be computed.
In our simulations we have chosen a system of units such that q = 1, m = 1, and a = 1. With this choice the unit of time, for instance is t 0 = (ma 3 /q 2 ) 1/2 = 1. We computed the potential energy V and the kinetic energy KIN = (m/2)R 2 αξ 2 α as well as their sum E tot = KIN + V . From the mean value of the kinetic energy we obtained the parameter Γ = 3N/(2 < KIN >). We present here only an example to illustrate the validity of algorithm ROLL-2.  Instantaneous and cumulated potential (left) and total (right) energies per particle for a sub-run of 2 × 10 5 steps. From the mean kinetic energy one has Γ ∼ 30.
A more extensive numerical study will be published elsewhere. A sample of N = 2000 particles was prepared in the state Γ ∼ 30. . We chose for the time increment h = 0.01 (in reduced units). The initial state was prepared so that each of the 6 components of the total angular momentum bivector per particle was ∼ 10 −16 . ROLL-2 ensures a perfect conservation of the angular velocity. With this respect ROLL-1 is slightly inferior because the angular velocities must be computed from the positions and the momenta yielding a slight loss in precision. Figure 1 shows data for the potential and total energies for a sub-run of 2 × 10 5 time steps. Since ROLL-2 is symplectic the conservation of the total energy is excellent. For a complete MD simulation involving 7 × 10 6 time step, we obtained • Γ = 30.002 (4). The reported statistical error, which corresponds to two standard deviations, was computed by dividing the run in n b = 1400 blocs of 5000 time steps and was obtained by a standard block analysis [24] • < V > /N = −0.848137(5), i.e. a relative error of 6 × 10 −6 . This value should be compared with the thermodynamic value < V > /N = −0.8481197(42) obtained in reference [23] by a finite size studies of MC simulations of the OCP at Γ = 30.
• E/N = −0.7981404542(4), i.e. a relative error of 5 × 10 −10 . This error is of course strongly dependent on the value of the time increment h = 0.01. Decreasing h to improve on the conservation of the total energy is not a good idea since then the phase space is not properly sampled. A more appealing analysis is probably to keep track of the extreme values of E/N during the run. We found −0.79814057 < E/N < −0.79814033 exemplifying the absence of drift and the excellent energy conservation for a quite long run of 7 × 10 6 time steps.

VI. CONCLUSION
We have presented in this article the symplectic integrator ROLL for MD simulations of simple fluids or plasmas on the hypersphere. It extends to hypersphere S d of arbitrary dimensions, notably d = 3, previous attempts made for 2-spheres [3,4]. The use of GA theory plays an important role in this extension because, otherwise, it would be impossible to introduce angular velocities without resorting to awkward matricial representations. In the special case d = 3, an alternative to GA could be the use of quaternions. However, this would be only but a peculiar case of the more general approach given here, which is valid even for d > 3.
Preliminary tests in the case of the OCP demonstrated the validity of ROLL and applications to extensive numerical experiments on OCP and related models are under way.
Since GA is also a very powerful approach for the classical mechanics of rigid bodies [10], the extension of ROLL to molecular fluids is clearly possible. The rotations of bodies in the tangent plane E d ξ could be described by the component Ω ⊥ of the angular velocity set to zero for point particles. Applications to polar fluids and ionic solutions is appealing.