On the quantum mechanical potential of mean force. II. Constrained path integral molecular dynamics integrators
OOn the quantum mechanical potential of mean force. II. Constrained pathintegral molecular dynamics integrators
Dmitri Iouchtchenko, Kevin P. Bishop, and Pierre-Nicholas Roy a) Department of Chemistry, University of Waterloo, Waterloo, Ontario, N2L 3G1,Canada
Building on Paper I of this series, which introduced path integral Monte Carlo (PIMC) estimators for thederivative of the potential of mean force (PMF), we propose two path integral molecular dynamics (PIMD)integrators that can make use of these estimators. These integrators, c-OBABO and c-BAOAB, are basedon the path integral Langevin equation (PILE) integrator, which has seen widespread success in PIMDapplications, but they include support for holonomic constraints. When the reaction coordinate is the distancebetween two centers of mass, we find that several exact expressions are accessible: the Fixman correction,the position constraint Lagrange multiplier, and various derivatives with respect to the reaction coordinate.It is observed that c-BAOAB tends to have a smaller time step error than c-OBABO, which is consistentwith previous studies on integrator step ordering in molecular dynamics with holonomic constraints and inPIMD. Further, we show that both the PMF of a water dimer and its derivative may be obtained from PIMDsimulations using c-BAOAB, yielding results in agreement with the path integral umbrella sampling methodpreviously used for this system.
I. INTRODUCTION
In Paper I of this series, we demonstrated that thequantum potential of mean force (PMF) may be obtainedfrom constrained path integral Monte Carlo (PIMC) sim-ulations using novel estimators for the derivative of thePMF. However, devising efficient configuration updateschemes for Monte Carlo methods is a nontrivial endeavor,as the volume of configuration space grows exponentiallyboth with the number of particles and the number ofbeads per particle. For most systems, it is not feasibleto simply perturb the position of each bead by a randomamount, since the magnitude of the perturbations mustbe decreased substantially when increasing the numberof beads, in order to avoid rejecting all proposed con-figurations. Thus, it is often beneficial to turn to pathintegral molecular dynamics (PIMD), in which configu-rations are updated based on the force arising from aneffective Hamiltonian. In the unconstrained case, the pathintegral Langevin equation (PILE) integrator providesa straightforward implementation of thermostatted ringpolymer time evolution for PIMD. In the present work, we incorporate bead-local holo-nomic constraints into the PILE in order to calculate thequantum PMF of molecular systems. We also show thatwhen the reaction coordinate for the PMF is the radialseparation between two centers of mass, the additionalcomputational effort necessary to enforce the constraintis negligible, as the Lagrange multiplier may be computeddirectly without resorting to an iterative scheme.A recent study has found that an approach which com-bines PIMD with umbrella sampling and histogram un-biasing is sufficient to calculate the quantum PMF of awater dimer. While it may be advantageous to stitch a) Electronic mail: [email protected] together multiple histograms, as each simulation will con-tribute data over a range of the reaction coordinate, it canbe quite challenging to obtain and converge a collectionof smooth histograms. On the other hand, independentvalues obtained from an estimator at single points aresimple to examine and improve in a systematic fashion.Indeed, we observe that when the PMF has a very rapidlychanging slope, as in the water dimer at low tempera-ture, it can be favorable to compute its derivative usingconstrained PIMD and our estimators rather than bynumerical differentiation of an unbiased PMF.The remainder of this article is organized as follows: inSec. II, we describe our notation; in Sec. III, we developtwo integrators for constrained path integral Langevindynamics; in Sec. IV, we derive exact expressions for aspecial case of the reaction coordinate; in Sec. V, weapply the integrators to a water dimer; and in Sec. VI,we summarize our findings.
II. BACKGROUND
As in Paper I, we consider Hamiltonians of the formˆ H = 12 ˆ˜ p · ˜ M − · ˆ˜ p + V (ˆ˜ q ) (1)(note the addition of tildes to distinguish these quanti-ties from the fictitious ones in the molecular dynamicssimulations), and work with the classical potential V cl ( q ) = f (cid:88) i =1 ˜ m i P (cid:126) ˜ β P (cid:88) j =1 (cid:16) q ( j ) i − q ( j +1) i (cid:17) + 1 P P (cid:88) j =1 V ( q ( j ) ) , (2)which is extracted from the discrete imaginary time pathintegral of the quantum partition function Z = Tr e − ˜ β ˆ H at reciprocal temperature ˜ β with P beads. a r X i v : . [ phy s i c s . c h e m - ph ] J a n To sample from the unconstrained path density e − ˜ βV cl ( q ) using molecular dynamics, one may associate afictitious mass m ( j ) i and momentum p ( j ) i with each Carte-sian bead coordinate q ( j ) i . The masses form the
P f × P f diagonal mass matrix M , while the momenta are collectedinto the vector p of length P f . We require that all thefictitious masses for a single degree of freedom be equal( m ( j ) i = m i ), so that each corresponding block of M isguaranteed to be invariant under every similarity transfor-mation. We consider only the typical case when the ratio˜ m i /m i is the same for all degrees of freedom; we refer tosaid ratio as ˜ m/m and define the single-bead fictitiousmass matrix as M = m ˜ m ˜ M . (3)For a simulation at a fictitious reciprocal temperature β , the momentum partition function (cid:90) d p e − β p · M − · p = (cid:18) πβ (cid:19) Pf | M | (4)may be inserted into configurational averages to expressthem as phase space averages: (cid:82) d q e − ˜ βV cl ( q ) f ( q ) (cid:82) d q e − ˜ βV cl ( q ) = (cid:82) d p (cid:82) d q e − βH cl ( p , q ) f ( q ) (cid:82) d p (cid:82) d q e − βH cl ( p , q ) , (5)with the classical Hamiltonian given by H cl ( p , q ) = 12 p · M − · p + ˜ ββ V cl ( q ) . (6)Although it appears that we’ve increased the complex-ity of the integrals, the phase space formulation of theexpectation value is readily evaluated using moleculardynamics techniques. A. Path integral Langevin equation (PILE) integrator
The PILE integrator is a combination of a white noiseLangevin thermostat with a generalized velocity Verletscheme that operates on path normal modes. The normalmode transformation is done via an orthogonal matrix C with elements C jk = (cid:114) P if k = 0 (cid:114) P cos 2 πjkP if 1 ≤ k < P (cid:114) P ( − j if k = P (cid:114) P sin 2 πjkP if P < k ≤ P −
1, (7) where j labels the beads (1 to P ) and k labels the modes (0to P − P and Q for the transformedcoordinates.This integrator is composed of three exact sub-integrators. The propagation of harmonic oscillators innormal mode coordinates is A P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ; Q ( k ) ← (cid:80) Pj =1 q ( j ) C jk ¯ P ( k ) ← S PP k (∆ t ) P ( k ) + S PQ k (∆ t ) M Q ( k ) ¯ Q ( k ) ← S QP k (∆ t ) M − P ( k ) + S QQ k (∆ t ) Q ( k ) ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) ; ¯ q ( j ) ← (cid:80) P − k =0 C jk ¯ Q ( k ) , (8) where we have made explicit the transformations to andfrom normal modes. The propagation coefficients S PP k (∆ t ) = S QQ k (∆ t ) = cos( ω k ∆ t ) , (9a) S PQ k (∆ t ) = − ω k sin( ω k ∆ t ) , (9b)and S QP k (∆ t ) = ∆ t sinc ( ω k ∆ t ) (9c)arise from the exact solution of Hamilton’s equations ofmotion for a harmonic oscillator with angular frequency ω k = 2 (cid:115) ˜ mP (cid:126) m ˜ ββ sin πkP . (10)Application of the remaining force may be done inCartesian coordinates: B (cid:26) ¯ p ( j ) ← p ( j ) + ∆ t ˜ βP β F ( q ( j ) ) , (11) where F ( q ( j ) ) = − ∇ V ( q ( j ) ) . (12)The normal mode degrees of freedom are independentlythermostatted using a Langevin thermostat: O P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ¯ P ( k ) ← T k (∆ t ) P ( k ) + U k (∆ t ) M η ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) , (13) where the coefficients are T k (∆ t ) = e − ∆ tγ k (14a)and U k (∆ t ) = (cid:114) β (1 − e − tγ k ) , (14b) γ k is a friction coefficient (the same for all degrees offreedom), and η is a vector of f pseudorandom numberssampled from a standard normal distribution. When k ≥
1, the value of γ k that minimizes the energy auto-correlation time in the free ( V (˜ q ) = 0) case is knownanalytically to be 2 ω k , and this value is typically usedeven when interactions are present; the centroid friction γ is a tunable simulation parameter.If these sub-integrators are referred to as A, B, andO, respectively, then they may be combined in the orderOBABO to form the PILE integrator. It is implied by thisnotation that all the steps other than the central (thatis, both repetitions of B and O) have a halved durationof ∆ t/
2, as required by the symmetric splitting of theFokker–Planck operator.
B. Constrained Hamiltonian integrators
The addition of holonomic constraints to a symplecticintegrator for Hamilton’s equations of motion may beaccomplished by a straightforward scheme, in which non-dynamical momentum perturbations ensure that both theposition and velocity constraints are satisfied at the endof each step. For example, this may be used to obtainthe well-known RATTLE algorithm from velocity Verlet.In this scheme, a generic integration step of the form (cid:40) ¯ p ← f p ( p , q )¯ q ← f q ( p , q ) (15) becomes the two-step sequence p (cid:48) ← p + ∇ ξ ( q ) · Λ ¯ p ← f p ( p (cid:48) , q )¯ q ← f q ( p (cid:48) , q ) ξ (¯ q ) = z (16a) (cid:40) ¯ p ← p + ∇ ξ ( q ) · Λ ˙ ξ ( q ) = 0 , (16b) where ξ ( q ) = z is the holonomic constraint to be main-tained, ˙ ξ ( q ) = 0 is the implicit velocity constraint, andeach Λ is a vector of Lagrange multipliers that results inthe final line of the corresponding step being valid. Notehow the first step begins by perturbing the momentumaway from the constraint manifold in order to ensurethat the position constraint is satisfied, while the secondstep projects the momentum back onto the constraintmanifold.In principle, the Lagrange multipliers may be found ateach step of the simulation by integrating their equationsof motion explicitly. However, this will result in a growingdiscrepancy between the calculated values and the valuesthat are needed to correctly enforce the constraints. Instead, one should solve a system of equations for theLagrange multipliers at each step; since these equationsare generally nonlinear, they are most often solved byiteration.
III. INTEGRATORS
Our aim is to compute the derivative of the PMF, − ˜ βA (cid:48) ( ξ (cid:63) ) = (cid:82) d q δ (cid:0) ξ ( q (1) ) − ξ (cid:63) (cid:1) e − ˜ βV cl ( q ) E i ( q ) (cid:82) d q δ (cid:0) ξ ( q (1) ) − ξ (cid:63) (cid:1) e − ˜ βV cl ( q ) , (17)in a PIMD setting using the two path integral estimators E ( q ) = ∂∂ξ log | J ( q (1) ) | + ˜ βP P (cid:88) j =1 F ( q ( j ) ) · ∂ q (1) ∂ξ (18a)and E ( q ) = ∂∂ξ log | J ( q (1) ) | + ˜ β F (1)cl ( q ) · ∂ q (1) ∂ξ (18b)from Paper I. Although the molecular dynamics sim-ulations will take place in Cartesian coordinates, it issubstantially more convenient to develop the theoryusing the generalized coordinates ( X , ξ , u ) which in-clude the reaction coordinate. The transformation tothese coordinates has non-zero Jacobian determinant J ( q ) = J ( X , ξ, u ); since the unconstrained beads arenot transformed ( u = q (2) , . . . , q ( P ) ), J has no depen-dence on their values and we write simply J ( X , ξ ), notingthat it has the same value as the Jacobian determinantof the transformation on just the first bead.By taking the time derivative of the explicit constraintequation ξ ( q (1) ) = ξ (cid:63) , we find the implicit velocity con-straint ˙ ξ ( q (1) ) = 0, which prevents us from using Eq. (4)directly to obtain the necessary phase space integral. In-stead, we find that the appropriate expression is (cid:90) d p U e − β p U · A − · p U = (cid:18) πβ (cid:19) Pf − | A | = (cid:18) πβ (cid:19) Pf − | Γ | Z ξ (19a)= (cid:18) πβ (cid:19) Pf − | M | | J ( X , ξ (cid:63) ) | Z ξ ( X , ξ (cid:63) ) , (19b)where p U are the momenta conjugate to all the uncon-strained bead coordinates U = X , u , Γ = J T · M · J (20)is the generalized mass matrix, J is the Jacobian matrix(whose determinant is J ), A is Γ without the row andcolumn corresponding to ξ , and Z ξ = ∇ ξ · M − · ∇ ξ = ∇ ξ · M − · ∇ ξ (21)is the element of Γ − at the row and column correspondingto ξ . The requisite determinant identity is proved inAppendix A.Therefore, we have − ˜ βA (cid:48) ( ξ (cid:63) ) = (cid:82) d p U (cid:82) d U e − βH ccl ( p U , U ; ξ (cid:63) ) Z − ξ E i ( U , ξ (cid:63) ) (cid:82) d p U (cid:82) d U e − βH ccl ( p U , U ; ξ (cid:63) ) Z − ξ , (22)where the constrained classical Hamiltonian H ccl ( p U , U ; ξ (cid:63) ) = 12 p U · A − · p U + ˜ ββ V cl ( U , ξ (cid:63) ) (23)can be obtained from Eq. (6) by setting ˙ ξ = 0. Wemay write the above expression in terms of averages overconstrained molecular dynamics simulations as − ˜ βA (cid:48) ( ξ (cid:63) ) = (cid:104) Z − ξ E i (cid:105) ξ (cid:63) (cid:104) Z − ξ (cid:105) ξ (cid:63) , (24)where the dependence of the constrained momentum par-tition function on the coordinates has given rise to theFixman correction. These constrained molecular dy-namics simulations are in practice carried out in Cartesiancoordinates with the constraint enforced via the standardmethod of Lagrange multipliers.
A. Constrained OBABO (c-OBABO) integrator
The method described in Sec. II B is applicable toconstrained Hamiltonian systems, but unfortunately, thesimulations of interest are to be run at constant tempera-ture using a Langevin thermostat. However, as Leli`evre et al. have shown, the thermostat step may be adjusted inexactly the same manner for Langevin dynamics.
Thisallows us to transform the PILE (OBABO) integratorinto the constrained version, c-OBABO. Although forour present purposes (computing the quantum PMF), weonly require a single constraint on one bead, we derivehere a more general integrator that supports multipleindependent holonomic constraints on all beads.Consider P functions ξ j ( q ( j ) ) and constant vectors z j (not necessarily of the same length for each bead),making up the constraint equations ξ j ( q ( j ) ) = z j . Theimplicit velocity constraints, when expressed in momen-tum form, are ∇ j ξ j ( q ( j ) ) T · M − · p ( j ) = . Upon addingthese constraints to the OBABO integrator and remov-ing redundant constraint steps, we obtain the c-OBABOintegrator for PIMD: O P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ¯ P ( k ) ← T k (∆ t/ P ( k ) + U k (∆ t/ M η ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) (25a) B (cid:26) ¯ p ( j ) ← p ( j ) + ∆ t ˜ β P β F ( q ( j ) ) (25b) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = (25c) ˜ A P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ; Q ( k ) ← (cid:80) Pj =1 q ( j ) C jk ¯ P ( k ) ← S PP k (∆ t ) P ( k ) + S PQ k (∆ t ) M Q ( k ) ¯ Q ( k ) ← S QP k (∆ t ) M − P ( k ) + S QQ k (∆ t ) Q ( k ) ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) + (cid:80) P(cid:96) =1 ˜ S PP j(cid:96) (∆ t ) ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ¯ q ( j ) ← (cid:80) P − k =0 C jk ¯ Q ( k ) + (cid:80) P(cid:96) =1 ˜ S QP j(cid:96) (∆ t ) M − ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ξ j (¯ q ( j ) ) = z j (25d) B (cid:26) ¯ p ( j ) ← p ( j ) + ∆ t ˜ β P β F ( q ( j ) ) (25e) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = (25f) O P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ¯ P ( k ) ← T k (∆ t/ P ( k ) + U k (∆ t/ M η ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) (25g) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = , (25h) where ˜ S XX j(cid:96) (∆ t ) = P − (cid:88) k =0 C jk S XX k (∆ t ) C (cid:96)k (26)are the normal mode propagation coefficients transformedto real space. For the velocity constraints, the Lagrangemultipliers may be computed directly as Λ j = − (cid:16) ∇ j ξ j ( q ( j ) ) T · M − · ∇ j ξ j ( q ( j ) ) (cid:17) − · ∇ j ξ j ( q ( j ) ) T · M − · p ( j ) , (27)without an iterative scheme.Because both the normal mode transformations andthe harmonic oscillator equations of motion are linear, ingoing from A in Eq. (8) to ˜A in Eq. (25d), the constraintforce was threaded through the sub-integrator, and isapplied only at the very end of ˜A, after the inverse normalmode transformation. Although all the beads are coupledby the constraint and the Lagrange multipliers cannotbe obtained directly in the general case, in Sec. IV weshow that this form can enable direct evaluation of theLagrange multiplier when only one constraint is needed.Additionally, because the constraints couple the normalmodes even in the absence of interactions, the standardderivation for the optimal friction of a thermostattedharmonic oscillator is not applicable. Despite this, wefind that using the unmodified friction values from thePILE ( γ k = 2 ω k for k ≥
1) is a valid strategy in practice.
B. Constrained BAOAB (c-BAOAB) integrator
It has been demonstrated that the alternate integratorstep order BAOAB may result in a smaller time steperror for PIMD. While Leli`evre et al. only consider the“side” scheme (which places the thermostat on the outersides of the integrator), the feasibility and benefits ofthe “middle” scheme (which has the thermostat centeredin the integrator) have also been recently establishedfor molecular dynamics with holonomic constraints. The c-BAOAB integrator for PIMD may therefore bewritten as the following arrangement of sub-integratorsand constraint steps: B (cid:26) ¯ p ( j ) ← p ( j ) + ∆ t ˜ β P β F ( q ( j ) ) (28a) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = (28b) ˜ A P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ; Q ( k ) ← (cid:80) Pj =1 q ( j ) C jk ¯ P ( k ) ← S PP k (∆ t/ P ( k ) + S PQ k (∆ t/ M Q ( k ) ¯ Q ( k ) ← S QP k (∆ t/ M − P ( k ) + S QQ k (∆ t/ Q ( k ) ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) + (cid:80) P(cid:96) =1 ˜ S PP j(cid:96) (∆ t/ ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ¯ q ( j ) ← (cid:80) P − k =0 C jk ¯ Q ( k ) + (cid:80) P(cid:96) =1 ˜ S QP j(cid:96) (∆ t/ M − ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ξ j (¯ q ( j ) ) = z j (28c) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = (28d) O P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ¯ P ( k ) ← T k (∆ t ) P ( k ) + U k (∆ t ) M η ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) (28e) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = (28f) ˜ A P ( k ) ← (cid:80) Pj =1 p ( j ) C jk ; Q ( k ) ← (cid:80) Pj =1 q ( j ) C jk ¯ P ( k ) ← S PP k (∆ t/ P ( k ) + S PQ k (∆ t/ M Q ( k ) ¯ Q ( k ) ← S QP k (∆ t/ M − P ( k ) + S QQ k (∆ t/ Q ( k ) ¯ p ( j ) ← (cid:80) P − k =0 C jk ¯ P ( k ) + (cid:80) P(cid:96) =1 ˜ S PP j(cid:96) (∆ t/ ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ¯ q ( j ) ← (cid:80) P − k =0 C jk ¯ Q ( k ) + (cid:80) P(cid:96) =1 ˜ S QP j(cid:96) (∆ t/ M − ∇ (cid:96) ξ (cid:96) ( q ( (cid:96) ) ) · Λ (cid:96) ξ j (¯ q ( j ) ) = z j (28g) B (cid:26) ¯ p ( j ) ← p ( j ) + ∆ t ˜ β P β F ( q ( j ) ) (28h) C ¯ p ( j ) ← p ( j ) + ∇ j ξ j ( q ( j ) ) · Λ j ∇ j ξ j ( q ( j ) ) T · M − · ¯ p ( j ) = . (28i) Our rather unsophisticated implementation of c-BAOAB suffers an increase in run time on the orderof 10 % compared to c-OBABO due to the additionalwork required. Despite this, as we show in Sec. V, theformer may still outperform the latter, since it can allowa much larger time step to be used. IV. EXACT RELATIONS FOR A CENTER OF MASSDISTANCE CONSTRAINT
When the reaction coordinate ξ is the radial separationbetween two centers of mass, several useful expressionsmay be derived. From this point onward, we work ex-plicitly with N particles in 3-dimensional space, withcoordinates x i (where x describes all the degrees of free-dom at the first bead, and can be thought of as a morestructured reinterpretation of q (1) ). The convex sums x α = N α (cid:88) i =1 m αi m α x αi (29a)and x β = N β (cid:88) i =1 m βi m β x βi (29b)are the centers of mass of α and β (which we requireto be non-empty and disjoint), where m αi is the masscorresponding to x αi and m α = N α (cid:88) i =1 m αi , (30)and likewise for β . We also use x γi to refer to the remain-ing particles, if any, which do not participate in eithercenter of mass, with N α + N β + N γ = N .We constrain the distance between the centers of massas ξ ( x ) = | r | = | x α − x β | = z, (31)for an arbitrary positive z . The derivatives of this functionare ∇ αi ξ ( x ) = m αi r m α ξ ( x ) , (32a) ∇ βi ξ ( x ) = − m βi r m β ξ ( x ) , (32b)and ∇ γi ξ ( x ) = 0 . (32c)Immediately, it follows that Eq. (21) evaluates to Z ξ ( x ) = N α (cid:88) i =1 m αi | ∇ αi ξ ( x ) | + N β (cid:88) i =1 m βi | ∇ βi ξ ( x ) | (33a)= 1 m α + 1 m β = 1 µ αβ , (33b)which is a constant (with µ αβ being the reduced mass),so the Fixman correction cancels from Eq. (24) leavingus with just − ˜ βA (cid:48) ( ξ (cid:63) ) = (cid:104)E i (cid:105) ξ (cid:63) . (34) A. Lagrange multiplier
It is also straightforward in this circumstance to obtaina closed-form solution for the Lagrange multiplier λ forthe constrained harmonic oscillator propagation step ˜A.Because the momentum perturbation was propagatedthrough the linear equations of motion, we know that theequation to be satisfied is ξ (cid:16) ¯ x + λ ˜ S QP M − ∇ ξ ( x ) (cid:17) = z, (35)where x is the value of q (1) before the propagation, and¯ x is its value after propagation by the unconstrained stepA. After some algebraic manipulations, we get (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:32) λ ˜ S QP µ αβ ξ ( x ) (cid:33) r + ∆ r (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = z, (36)where ∆ r = ¯ r − r is the change in relative position ofthe centers of mass due to the unconstrained step. Weassume that ˜ S QP is not zero.From geometric considerations, we see that solutionsexist only when the component of ∆ r orthogonal to r does not surpass z : | ∆ r | − (∆ r · ˆr ) ≤ z , (37)where ˆr = r / | r | is the unit vector in the direction of r .When this condition is met, the Lagrange multiplier isgiven by λ = − µ αβ ˜ S QP (cid:34) | r | + ∆ r · ˆr − (cid:114) z − (cid:16) | ∆ r | − (∆ r · ˆr ) (cid:17)(cid:35) , (38)which may be computed directly for very little cost. Whenthe square root is not zero, there is another solution for λ ,which has the root being added rather than subtracted,but this corresponds to the interchange of the two centersof mass; we assume that the time step is sufficiently smallthat this will never be the desired outcome. The interpretation of the expression for λ is simple. Theresulting shifts of the particles in α and β are, respectively, − m β ˆr m α + m β (cid:34) | r | + ∆ r · ˆr − (cid:114) z − (cid:16) | ∆ r | − (∆ r · ˆr ) (cid:17)(cid:35) (39a)and m α ˆr m α + m β (cid:34) | r | + ∆ r · ˆr − (cid:114) z − (cid:16) | ∆ r | − (∆ r · ˆr ) (cid:17)(cid:35) . (39b)The first two terms of each are responsible for completelyremoving any separation between the centers of massalong r : the first term takes care of the original vector,while the second term handles the component of ∆ r thatis parallel to r . The prefactors ensure that the groupsmove toward one another, and that the lighter groupmoves farther. The square root term then restores someof the separation along r , with the exact amount beingdetermined by the excess allowance from Eq. (37).It is important to note that Eq. (38) is not an equationof motion for the constraint. Rather, this is an exact so-lution to the Lagrange multiplier optimization, equivalentto one that would be obtained by an iterative scheme.As such, it is not susceptible to numerical drift of theconstrained coordinate. B. Derivatives for estimators
Evaluation of the estimators E and E requires knowl-edge of ∂∂ξ log | J ( x ) | and ∂ x ∂ξ , which we find by performingseveral coordinate transformations. We begin by trans-forming x into the Jacobi coordinates y α(cid:96) = (cid:96) (cid:88) i =1 m αi m → (cid:96)α x αi − (cid:26) x α,(cid:96) +1 if 1 ≤ (cid:96) ≤ N α −
10 if (cid:96) = N α , (40)with m → (cid:96)α = (cid:96) (cid:88) i =1 m αi , (41)and similarly for β . The remaining coordinates x γi are leftunmodified. As shown in Appendix B, this transformationhas unit Jacobian determinant.The coordinates y αN α and y βN β are the centers of mass x α and x β , so we further transform them to R = m α y αN α + m β y βN β m α + m β (42a)and r = y αN α − y βN β , (42b)and this change of variables also has unit Jacobian deter-minant.Finally, the transformation of r to the spherical coor-dinates ( ξ , cos θ , ϕ ) is known to have a Jacobian deter-minant whose absolute value is ξ . The overall transfor-mation from x to ( y αi , y βi , R , ξ , cos θ , ϕ , x γi ) thereforehas | J ( x ) | = ξ , so ∂∂ξ log | J ( x ) | = 2 ξ . (43)The original Cartesian coordinates may be written as x αi = R + m β r m α + m β + N α − (cid:88) (cid:96) = i m α,(cid:96) +1 m → (cid:96) +1 α y α(cid:96) − i = 1 m → i − α m → iα y α,i − if 2 ≤ i ≤ N α (44a)and x βi = R − m α r m α + m β + N β − (cid:88) (cid:96) = i m β,(cid:96) +1 m → (cid:96) +1 β y β(cid:96) − i = 1 m → i − β m → iβ y β,i − if 2 ≤ i ≤ N β . (44b)Hence, the derivatives of x with respect to ξ are ∂ x αi ∂ξ = m β m α + m β ˆr , (45a) ∂ x βi ∂ξ = − m α m α + m β ˆr , (45b)and ∂ x γi ∂ξ = 0 . (45c)As observed in the shifts to enforce the constraints inEq. (39), the prefactor for the lighter group is largerin magnitude than for heavier group: the former mustmove faster under a changing separation distance thanthe latter. V. RESULTS
To demonstrate the effectiveness of the constrainedPIMD integrators, we apply them to a q-SPC/Fw wa-ter dimer and obtain its quantum PMF as a functionof the distance ξ ( x ) = | x α − x β | between the molecules’centers of mass x α and x β . To generate reference re-sults, we use the path integral umbrella sampling method(US/WHAM), which requires a histogram unbiasing stepin order to stitch together the obtained histograms. − − − − ˜ β A / n m − T = 10 K , ξ ? = 0 . c-OBABO, estimator 1c-OBABO, estimator 2 c-BAOAB, estimator 1c-BAOAB, estimator 2 − − ∆ t/ ps − ˜ β A / n m − T = 300 K , ξ ? = 0 .
25 nm
FIG. 1. Comparison of the time step convergence of the PMFderivative of a water dimer using the integrators c-OBABOand c-BAOAB in Eqs. (25) and (28) at T = 10 K, ξ (cid:63) = 0 . T = 300 K, ξ (cid:63) = 0 .
25 nm (bottom panel). Linesegments added to guide the eye.
Before attempting to generate a PMF, we first identifywhich integrator has better time step error characteristics.In Fig. 1, we show two combinations of temperature andconstraint distance that have an appreciable differencein behavior between c-OBABO (solid) and c-BAOAB(dashed). In both cases, the time step error decreasesfaster for the latter, which is to be expected from pastwork with this ordering of integrator steps.
For theremainder of the calculations, we only use c-BAOAB, asit allows us to save some computational effort by using alarger time step to achieve the same level of error.As shown in Fig. 2, direct estimation of the derivative ofthe PMF is successful using the c-BAOAB integrator andboth estimators. The reference curves are obtained bynumerically differentiating the US/WHAM PMF, whichcauses a slight wobble that is noticeable at 300 K. At 10 K,the PIMD simulations show that the derivative changesrather sharply near 0 . . . . . ξ ? / nm − − − − ˜ β A / n m − T = 300 K T = 10 K Estimator 1Estimator 2
FIG. 2. Derivative of the PMF of a water dimer at T = 300 K(top curve, less saturated) and 10 K (bottom curve, moresaturated). Error bars are not visible, because they are smallerthan the symbols. The solid curves show the US/WHAMresults. . . . . ξ ? / nm σ / n m − T = 10 K T = 300 K Estimator 1Estimator 2
FIG. 3. Comparison of the standard error of the mean ofthe PMF derivatives in Fig. 2 at T = 300 K (lower points, lesssaturated) and 10 K (upper points, more saturated). without any interesting features. However, the picture isnot quite as simple at 10 K, which has a crossing preciselyin the problematic region near ξ (cid:63) = 0 . − ˜ βA (cid:48) values and plot the renormalized PMF ˜ A in Fig. 4, usingthe same procedure as in Paper I. There is excellentagreement with the reference results at both temperatures,suggesting that this approach of integrating the derivativeof the PMF is a viable alternative to umbrella samplingand histogram unbiasing. VI. CONCLUSIONS
We have augmented the PILE integrator with bead-local holonomic constraints to produce the c-OBABO . . . . ξ ? / nm − − − ˜ A / k J m o l − T = 300 K T = 10 K Estimator 1Estimator 2
FIG. 4. Potential of mean force of a water dimer at T = 300 K(top curve, less saturated) and 10 K (bottom curve, moresaturated). Error bars are not visible, because they are smallerthan the symbols. Additional points extending to ξ (cid:63) = 1 . and c-BAOAB integrators for PIMD. These integratorsallow our estimators from Paper I to be used with molec-ular dynamics in addition to Monte Carlo, which greatlyexpands their scope of applicability. The constrained har-monic oscillator propagation step ˜A of these integratorshas the constraint force applied at the end rather thanthe beginning of the step, which allows us to find an exactexpression for the Lagrange multiplier for some reactioncoordinates, such as the distance between two centers ofmass.Using one of the constrained integrators, we have com-puted the derivative of the PMF of a water dimer at10 K and 300 K. We observe that the constrained PIMDmethod captures the sharp step in the derivative bet-ter than umbrella sampling. We have also been able tosuccessfully integrate the derivative to obtain the PMFitself.Together, these novel estimators and integrators maybe utilized to find free energy differences and potentialsof mean force of molecular clusters using implementa-tions based on existing PIMD simulation software. Theconventional umbrella sampling and histogram unbiasingapproach has several practical drawbacks: the user mustchoose a force constant for the restraint; the ends of thehistogram have a tendency to be noisy, so histogramsmust be generated well past the region of interest; in ad-dition to requiring a grid of points for the unbiasing step,the user must also select the locations of the restraintwindows and make sure that the resulting histogramshave sufficient overlap; the iterative unbiasing procedurerequires a stopping criterion. Our method, on the otherhand, allows the result of each PIMD simulation to beconverged independently of all the others, and there areno additional tunable parameters introduced into thesimulations besides the choice of integration grid. ACKNOWLEDGMENTS
This research was supported by the Natural Sciencesand Engineering Research Council of Canada (NSERC)(RGPIN-2016-04403), the Ontario Ministry of Researchand Innovation (MRI), the Canada Research Chair pro-gram (950-231024), and the Canada Foundation for Inno-vation (CFI) (project No. 35232).
Appendix A: Block matrix determinant identity
Consider an invertible block matrix and its inverse withidentical block structure, (cid:18)
A BC D (cid:19) = (cid:18) W XY Z (cid:19) − , (A1)whose blocks A and Z are themselves invertible. Uponcombining two manifestations of the Schur complement, (cid:12)(cid:12)(cid:12)(cid:12) A BC D (cid:12)(cid:12)(cid:12)(cid:12) = | A | (cid:12)(cid:12) D − CA − B (cid:12)(cid:12) (A2)and Z − = D − CA − B , (A3)we conclude that (cid:12)(cid:12)(cid:12)(cid:12) A BC D (cid:12)(cid:12)(cid:12)(cid:12) = | A | | Z | − . (A4) Appendix B: Jacobian determinant of transformation toJacobi coordinates
The transformation to Jacobi coordinates described inSec. IV B couples neither α and β nor the spatial degreesof freedom, and also does not involve γ , so its Jacobianmatrix will be block diagonal with 7 independent blocks.The entire γ block is irrelevant, as its determinant is1. The remaining 6 blocks are identical in form, so wemay scrutinize only one of them, using x i to label anarbitrary Cartesian component of either x αi or x βi , withcorresponding masses m i and transformed coordinates y (cid:96) .Trivially, the derivative of y (cid:96) with respect to x i is ∂y (cid:96) ∂x i = m i m → (cid:96) if i ≤ (cid:96) − i = (cid:96) + 10 otherwise, (B1) so the Jacobian matrix J with elements J (cid:96)i = ∂y (cid:96) ∂x i (B2)is a lower Hessenberg matrix (all elements above thesuperdiagonal i = (cid:96) + 1 are zero). Thus, we may use arecurrence relation to find | J | . In the following, J ( n ) denotes the n × n leading principalsubmatrix of J , and we see that | J (0) | = 1 and | J (1) | = 1.For subsequent n , | J ( n ) | = J nn | J ( n − | + n − (cid:88) i =1 ( − n − i J ni n − (cid:89) j = i J j,j +1 | J ( j − | (B3a)= m n m → n | J ( n − | + n − (cid:88) i =1 m i m → n n − (cid:89) j = i | J ( j − | . (B3b)By induction, it stands to reason that | J ( n ) | = 1 for all n , and therefore | J | = 1, so the entire transformationto Jacobi coordinates for both centers of mass has unitJacobian determinant. This is consistent with the view ofJacobi coordinates as iterated pairwise transformationsto center of mass and relative distance coordinates, whichindividually have unit Jacobian determinant. M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopou-los, J. Chem. Phys. , 124104 (2010). K. P. Bishop and P.-N. Roy, J. Chem. Phys. , 102303 (2018). M. Parrinello and A. Rahman, J. Chem. Phys. , 860 (1984). S. Reich, SIAM J. Numer. Anal. , 475 (1996). H. C. Andersen, J. Comput. Phys. , 24 (1983). J.-P. Ryckaert, G. Ciccotti, and H. J. Berendsen, J. Comp. Phys. , 327 (1977). M. Fixman, J. Chem. Phys. , 1527 (1978). W. K. den Otter, J. Chem. Theory Comput. , 3861 (2013). T. Leli`evre, M. Rousset, and G. Stoltz,
Free Energy Compu-tations: A Mathematical Perspective (Imperial College Press,2010). T. Leli`evre, M. Rousset, and G. Stoltz, Math. Comput. , 2071(2012). J. Liu, D. Li, and X. Liu, J. Chem. Phys. , 024103 (2016). Z. Zhang, X. Liu, K. Yan, M. E. Tuckerman, and J. Liu, J. Phys.Chem. A , 6056 (2019). D. Iouchtchenko, “WaterMeanMeanForce,” https://github.com/0/WaterMeanMeanForce.jl , accessed 2021-01-03. F. Paesani, W. Zhang, D. A. Case, T. E. Cheatham III, andG. A. Voth, J. Chem. Phys. , 184507 (2006). R. W. Cottle, Lin. Algebra Appl. , 189 (1974). N. D. Cahill, J. R. D’Errico, D. A. Narayan, and J. Y. Narayan,Coll. Math. J.33