{'fraction_non_alphanumeric': 0.037422201213267156, 'fraction_numerical': 0.017332387930355315, 'mean_word_length': 4.386590282198175, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
NVU dynamics. III. Simulating molecules at constant potential energy

This is the final paper in a series that introduces geodesic molecular dynamics at constant potential energy. This dynamics is entitled NVU dynamics in analogy to standard energy-conserving Newtonian NVE dynamics. In the first two papers [Ingebrigtsen et al., J. Chem. Phys. 135, 104101 (2011); ibid, 104102 (2011)], a numerical algorithm for simulating geodesic motion of atomic systems was developed and tested against standard algorithms. The conclusion was that the NVU algorithm has the same desirable properties as the Verlet algorithm for Newtonian NVE dynamics, i.e., it is time-reversible and symplectic. Additionally, it was concluded that NVU dynamics becomes equivalent to NVE dynamics in the thermodynamic limit. In this paper, the NVU algorithm for atomic systems is extended to be able to simulate geodesic motion of molecules at constant potential energy. We derive an algorithm for simulating rigid bonds and test this algorithm on three different systems: an asymmetric dumbbell model, Lewis-Wahnstrom OTP, and rigid SPC/E water. The rigid bonds introduce additional constraints beyond that of constant potential energy for atomic systems. The rigid-bond NVU algorithm conserves potential energy, bond lengths, and step length for indefinitely long runs. The quantities probed in simulations give results identical to those of Nose-Hoover NVT dynamics. Since Nose-Hoover NVT dynamics is known to give results equivalent to those of NVE dynamics, the latter results show that NVU dynamics becomes equivalent to NVE dynamics in the thermodynamic limit also for molecular systems.


I. INTRODUCTION
In two recent papers 1,2 (henceforth: Papers I and II) molecular dynamics at constant potential energy was introduced, tested, and compared to well-known molecular dynamics algorithms. This new molecular dynamics is entitled NVU dynamics in analogy to standard energy-conserving Newtonian NVE dynamics. The conclusion was that NVU dynamics is a fully valid molecular dynamics, which for sufficiently large systems can be used interchangeably with NVE dynamics for calculating most quantities of interest. NVU dynamics is not faster than standard NVE or NVT dynamics, but introduces a new way of thinking about molecular dynamics. Molecular dynamics at constant potential energy was previously considered by Cotterill and co-workers 3-6 , by Scala et al. 7 , and most recently by Stratt and co-workers [8][9][10][11] , who actually allowed also lower potential energy values. Our motivation for studying NVU dynamics derived from recent work on strongly correlating liquids and their isomorphs [12][13][14][15][16][17][18][19] (see the Introduction of Paper I).
NVU dynamics is defined by geodesic motion on the constant-potential-energy hypersurface Ω defined by Here R ≡ {r (1) , ..., r (N ) } in which r (k) is the position vector of the k'te particle (we follow here the notation of the Appendix of Paper II), and U is the potential-energy function of an N -particle classical system. A geodesic on Ω is a curve that satisfies the condition of stationary length for fixed endpoints R A and R B , i.e., where dl is the line element of the metric. The shortest path between any two points is a geodesic. On a sphere geodesics are great circles, the "straightest lines" of the surface. Traversing a geodesic at constant velocity thus corresponds to a generalization of Newton's first law to a curved space (the surface itself). In Paper I the NVU algorithm was developed via a discretization of Eq. (2), subsequently carrying out the variation. This technique, which is known as variational integration [20][21][22][23] , resulted in a "basic" NVU algorithm that is similar to the well-known Verlet algorithm R i+1 = 2R i − R i−1 + (∆t) 2 F i /m for Newtonian (NVE ) dynamics (m is the particle mass which is assumed identical in this section, and F i ≡ −∇ R i U is the 3N -dimensional force vector); the index i refers to step i of the integration sequence. In the Verlet algorithm ∆t is a fixed time step length. In comparison, the basic NVU algorithm is given by (Paper I) If the number of particles N increases, the relative variation of the term −2F i · (R i − R i−1 )/F 2 i decreases, and this is why equivalence with Newtonian NVE dynamics is established in the thermodynamic limit. This equivalence should be understood in the sense that the relative deviations between, for instance, NVE and NVU auto-correlation functions go to zero as N → ∞.
Paper I additionally developed a "stabilized" version of the basic NVU algorithm to prevent accumulation of numerical errors. This version of the algorithm is given by (defining the position changes where l 0 is the step length and All simulations in Papers I and II were performed with the stabilized algorithm. The basic algorithm was used, however, for theoretical considerations.
In this article we extend the stabilized NVU algorithm to deal with simulations of molecular systems. Molecular systems are simulated by introducing rigid and/or flexible bonds between the atoms in the modelling. Flexible bonds introduce merely an additional contribution to U , for instance, harmonic spring potentials. The basic NVU algorithm conserves the total potential energy and can readily simulate flexible bonds. The focus in this paper is thus on implementing rigid bonds in the framework of NVU dynamics.
Section II considers NVU dynamics with rigid bonds. Introducing rigid bonds in the simulations leads to Lagrangian multipliers in addition to those introduced in order to keep the potential energy constant (Paper I). Section II is fairly technical and easiest to read after reading Paper I. Section III gives simulation and model details. Section IV tests the rigid-bond NVU algorithm, and Sec. V investigates the NVU sampling properties by comparing the NVU results to Nosé-Hoover NVT results 24,25 on three different systems: the asymmetric dumbbell model 26 , Lewis-Wahnström OTP 27 , and rigid SPC/E water 28 . Nosé-Hoover NVT dynamics is known to give results equivalent to NVE dynamics in the thermodynamic limit 29 , and we refer to these dynamics interchangeably in the forthcoming sections. Finally, Sec. VI concludes.

II. RIGID-BOND N V U ALGORITHM
The rigid bonds 30,31 introduce constraints among the particle coordinates of the system. Each constraint α = 1, ..., G is of the form it expresses that the distance between particles k α and l α is a constant, C α . In Papers I and II the integral of Eq.
(2) was merely restricted to the constant-potential-energy hypersurface Ω. Each rigid bond constraint introduces a function σ α to be kept constant, and thus the integral of Eq. (2) is now further restricted to the sub-manifold ω of Ω where the bond constraints are satisfied, If the bond constraints are independent, as assumed throughout the paper, ω is a (3N -G -1)-dimensional compact Riemannian manifold. The variational principle defining NVU dynamics with rigid bonds is given by Most of Papers I and II dealt with the case of identical particle masses, but we wish here to develop a completely general molecular NVU algorithm. The line element dl is defined by wherem k = m k / m is the "reduced" mass of particle k. Equation (10) is not the standard Euclidean line element, but a mass-weighted line element that goes back to Hertz 32,33 . We shall refer to this metric as the "Hertzian" metric. The point of this particular metric is that it ensures equivalence between NVU and NVE dynamics for systems of atoms and molecules of varying mass. In appendix A we derive the variable-mass atomic NVU algorithm applying the Hertzian metric (correcting also a typo of the Appendix of Paper II).
Applying the variational integration technique to Eq. (9) gives In Eq. (11) the path is divided into a number of discrete points and one Lagrangian multiplier Λ αi is introduced for each constraint α at every point i. Following standard notation for constraint molecular dynamics 30,31 , the Lagrangian multipliers of the bond constraints are chosen with a positive sign. As in Papers I and II we now make the Ansatz of constant step length l 0 , i.e., km k r Carrying out the variation of Eq. (11) using Eq. (12) leads to (compare the derivation in Paper I) where f U is the force on particle k at step i. This equation constitutes the NVU algorithm with rigid bonds. It has a close resemblance to the Lagrangian equations of motion with holonomic constraints 31 , i.e., rigidbond NVE dynamics 30 . Equation (13) contains G + 1 Lagrangian multipliers for each integration step, which must be determined to complete the algorithm.

A. Determining the NVU Lagrangian multipliers
This section shows how to calculate the Lagrangian multipliers. Since the algorithm is to be implemented on a computer (with finite-precision), we shall proceed directly to a "stabilized" algorithm conserving for indefinitely long runs potential energy, bond lengths, and step length (in 3N -dimensions). The resulting algorithm reduces to the stabilized atomic NVU algorithm of Eqs. (4)- (6) in the case of no bonds constraints.
Some notation used in the following derivation is now introduced (the nomenclature of text is summarized in Table  I).

3-dimensional vectors r
Displacement of the position of particle k between step i and i + 1. (δ Displacement of the positions of particles kα and lα at step i.
Displacement of the velocities of particles kα and lα Sum of displacements of positions and velocities of particles kα and lα at, respectively, step i and Displacement of the forces on particles kα and lα at step i divided by their reduced particle mass.
Displacement of the constraint forces on particles kα and lα at step i divided by their reduced particle mass.
Position of all particles at step i.
Displacement of the positions between step i and i + 1.
Force on all particles at step i.
Force on all particles at step i divided by the reduced particle mass.
Constraint force on all particles at step i divided by the reduced particle mass. Defining δ In analogy to rigid-bond NVE dynamics we call g (k) i the "constraint force" on particle k at step i. Introducing the notationF i ≡ {f The Lagrangian multipliers are calculated by combining a result derived in Paper I with the method applied in the SHAKE algorithm 30 for rigid bonds in NVE dynamics 30,35,36 . The SHAKE algorithm calculates the Lagrangian multipliers from the equations (r α i+1 ) 2 = C 2 α . In doing so, the target value of the constraints C α appears explicitly in the algorithm, making the bond lengths insensitive to numerical error. The expression for r α i+1 is supplied by the integration algorithm containing herein the Lagrangian multipliers. In our case, this gives G equations with G + 1 unknowns. The missing equation is supplied by an expression derived in Paper I, namely that to third order in the step length. In the discrete sequence of points U i+1 is set equal to U 0 (the constant defining Ω), making the constraint of constant potential energy also insensitive to numerical error. We thus have the following G + 1 equations for calculating the Lagrangian multipliers By Eqs. (16) and (17); /m lα , since by Eqs. (14) and (15); The above coupled quadratic equations for the Lagrangian multipliers are now solved following the produce of the MILC-SHAKE algorithm 37 , which starts by neglecting the second order terms in the Lagrangian multipliers and solving the resulting linear equations. Afterwards, the second order terms are taken into account in an iterative manner -the details of which are described below. For each integration step i, the linearized equations are given by where .., Λ Gi }, and b i a G + 1 column vector. We start by calculating explicitly the first few elements of the matrix A i . A 11 consists merely of the factor infront of λ i in Eq. (20), i.e., A 11 = −l 0Fi · F i . The second element A 12 appears after expansion of the dot product The last equation follows as the Lagrangian multipliers appear in pairs, differing only by the sign from ∇ r (kα ) i σ α and the term In the second row of A i , the short-hand notation , the factor infront of λ i after squaring of the parentheses. The next element A 22 appears after expanding s 1 In this sum, we identify the factor in front of Λ 1i , giving , and similarly for the remaining elements of the second row.
Altogether, the elements of A i are thus given by The column vector b i consists of all zeroth-order terms in Eqs. (20) and (21) b Turning now to the iteration procedure, the second-order terms in the Lagrangian multipliers (Eq. (21)) are taken into account by iterating the right-hand side of Eq. (22) via the scheme (α = 1, ..., G) b j+1 The superscript j refers here to iteration j, and (r α i+1 ) 2 j are the positions associated with iteration j. The element b 0 is not updated as it derives from the constraint of constant potential energy. For each iteration j the term C 2 α − (r α i+1 ) 2 j is expected to become smaller as the bonds are satisfied better and better, and convergence was achieved within a few iterations 37 .
For each integration step i, the algorithm for determining the NVU Lagrangian multipliers thus proceeds as follows 1. The Lagrangian multipliers of iteration j, (λ i ) j , are calculated from Eq. (22).
How is constant step length l 0 ensured numerically? Generalizing the approach of Paper I we introduce a normalizing factor such that where i .
The normalizing factor is close to unity 1 and ensures trivially km k (δ (k) i+1/2 ) 2 = l 2 0 , i.e., that the step length is conserved. The algorithm is now absolutely stable, conserving potential energy, bond lengths, and step length for indefinitely long runs. The stability of the NVU algorithm is tested numerically in Sec. IV.

B. Alternative determination of the NVU Lagrangian multipliers
The previous section followed the traditional way of calculating the Lagrangian multipliers. The NVU Lagrangian multipliers may also be calculated by Taylor expanding the constraints σ α in analogy to the method sketched above for the potential energy. In this way, the constraints of constant potential energy and constant bond lengths are treated on equal footing. The set of equations to be solved is the following (recall that R i+1 −R i−1 = 2∆ i−1/2 +l 0 λ iFi +l 0Gi ) The determination of the Lagrangian multipliers is linear and thus no iterations are needed. The bond constraints σ α are obeyed to the same order O(l 3 0 ) as the constraint of constant potential energy. The sampling properties of this novel, alternative determination method is tested briefly in Sec. V. It appears to be a promising new way of determining the Lagrangian multipliers in connection with rigid bonds, which might also be useful for standard bond-constraint NVE or NVT simulations.

III. SIMULATION DETAILS AND MODEL SYSTEMS
We investigated three systems: The asymmetric dumbbell model, the Lewis-Wahnström OTP model, and rigid SPC/E water. For all simulated pair potentials the shifted-force truncation scheme was applied at a cut-off radius r c . If the pair potential is v(r) and the pair force is f (r) = −v ′ (r), the shifted force is given by 34,38 This corresponds to using the following pair potential below r c : v SF (r) = v(r) − v ′ (r c )(r − r c ) − v(r c ). All simulations were performed with the NVT and NVU algorithms. Recall that NVE and NVT dynamics give equivalent results 29 ; for this reason no simulations are presented for NVE dynamics. The RUMD code 39 was used for molecular dynamics simulations (an optimized open-source GPU code). The NVT ensemble is generated via the Nosé-Hoover algorithm 24,25,40 , and the bonds held fixed using the time-reversible constraint algorithm of Refs. 35and36 . The NVU algorithm is described in Sec. II. The starting files for NVU dynamics were taken from an equilibrated NVT simulation. The positions and velocities of the NVT configuration do not correspond perfectly to motion on ω, since the potential energy and step length are not those of U 0 and l 0 , respectively. As all three constraints are to be satisfied simultaneously, this results in numerical problems when starting the simulation from the particular NVT configuration. A more gentle procedure is thus applied, where the atomic NVU algorithm is used for a couple of integration steps to ensure the values of U 0 and l 0 . Afterwards, the rigid-bond NVU algorithm is used.
A. NVU iteration procedure The quadratic equations (Eq. (25)) were iterated with a fixed number of iterations (between 3 and 5). The linear systems were solved utilizing CUSP 41 , a library for solving systems of linear equations on the GPU. More specifically, the stabilized biconjugate gradient algorithm with a Jacobi preconditioner 42 was used with the initial value λ i = 0. The relative tolerance τ for the asymmetric dumbbell and Lewis-Wahnström OTP models was chosen as τ = 10 −7 and for rigid SPC/E water as τ = 3 · 10 −7 . A larger tolerance was chosen for rigid SPC/E water due to convergence issues in connection with shifted-force Coulomb interactions (see below).
The maximum number of allowed iterations was 50. A restart scheme was applied when the solver did not converge within the chosen tolerance. In this case the solver (and quadratic iteration) was restarted from the partially estimated "solution" adding 2 · 10 −7 to the tolerance. It should be noted that the stabilized biconjugate gradient algorithm may get "trapped", resulting in a break-down of the CUSP linear solver. If this happens, it is detected by our program, and the solver and quadratic iteration are restarted, with a smaller number (10) of maximum allowed iterations for the solver.

B. The asymmetric dumbbell
The asymmetric dumbbell model 26 consists of a large (A) and a small (B) Lennard-Jones (LJ) particle, rigidly bonded with bond distance of r AB = 0.29/0.4963 (here and henceforth units are given in LJ units referring to the A particle such that σ AA = 1, ǫ AA = 1, and m A = 1). The asymmetric dumbbell model has σ BB = 0.3910/0.4963, ǫ BB = 0.66944/5.726, and m B = 15.035/77.106. The AB interaction between different molecules is determined by the Lorentz-Berthelot mixing rule 34 . n = 500 molecules (here and henceforth n denotes the number of molecules and N the number of atoms) were used in the simulations with a pair-potential cut-off of r c = 2.5. The step length l 0 was fixed in the range 0.125-0.138 depending on the state point.
Simulations were also performed where the rigid bonds were replaced with stiff harmonic springs. The spring constant was k = 3000, while all other model parameters remained unchanged.

C. Lewis-Wahnström OTP
The Lewis-Wahnström OTP model 27 consists of three identical LJ particles rigidly bonded in an isosceles triangle with sides of r AA = 1 and top angle of 75 • . All parameters (including the masses) are unity for the OTP model. n = 320 molecules were simulated and a pair-potential cut-off of r c = 2.5 was used. The step length was 0.100.

D. SPC/E Water
The SPC/E water model 28 is an isosceles triangle with sides r OH = 1/3.166 and top angle 109.47 • . The OO intermolecular interactions are given by the LJ pair potential (ǫ OO = 1, σ OO = 1, and m O = 15.9994/1.00794). The three particles are charged with q O = −22.0 and q H = |q O |/2. n = 2000 molecules were simulated and a pair-potential cut-off of r c = 6.28 for both LJ and Coulomb interactions was applied 43,44 . The step length was fixed in the range 0.06-0.07 depending on the state point. For this system the numerical stability is surprisingly sensitive to the cut-off used in the Coulomb interactions, but a larger shifted-force cut-off improves this behavior 44 .

IV. TESTING THE STABILITY OF THE RIGID-BOND N V U ALGORITHM
This section tests the conservation properties of the rigid-bond NVU algorithm. Table II shows the potential energy, the deviation of bond lengths, and step length as functions of integration step number for Lewis-Wahnström OTP at ρ = 0.329 and T = 0.700. It is clear that these quantities are conserved by the algorithm and that no drift occurs. The step length is conserved to the highest accuracy since it is not prone to numerical error in determining the Lagrangian multipliers.  In NVU dynamics there is, as such, no notation of time; a geodesic on the manifold can be traversed with any velocity. Comparing the NVU algorithm of Eq. (13) to the rigid-bond Verlet algorithm 30 r

Integration steps
i ], we can define the term l 0 λ i m as a varying "time step" length of the NVU algorithm (see also Paper II), i.e., The integration steps of the NVU algorithm are thus henceforth referred to as "time steps". The average of Eq. (32) is used in Sec. V when comparing to NVT dynamics. As was the case for the atomic NVU algorithm (Paper I), l 0 λ i m is Gaussian distributed for large systems and its relative variation decreases as the number of particles increases. It thus becomes a better and better approximation to treat this term as constant, implying equivalent sampling properties of NVU and NVE dynamics also when rigid bonds are included in the simulations. The NVU algorithm is now compared to NVT dynamics for the three different models. First, we consider the asymmetric dumbbell model 26 , both rigid and flexible. Afterwards, the Lewis-Wahnström OTP model 27 , and finally rigid SPC/E water 28 .
A. The asymmetric dumbbell model In Figs. 2(a) and (b) are shown, respectively, the molecular center-of-mass (CM) radial distribution functions and the CM incoherent intermediate scattering functions for the rigid asymmetric dumbbell model 26 for different temperatures at ρ = 0.932. The black circles and curves give NVT simulation results while the red crosses give the NVU simulation results. The two radial distribution functions in Fig. 2(a) agree very well, and this is also the case for the dynamics in Fig. 2(b). For reference, we also simulated (Fig. 3) the corresponding quantities for the flexible-bond asymmetric dumbbell model at the state points of Fig. 2. Again, there is a very good agreement between NVU and NVT dynamics. For comparison, we also show in Fig. 5 the corresponding particle quantities for the OTP model.

C. SPC/E Water
Finally, we consider in Fig. 6 the same quantities as above for the rigid SPC/E water model 28 . Again, full equivalence between NVU and NVT dynamics is found. The linear algorithm for determining the Lagrangian multipliers presented in Sec. II B (Eqs. (29) and (30)) is tested in Fig. 7 by probing the molecular CM radial distribution functions. NVU and NVT dynamics also here give identical results. Comparison of structure in NVU and NVT simulations of rigid SPC/E water at ρ = 1.000 and T = 3.800 applying the linear method to determine the Lagrangian multipliers (Eqs. (29) and (30)). The bond lengths are here conserved to order 10 −6 in the standard deviation of the bonds (using single-precision).
We conclude from the presented results that for sufficiently large molecular systems with flexible and/or rigid bonds, NVU dynamics is equivalent to Nosé-Hoover NVT dynamics (and, by implication, Newtonian NVE dynamics).

VI. SUMMARY
NVU dynamics is molecular dynamics at constant potential energy realized by tracing out a geodesic on the constant-potential-energy hypersurface Ω (Eq. (1)). In Papers I and II 1,2 , a "basic" and a "stabilized" atomic NVU algorithm for simulating geodesics on Ω was developed. The basic NVU algorithm has excellent stability; it is timereversible and symplectic, and the stabilized algorithm was developed only to prevent accumulation of numerical error. It was found that atomic NVU dynamics becomes equivalent to atomic NVE dynamics in the thermodynamic limit.
In this paper the stabilized NVU algorithm has been extended to simulate molecules at constant potential energy. Molecules are simulated by introducing rigid and/or flexible bonds in the models. The atomic NVU algorithm keeps the potential energy constant and can thus right away simulate flexible bonds. The focus here was on incorporating rigid bonds in the framework of NVU dynamics, which leads to the introduction of additional Lagrangian multipliers beyond those of the constraint of constant potential energy. This is completely analogous to the approach for simulating rigid bonds in standard Newtonian NVE dynamics 30,35,36 . In the NVU algorithm, a set of coupled quadratic equations was constructed for calculating the Lagrangian multipliers and solved in an iterative manner as a linear system, a procedure developed for rigid-bond NVE dynamics in the MILC-SHAKE algorithm 37 . In addition, a set of linear equations was presented for calculating the Lagrangian multipliers, and appears to be a promising new way of simulating rigid bonds.
The rigid-bond NVU algorithm reduces to the atomic NVU algorithm when there are no rigid bonds. The algorithm was tested on three different model systems: the asymmetric dumbbell model, Lewis-Wahnström OTP, and rigid SPC/E water. The probed quantities in the simulation gave identical results to those of Nosé-Hoover NVT dynamics. We conclude that also for molecular systems do NVU dynamics become equivalent to NVE dynamics in the thermodynamic limit (since NVE and NVT dynamics are known to give equivalent results 29 ).