On the quantum mechanical potential of mean force. I. A path integral perspective
OOn the quantum mechanical potential of mean force. I. A path integralperspective
Dmitri Iouchtchenko, Kevin P. Bishop, and Pierre-Nicholas Roy a) Department of Chemistry, University of Waterloo, Waterloo, Ontario, N2L 3G1,Canada
We derive two path integral Monte Carlo estimators for the derivative of the potential of mean force (PMF),which may be numerically integrated to yield the PMF. For the first estimator, we perform the differentiationon the exact path integral, and for the second, we perform the differentiation on the path integral afterdiscretization. These estimators are successfully validated against reference results for the harmonic oscillatorand Lennard-Jones dimer systems using constrained path integral Monte Carlo (PIMC) simulations. Specifically,the estimators reproduce both the derivative of the PMF, as well as the PMF itself, for the model systemsat multiple temperatures. In Paper II, these estimators are implemented alongside path integral moleculardynamics (PIMD) with a constrained path integral Langevin equation thermostat for use with more generalsystems and potentials.
I. INTRODUCTION
Free energy calculations provide vital theoretical in-sights into the behaviour of chemical systems (such asequilibrium constants and stability of molecular confor-mations) and are commonly used to make comparisonswith experimental results. The classical potential of meanforce (PMF), also referred to as a free energy profile,may be obtained from classical molecular dynamics sim-ulations in a number of ways: umbrella sampling withthe weighted histogram analysis method (WHAM), blue moon sampling, metadynamics, and potentialof mean constraint force, among others.The quantum mechanical analogue of the PMF hasnot been as thoroughly studied, but there are many in-teresting systems where nuclear quantum effects play anintegral role. These effects are essential within simulationsperformed at low temperature or containing light atoms.Crucially, the inclusion of nuclear quantum effects hasbeen demonstrated to be indispensable in the accuratedetermination of properties of the water dimer as wellas other small water clusters due to the presence of thelight hydrogen atoms.
Recent work has combined theexisting umbrella sampling method with path integralmolecular dynamics simulations to study the free energyprofile of water–water and water–methanol dimers while accounting for such nuclear quantum effects.Many existing efforts to compute the quantumPMF have used the path integral centroidcoordinate, but this is known to produce deviationsfrom the exact quantum mechanical result. Our cur-rent focus is therefore to obtain the PMF as a function ofthe true reaction coordinate rather than the path centroid.In the present work, we formally derive two path inte-gral Monte Carlo (PIMC) estimators for the derivativeof the PMF, which can be integrated to determine thePMF. The first estimator is obtained by performing an a) Electronic mail: [email protected] analytical differentiation of the exact path integral, fol-lowed by its discretization over the path. Conversely, thesecond estimator is derived by initially discretizing thepath integral before performing the analytical differen-tiation. Theoretically, these estimators should providethe same numerical results, and to verify this, they arebenchmarked against known model systems.Both estimators may be used in conjunction with thepath integral Langevin equation (PILE), as we showin Paper II of this series, titled “Constrained path in-tegral molecular dynamics integrators”. Together withthe constrained PILE formulation, we aim to use theseestimators for systems that are not as easily studied withPIMC, such as low-temperature molecular clusters.The remainder of this article is organized as follows: inSec. II, we describe our notation; in Sec. III, we developtwo estimators for the derivative of the PMF; in Sec. IV,we apply the estimators to model systems; and in Sec. V,we summarize our findings. II. BACKGROUND
We consider systems with f Cartesian degrees of free-dom, that we label q , q , . . . , q f ; commonly, there are N particles in three spatial dimensions, in which case f = 3 N . For convenience, we group them into a singlevector q .We restrict the Hamiltonian to have the formˆ H = ˆ K + ˆ V = f (cid:88) i =1 ˆ p i m i + V (ˆ q ) = 12 ˆ p · M − · ˆ p + V (ˆ q ) , (1)where ˆ p i is the momentum operator conjugate to theposition operator ˆ q i , m i is the mass corresponding to q i ,and M is the diagonal mass matrix whose elements are m i . The restriction on the kinetic energy allows us to a r X i v : . [ phy s i c s . c h e m - ph ] J a n write the exact free particle propagator (cid:104) q (cid:48) | e − τ ˆ K | q (cid:105) = (cid:115) | M | (2 π (cid:126) τ ) f e − (cid:126) τ ( q (cid:48) − q ) · M · ( q (cid:48) − q ) (2)for an imaginary time duration τ . We require that thepotential energy be diagonal in the position representationso that (cid:104) q (cid:48) | e − τ ˆ V | q (cid:105) = δ ( q (cid:48) − q ) e − τV ( q ) . (3)Despite these limitations, such Hamiltonians are generalenough to describe many diverse systems of itinerantparticles.The partition function of a system with Hamiltonianˆ H at reciprocal temperature β = 1 /k B T is Z = Tr e − β ˆ H = (cid:90) d q (cid:104) q | e − β ˆ H | q (cid:105) , (4)and the thermal expectation value of an operator ˆ O is (cid:104) ˆ O (cid:105) β ˆ H = 1 Z Tr e − β ˆ H ˆ O = (cid:82) d q (cid:104) q | e − β ˆ H ˆ O | q (cid:105) (cid:82) d q (cid:104) q | e − β ˆ H | q (cid:105) . (5)As a means of evaluating (cid:104) ˆ O (cid:105) β ˆ H , we may construct adiscretized imaginary time path integral for the partitionfunction. To that end, we first rename q to Q (1) and theninsert P − = (cid:90) d Q ( j ) | Q ( j ) (cid:105)(cid:104) Q ( j ) | , (6)which introduce the additional Cartesian coordinates Q (2) ,. . . , Q ( P ) along the imaginary time path; we combinethem all into the vector Q and refer to them as “beads”,picturing the path as a necklace. This results in Z = (cid:90) d Q P (cid:89) j =1 (cid:104) Q ( j ) | e − βP ˆ H | Q ( j +1) (cid:105) , (7)where it should be understood that the path is cyclic inimaginary time (that is, Q ( P +1) is an alias for Q (1) ).To evaluate each high-temperature propagator, since[ ˆ K, ˆ V ] (cid:54) = 0, we rely on the Trotter factorization e − β ˆ H = lim P →∞ (cid:16) e − βP ˆ K e − βP ˆ V (cid:17) P , (8)which allows us to start with the approximation (cid:104) q (cid:48) | e − βP ˆ H | q (cid:105) ≈ (cid:115) | M | P f (2 π (cid:126) β ) f e − P (cid:126) β ( q (cid:48) − q ) · M · ( q (cid:48) − q ) − βP V ( q ) (9)and systematically improve the error in the product ofthese approximate factors by increasing P . For any finite P , we may construct the approximate path density π ( Q ) = (cid:18) | M | P f (2 π (cid:126) β ) f (cid:19) P e − βV cl ( Q ) (10) 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 ) ) , (11)so that Z = lim P →∞ (cid:90) d Q π ( Q ) . (12)For the remainder of this article, we drop the P → ∞ limit for the sake of brevity.In order to use PIMC sampling to calculate (cid:104) ˆ O (cid:105) β ˆ H , itis necessary to procure an estimator function E ˆ O ( Q ), thedetails of which depend on the nature of the operator.The operator expression in Eq. (5) is then replaced by aratio of integrals containing only regular functions: (cid:104) ˆ O (cid:105) β ˆ H = (cid:104)E ˆ O (cid:105) π = (cid:82) d Q π ( Q ) E ˆ O ( Q ) (cid:82) d Q π ( Q ) . (13)This ratio is commonly evaluated as (cid:104)E ˆ O (cid:105) π ≈ N MC N MC (cid:88) i =1 E ˆ O ( Q [ i ] ) (14)by drawing the samples { Q [ i ] } N MC i =1 from π ( Q ) usingMarkov chain Monte Carlo. III. ESTIMATORS
It is generally more convenient to work with path inte-grals in Cartesian coordinates q , but the PMF A ( ξ (cid:63) ) isexpressed in terms of an arbitrary curvilinear coordinate ξ at some value ξ (cid:63) . To connect the two, we introduce aninvertible coordinate transformation to the generalizedcoordinates X , X , . . . , X f − , ξ , where the first f − X . This transformationhas non-zero Jacobian determinant J ( q ) = J ( X , ξ ). Thespecial coordinate ξ is referred to as the reaction coordi-nate ; for example, it may be the distance between twospecific centers of mass in a cluster.Using the diagonal reduced density (cid:37) ( ξ (cid:63) ) = 1 Z (cid:104) ξ (cid:63) | Tr X e − β ˆ H | ξ (cid:63) (cid:105) (15a)= 1 Z (cid:90) d X (cid:104) X ξ (cid:63) | e − β ˆ H | X ξ (cid:63) (cid:105) (15b)at reciprocal temperature β , we may construct the overallobject of interest: the potential of mean force A ( ξ (cid:63) ) = − β log (cid:37) ( ξ (cid:63) ) (cid:37) , (16)where (cid:37) is an arbitrary constant with the same physicaldimension as ξ − . Choosing a value for (cid:37) sets the zero ofenergy for the PMF. Note that our definitions imply that (cid:37) − = (cid:82) d ξ (cid:63) e − βA ( ξ (cid:63) ) , which does not contain any explicitvolume element factors; instead, we encounter a geometricterm in the estimators. Even though the momentumoperator ˆ p ξ conjugate to the reaction coordinate operatorˆ ξ satisfies (cid:104) ξ | ˆ p ξ = − i (cid:126) ∂∂ξ (cid:104) ξ | , (17)in general we find that (cid:104) q | ˆ p ξ (cid:54) = − i (cid:126) ∂∂ξ (cid:104) q | , (18)and the missing portion is directly responsible for thegeometric term.We wish to compute A ( ξ (cid:63) ) via its derivative A (cid:48) ( ξ (cid:63) ) = − β dd ξ (cid:63) log (cid:37) ( ξ (cid:63) ) (cid:37) = − β (cid:37) (cid:48) ( ξ (cid:63) ) (cid:37) ( ξ (cid:63) ) . (19)As shown in Appendix A, we may write the diagonalreduced density in Cartesian coordinates as (cid:37) ( ξ (cid:63) ) = 1 Z (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | e − β ˆ H | q (cid:105) , (20)so − βA (cid:48) ( ξ (cid:63) ) = dd ξ (cid:63) (cid:82) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | e − β ˆ H | q (cid:105) (cid:82) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | e − β ˆ H | q (cid:105) . (21)Because the denominator resembles a constrained versionof the partition function Z in Eq. (4), we use this asthe starting point to derive two path integral estimators E ( Q ) and E ( Q ) which satisfy − βA (cid:48) ( ξ (cid:63) ) = (cid:104)E i (cid:105) π,ξ (cid:63) = (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) E i ( Q ) (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) . (22)For E ( Q ), we first differentiate and then discretize thepath integral in the numerator, and for E ( Q ) we do thereverse. A. Estimator 1: Differentiate then discretize
The main quantity in question is
Z(cid:37) (cid:48) ( ξ (cid:63) ) = dd ξ (cid:63) (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | e − β ˆ H | q (cid:105) , (23)which we write using Appendix B as Z(cid:37) (cid:48) ( ξ (cid:63) ) = (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:20) J ξ ( q ) + ∂∂ξ (cid:21) (cid:104) q | e − β ˆ H | q (cid:105) . (24) The simpler of the two terms is the geometric one, whichstems from the coordinate transformation: (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | e − β ˆ H | q (cid:105) J ξ ( q ) . (25)The remaining term (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) f (cid:88) i =1 ∂q i ∂ξ G i ( q ) (26)is more involved, requiring the derivatives of the imaginarytime propagator: G i ( q ) = ∂∂q i (cid:104) q | e − β ˆ H | q (cid:105) = 1 i (cid:126) (cid:104) q | (cid:104) e − β ˆ H , ˆ p i (cid:105) | q (cid:105) , (27)where the derivative–commutator identity is derived inAppendix C.Using the Kubo formula for the commutator with theexponential of an operator, we find that G i ( q ) = − i (cid:126) (cid:90) β d λ (cid:104) q | e − ( β − λ ) ˆ H [ ˆ H, ˆ p i ] e − λ ˆ H | q (cid:105) . (28)Since the kinetic energy operator in Eq. (1) commuteswith ˆ p i , only the commutator with the potential energyremains: [ ˆ V , ˆ p i ] . It follows from [ˆ q i , ˆ p j ] = i (cid:126) δ ij that[ V (ˆ q ) , ˆ p i ] = − i (cid:126) F i (ˆ q ) , (29)where the component of the force vector F ( q ) on thecoordinate q i is given by F i ( q ) = − ∂∂q i V ( q ) . (30)Hence, G i ( q ) = (cid:90) β d λ (cid:104) q | e − ( β − λ ) ˆ H F i (ˆ q ) e − λ ˆ H | q (cid:105) . (31)Having obtained the necessary expressions, we maydiscretize the path integral in the usual fashion. Thegeometric term in Eq. (25) poses no difficulty, and we get (cid:90) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) J ξ ( Q (1) ) . (32)The integral from the Kubo formula is discretized into anaverage over the path, and because F i (ˆ q ) is diagonal inthe additional path coordinates, we only need to performthe substitution G i ( q ) → π ( Q ) βP P (cid:88) j =1 F i ( Q ( j ) ) , (33)turning Eq. (26) into (cid:90) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) βP P (cid:88) j =1 F ( Q ( j ) ) · ∂ Q (1) ∂ξ . (34)Thus, the PMF derivative may be written as − βA (cid:48) ( ξ (cid:63) ) = (cid:104)E (cid:105) π,ξ (cid:63) = (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) E ( Q ) (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) , (35)where E ( Q ) = ∂∂ξ log | J ( Q (1) ) | + βP P (cid:88) j =1 F ( Q ( j ) ) · ∂ Q (1) ∂ξ (36)is the first estimator. In the P = 1 case, it reduces to aform recognizable from classical mechanics: − β E ( q ) = − β ∂∂ξ log | J ( q ) | − F ( q ) · ∂ q ∂ξ (37a)= ∂∂ξ (cid:20) V ( q ) − β log | J ( q ) | (cid:21) . (37b) B. Estimator 2: Discretize then differentiate
To derive another estimator, we start from Eq. (24),but first discretize the path into P imaginary time stepsto find Z(cid:37) (cid:48) ( ξ (cid:63) ) = (cid:90) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) × (cid:34) J ξ ( Q (1) ) + ∂ Q (1) ∂ξ · ∂∂ Q (1) (cid:35) π ( Q ) . (38)The geometric term will again be as in Eq. (32), but theother term is now straightforward to compute via ordinarycalculus, requiring only1 π ( Q ) ∂π ( Q ) ∂ Q ( j ) = − β ∂V cl ( Q ) ∂ Q ( j ) = β F ( j )cl ( Q ) = βP F ( Q ( j ) ) − P (cid:126) β M · (cid:104) Q ( j ) − Q ( j +1) − Q ( j − (cid:105) , (39)in which the classical force F ( j )cl ( Q ) on bead j is obtainedfrom the classical potential. The PMF derivative maytherefore also be written as − βA (cid:48) ( ξ (cid:63) ) = (cid:104)E (cid:105) π,ξ (cid:63) = (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) E ( Q ) (cid:82) d Q δ (cid:16) ξ ( Q (1) ) − ξ (cid:63) (cid:17) π ( Q ) , (40)where E ( Q ) = ∂∂ξ log | J ( Q (1) ) | + β F (1)cl ( Q ) · ∂ Q (1) ∂ξ (41)is the second estimator.It is perhaps a little surprising that the sum βP (cid:80) Pj =2 F ( Q ( j ) ) from E , which involves all coordinatesexcept the constrained one, appears to have been replacedby − P (cid:126) β M · (cid:16) Q (1) − Q (2) − Q ( P ) (cid:17) , which depends ononly three coordinates. This results in the peculiar iden-tity (cid:42) F (1)cl ( Q ) · ∂ Q (1) ∂ξ (cid:43) π,ξ (cid:63) P →∞ = (cid:42) P P (cid:88) j =1 F ( Q ( j ) ) · ∂ Q (1) ∂ξ (cid:43) π,ξ (cid:63) , (42)which relates the classical force on the constrained coor-dinates to the average force over the path. C. Removal of geometric term
It is occasionally more convenient to work with˜ (cid:37) ( ξ (cid:63) ) = (cid:37) ( ξ (cid:63) ) f ( ξ (cid:63) ) , (43)for some function f , than with (cid:37) ( ξ (cid:63) ) directly. Consider,for example, the spherical coordinates ( ξ , cos θ , ϕ ), whichhave Jacobian determinant J ( ξ, cos θ, ϕ ) = − ξ . (44)The normalization1 = (cid:90) d ξ (cid:63) ( ξ (cid:63) ) ˜ (cid:37) ( ξ (cid:63) ) (45)is often more natural than1 = (cid:90) d ξ (cid:63) (cid:37) ( ξ (cid:63) ) . (46)Using ˜ (cid:37) leads to the modified PMF˜ A ( ξ (cid:63) ) = − β log ˜ (cid:37) ( ξ (cid:63) )˜ (cid:37) = A ( ξ (cid:63) ) + 1 β log ˜ (cid:37) f ( ξ (cid:63) ) (cid:37) , (47)where the arbitrary constant ˜ (cid:37) has the same physicaldimension as ( ξf ( ξ )) − . The corresponding derivative is − β ˜ A (cid:48) ( ξ (cid:63) ) = − βA (cid:48) ( ξ (cid:63) ) − ∂∂ξ log f ( ξ (cid:63) ) , (48)and it follows immediately that˜ E ( Q ) = E ( Q ) − ∂∂ξ log f ( ξ (cid:63) ) (49a)and ˜ E ( Q ) = E ( Q ) − ∂∂ξ log f ( ξ (cid:63) ) (49b)may be used to estimate − β ˜ A (cid:48) ( ξ (cid:63) ).Whenever a transformation from q to X , ξ exists withJacobian determinant J ( X , ξ ) that is a function of only ξ (as in the above spherical coordinates example) thegeometric term may be exactly cancelled from E and E by setting f ( ξ ) = | J ( ξ ) | . In such a situation, ∂∂ξ log f ( ξ (cid:63) ) = ∂∂ξ log | J ( Q (1) ) | , (50)so we are left with just˜ E ( Q ) = βP P (cid:88) j =1 F ( Q ( j ) ) · ∂ Q (1) ∂ξ (51a)and ˜ E ( Q ) = β F (1)cl ( Q ) · ∂ Q (1) ∂ξ . (51b)This modification provides no substantial computationalbenefits, as the omitted expression will be a constant withrespect to the integration (for example, 2 /ξ (cid:63) for sphericalcoordinates). However, it may make sense to exclude theterm from the calculation entirely if it is destined to beexcised after the calculation is completed. D. Kubo formula in generalized coordinates
Starting from Eq. (15), which expresses the diagonalreduced density in terms of the generalized coordinates,application of Appendix C immediately yields
Z(cid:37) (cid:48) ( ξ (cid:63) ) = 1 i (cid:126) (cid:90) d X (cid:104) X ξ (cid:63) | (cid:104) e − β ˆ H , ˆ p ξ (cid:105) | X ξ (cid:63) (cid:105) . (52)This bypasses many of the convoluted steps found aboveand leaves us with a succinct expression, which takes onthe form − i (cid:126) (cid:90) d X (cid:90) β d λ (cid:104) X ξ (cid:63) | e − ( β − λ ) ˆ H [ ˆ H, ˆ p ξ ] e − λ ˆ H | X ξ (cid:63) (cid:105) (53)after treatment with the Kubo formula. Since[ V (ˆ q ) , ˆ p ξ ] = − i (cid:126) F ξ (ˆ q ) (54)involves only the force along the reaction coordinate, pro-ceeding in this direction seems like the obvious choice.However, ˆ p ξ is not guaranteed to commute with the Carte-sian momenta, and the commutator [ ˆ K, ˆ p ξ ] is not alwaysdiagonal in the position representation. Consequently,Eq. (53) does not lend itself well to discretization and wedo not pursue this approach to the derivation further. IV. RESULTS
As a proof of concept, we use the estimators E and E to compute derivatives of the PMF of two small systems,for which reference results (either exact or numerical) maybe calculated: a one-dimensional harmonic oscillator anda Lennard-Jones model of the Ar dimer. To perform thepath integral Monte Carlo sampling, we have implementeda basic Markov chain integrator using the Metropolis–Hastings acceptance criterion. The constraint is exactlyenforced by sampling in the generalized coordinates forthe first bead: updates are proposed for X (1) , but ξ (1) isheld fixed at ξ (cid:63) . A. Harmonic oscillator
The simplest non-trivial problem we can consider is thedependable harmonic oscillator, with the Hamiltonianˆ H = ˆ p m + 12 mω ˆ q (55)and the reaction coordinate ξ = q . Since the eigenstatesof this Hamiltonian are known analytically, we may writedown the normalized diagonal density (cid:37) ( ξ (cid:63) ) = 1 Z e − β (cid:126) ω α √ π e − ( αξ (cid:63) ) ∞ (cid:88) n =0 e − β (cid:126) ωn n n ! H n ( αξ (cid:63) ) , (56)where H n ( x ) is the order- n Hermite polynomial at x , α = (cid:112) mω/ (cid:126) , and the partition function is Z = 12 csch( β (cid:126) ω/ . (57)Using the identity ∞ (cid:88) n =0 k n n n ! H n ( x ) H n ( y ) = e k x y − kxyk − √ − k (58)for | k | <
1, which in our case simplifies to ∞ (cid:88) n =0 k n n n ! H n ( x ) = e kx k +1 √ − k , (59)we find that (cid:37) ( ξ (cid:63) ) = (cid:114) α tanh( β (cid:126) ω/ π e − α tanh( β (cid:126) ω/ ξ (cid:63) ) . (60)Thus, − βA (cid:48) ( ξ (cid:63) ) = − α tanh( β (cid:126) ω/ ξ (cid:63) , (61)which is proportional to ξ (cid:63) . . . . . . . T / K − − − − − β A / n m − ξ ? = − ξ ? = 0 nm ξ ? = 4 nm Estimator 1Estimator 2
FIG. 1. Comparison of the estimators E and E in Eqs. (36)and (41) for the computation of the PMF derivative of aharmonic oscillator at ξ (cid:63) = − The necessary quantities for the PIMC estimators are ∂∂ξ log | J ( q ) | = 0 , (62a) ∂q∂ξ = 1 , (62b)and F ( q ) = − mω q. (62c)For this example, we have arbitrarily chosen m =1 . − and ω = 2 . − . The derivative of the PMFas computed using the Monte Carlo estimators agrees verywell with the exact result over a range of temperaturesand constraint positions, as shown in Fig. 1. B. Lennard-Jones dimer
To demonstrate that these estimators are applicableto a curvilinear reaction coordinate, we study a diatomicmolecule with reduced mass µ and Lennard-Jones inter-actions. Without the term for translation of the center ofmass, its Hamiltonian isˆ H = ˆ p q µ + V LJ ( ˆ ξ ) , (63)where q is the radial separation vector between the atoms,whose magnitude ξ = | q | we use as the reaction coordinate,and V LJ ( ξ ) = 4 ε (cid:34)(cid:18) σξ (cid:19) − (cid:18) σξ (cid:19) (cid:35) (64)is the Lennard-Jones potential. Unlike the harmonic oscil-lator example, this system has a potential that vanishesat large separation, allowing the dimer to dissociate. . . . . . . . ξ ? / nm − − − β A / n m − T = 20 K T = 4 K T = 2 K Estimator 1Estimator 2
FIG. 2. Comparison of the estimators E and E in Eqs. (36)and (41) for the computation of the PMF derivative of aLennard-Jones dimer at T = 20 K (top curve, least saturated),4 K (middle curve), and 2 K (bottom curve, most saturated).Error bars are not visible, because they are smaller than thesymbols. The solid curves show the NMM results. Expressing q in spherical coordinates ( ξ, cos θ, ϕ ), wehave that the magnitude of the Jacobian determinant is | J ( X , ξ ) | = ξ . (65)In order to evaluate the PIMC estimators, we thereforerequire the following quantities: ∂∂ξ log | J ( q ) | = 2 ξ , (66a) ∂ q ∂ξ = q ξ , (66b)and F ( q ) = 24 ε q ξ (cid:34) (cid:18) σξ (cid:19) − (cid:18) σξ (cid:19) (cid:35) . (66c)Note that we retain the geometric term during the simula-tion and explicitly remove it in the subsequent numericalintegration. To perform a reference calculation, we usenumerical matrix multiplication (NMM), as described inAppendix D.For the Lennard-Jones parameters provided in Ref. 29for Ar ( ε = 119 . σ = 3 . (cid:6) A), the results inFig. 2 confirm that the estimators E and E functioncorrectly with radial distance as a reaction coordinate. Inparticular, the rapid change in the slope of the PMF iscaptured at the lower temperatures.It is possible to numerically integrate the derivative A (cid:48) ( ξ (cid:63) ) to recover the PMF A ( ξ (cid:63) ). We do so using themidpoint rule on a grid of points ξ (cid:63)i with spacing ∆ ξ (cid:63) (as in Fig. 2), and with ξ (cid:63) placed at the largest value of ξ (cid:63) . We also define the virtual point ξ (cid:63) = ξ (cid:63) + ∆ ξ (cid:63) and ashifted grid of points¯ ξ (cid:63)i = ξ (cid:63)i − ∆ ξ (cid:63) , (67) . . . . . . . ξ ? / nm − . − . . . . ˜ A / k J m o l − T = 20 K T = 2 K Estimator 1Estimator 2
FIG. 3. Comparison of the estimators E and E in Eqs. (36)and (41) for the computation of the PMF of a Lennard-Jonesdimer at T = 20 K (narrow curve, least saturated) and 2 K(wide curve, most saturated). Error bars are not visible, be-cause they are smaller than the symbols. Additional pointsextending to ξ (cid:63) = 1 . with ¯ ξ (cid:63) acting as a “point at infinity” (the dimer is con-sidered to have dissociated when the atoms are at least¯ ξ (cid:63) apart). Correspondingly, we set ˜ A ( ¯ ξ (cid:63) ) = 0 = ˜ A (cid:48) ( ξ (cid:63) ),using the normalization in Eq. (45).In Fig. 3, we show˜ A ( ¯ ξ (cid:63)j ) = ∆ ξ (cid:63) β j (cid:88) i =1 (cid:20) − βA (cid:48) ( ξ (cid:63)i ) − ξ (cid:63)i (cid:21) , (68)which is the renormalized PMF with the desired energyoffset. The matching NMM curves are calculated from (cid:37) ( ξ (cid:63) ) as ˜ A ( ξ (cid:63) ) = − β log (cid:37) ( ξ (cid:63) )( ¯ ξ (cid:63) ) (cid:37) ( ¯ ξ (cid:63) )( ξ (cid:63) ) (69)to ensure a compatible energy offset. Even though theintegration grid is rather sparse, especially where the slopeof the PMF changes suddenly for T = 2 K, the obtainedPMFs are consistent with the reference results. V. CONCLUSIONS
We have described two path integral Monte Carlo esti-mators for the calculation of the derivative of the quantummechanical potential of mean force. Notably, the curvesobtained from these estimators are in terms of the truequantum reaction coordinate, unlike other methods thatutilize the path centroid.The first estimator, Eq. (36), was obtained by initiallydifferentiating the exact path integral and then discretiz-ing the resulting path integral into imaginary time steps.Alternatively, the second estimator, Eq. (41), was ob-tained by discretizing the exact path integral first and then performing the differentiation after. In principle,these should be equivalent operations in the P → ∞ limit, and we have demonstrated that both estimatorsreproduce the correct derivative of the PMF for the one-dimensional harmonic oscillator and Lennard-Jones dimer.In contrast to existing histogram-based methods for theevaluation of free energies, these novel estimators can beused to ascertain information about the free energy profileat just a single point along the reaction coordinate.Furthermore, it is possible to numerically integrate thecomputed derivatives evaluated from these estimators toobtain the PMF itself. As shown in the argon dimer exam-ple, even when the integration grid is not very dense, thismethod successfully reproduces the known PMF obtainedfrom numerical matrix multiplication.In Paper II of this series, we show how these estima-tors may be used with path integral molecular dynamics.This is achieved by applying techniques from constrainedLangevin dynamics to the PILE integrator in order toconstrain one of the beads. The extension of these esti-mators to path integral molecular dynamics simulationswill allow for their application to more general systemsand potentials, such as small water clusters. ACKNOWLEDGMENTS
We thank Raymond Kapral for providing the initialdirection for the derivation of the discretized constrainedpath integral. This research was supported by the NaturalSciences and Engineering Research Council of Canada(NSERC) (RGPIN-2016-04403), the Ontario Ministry ofResearch and Innovation (MRI), the Canada ResearchChair program (950-231024), and the Canada Foundationfor Innovation (CFI) (project No. 35232).
Appendix A: Kets in curvilinear coordinates
A wavefunction ψ ( q ) may be thought of as the con-crete manifestation of an abstract ket | ψ (cid:105) in a continuousrepresentation: ψ ( q ) = (cid:104) q | ψ (cid:105) . (A1)Although the object | q (cid:105) (which represents a state withdefinite Cartesian position q ) is not an element of Hilbertspace, it is common to formally treat it as if it were.Given a change of variables from q to X , ξ with Jacobiandeterminant J ( q ) = J ( X , ξ ), it is useful to define | X ξ (cid:105) in a way that fulfills (cid:90) d q |(cid:104) q | ψ (cid:105)| = (cid:90) d X (cid:90) d ξ |(cid:104) X ξ | ψ (cid:105)| , (A2)which is analogous to the statement that the resolutionof the identity ˆ = (cid:90) d X (cid:90) d ξ | X ξ (cid:105)(cid:104) X ξ | (A3)should have the usual form, even in curvilinear coordi-nates.Since (cid:90) d q |(cid:104) q | ψ (cid:105)| = (cid:90) d X (cid:90) d ξ | J ( X , ξ ) | |(cid:104) q ( X , ξ ) | ψ (cid:105)| , (A4)it follows that the definition | X ξ (cid:105) = (cid:112) | J ( X , ξ ) | | q ( X , ξ ) (cid:105) = (cid:112) | J ( q ) | | q (cid:105) (A5)is sufficient. This is the approach described in Ref. 25, andthe one we use in the present work. Using this definition,we see that the diagonal matrix elements of the partialtrace of an operator ˆ O with respect to X may be expressedas (cid:104) ξ (cid:63) | Tr X ˆ O | ξ (cid:63) (cid:105) = (cid:90) d X (cid:104) X ξ (cid:63) | ˆ O | X ξ (cid:63) (cid:105) (A6a)= (cid:90) d X (cid:90) d ξ δ ( ξ − ξ (cid:63) ) (cid:104) X ξ | ˆ O | X ξ (cid:105) (A6b)= (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) X ξ | ˆ O | X ξ (cid:105)| J ( q ) | (A6c)= (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:104) q | ˆ O | q (cid:105) (A6d)in Cartesian coordinates. Appendix B: Derivative of a Dirac delta function integral
We wish to take the derivative D ( ξ (cid:63) ) = dd ξ (cid:63) (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) f ( q ) . (B1)We first obtain the one-dimensional resultdd ξ (cid:63) (cid:90) d ξ δ ( ξ − ξ (cid:63) ) f ( ξ ) = (cid:90) d ξ δ ( ξ − ξ (cid:63) ) dd ξ f ( ξ ) (B2)by noting thatdd ξ (cid:63) f ( ξ (cid:63) ) = (cid:90) d ξ δ ( ξ − ξ (cid:63) ) dd ξ f ( ξ ) . (B3)For the general case, we change coordinates to those inwhich ξ appears explicitly: D ( ξ (cid:63) ) = (cid:90) d X dd ξ (cid:63) (cid:90) d ξ δ ( ξ − ξ (cid:63) ) | J ( X , ξ ) | f ( X , ξ )(B4a)= (cid:90) d X (cid:90) d ξ δ ( ξ − ξ (cid:63) ) ∂∂ξ | J ( X , ξ ) | f ( X , ξ ) (B4b)= (cid:90) d X (cid:90) d ξ δ ( ξ − ξ (cid:63) ) (cid:20) ∂∂ξ | J ( X , ξ ) | (cid:21) f ( X , ξ )+ (cid:90) d X (cid:90) d ξ δ ( ξ − ξ (cid:63) ) | J ( X , ξ ) | ∂∂ξ f ( X , ξ )(B4c)= (cid:90) d q δ ( ξ ( q ) − ξ (cid:63) ) (cid:20) J ξ ( q ) + ∂∂ξ (cid:21) f ( q ) , (B4d) where J ξ ( q ) = ∂∂ξ log | J ( q ) | = ∂∂ξ | J ( X , ξ ) || J ( X , ξ ) | , (B5)and we formally apply the logarithmic derivative notationeven when the function is not dimensionless. Appendix C: Derivative–commutator identity for diagonalmatrix elements
It is well-known that momentum operators lead to dif-ferentiation in the position representation. For example, (cid:104) q | ˆ p i ˆ A | q (cid:48) (cid:105) = − i (cid:126) ∂∂q i (cid:104) q | ˆ A | q (cid:48) (cid:105) (C1)for an arbitrary operator ˆ A , where ˆ p i is the momentumoperator conjugate to ˆ q i . However, this relationship doesnot generally hold when q and q (cid:48) are the same variable: (cid:104) q | ˆ p i ˆ A | q (cid:105) (cid:54) = − i (cid:126) ∂∂q i (cid:104) q | ˆ A | q (cid:105) . (C2)Instead, for a Hermitian operator ˆ A with the eigenvalueequation ˆ A | a (cid:105) = a | a (cid:105) , we have that (cid:104) q | ˆ p i ˆ A | q (cid:105) = (cid:88) a (cid:104) q | ˆ p i | a (cid:105) (cid:104) a | ˆ A | q (cid:105) (C3a)= − i (cid:126) (cid:88) a a (cid:20) ∂∂q i (cid:104) q | a (cid:105) (cid:21) (cid:104) a | q (cid:105) (C3b)and (cid:104) q | ˆ A ˆ p i | q (cid:105) = i (cid:126) (cid:88) a a (cid:104) q | a (cid:105) (cid:20) ∂∂q i (cid:104) a | q (cid:105) (cid:21) . (C4)Thus, we conclude that ∂∂q i (cid:104) q | ˆ A | q (cid:105) = ∂∂q i (cid:88) a (cid:104) q | ˆ A | a (cid:105) (cid:104) a | q (cid:105) (C5a)= (cid:88) a a (cid:104) q | a (cid:105) (cid:20) ∂∂q i (cid:104) a | q (cid:105) (cid:21) + (cid:88) a a (cid:20) ∂∂q i (cid:104) q | a (cid:105) (cid:21) (cid:104) a | q (cid:105) (C5b)= 1 i (cid:126) (cid:104) q | ˆ A ˆ p i | q (cid:105) − i (cid:126) (cid:104) q | ˆ p i ˆ A | q (cid:105) (C5c)= 1 i (cid:126) (cid:104) q | (cid:104) ˆ A, ˆ p i (cid:105) | q (cid:105) . (C5d) Appendix D: Numerical matrix multiplication for a radialcoordinate
In Ref. 29, expressions for numerical matrix multipli-cation of the path integral of a system described by athree-dimensional relative coordinate are given, but notderived. In this section, we briefly explain why the radialpropagator has such a curious form.The operator in the kinetic energy of Eq. (63) may beexpressed as ˆ p q = ˆ p ξ + ˆ (cid:96) ˆ ξ , (D1)where ˆ p ξ is the radial momentum operator, and ˆ (cid:96) is thesquared angular momentum operator, whose eigenstatesare the spherical harmonics | (cid:96) m (cid:105) with eigenvalues (cid:126) (cid:96) ( (cid:96) +1). The radial momentum operator is not self-adjointand does not have a spectrum of eigenstates, so thespectral theorem does not apply to it and the appropriateresolution of identity is not given by (cid:90) d p ξ | p ξ (cid:105)(cid:104) p ξ | , (D2) despite the wavefunctions (cid:104) ξ | p ξ (cid:105) = e iξpξ (cid:126) √ π (cid:126) (D3)satisfying ˆ p ξ | p ξ (cid:105) = p ξ | p ξ (cid:105) . Thus, we must be carefulwhen rederiving Eq. (16c) of Ref. 29.We turn to the operator ˆ p ξ , which is well-behaved andhas the eigenstates (cid:104) ξ | p (2) ξ (cid:105) = e iξpξ (cid:126) − e − iξpξ (cid:126) i √ π (cid:126) = 1 √ π (cid:126) sin ξp ξ (cid:126) (D4)with eigenvalues p ξ . We may use these states to constructthe correct resolution of the identity,ˆ = (cid:90) d p ξ | p (2) ξ (cid:105)(cid:104) p (2) ξ | = (cid:90) d p ξ (cid:16) | p ξ (cid:105)(cid:104) p ξ | − | p ξ (cid:105)(cid:104)− p ξ | (cid:17) , (D5)which, as expected, results in (cid:104) ξ (cid:48) | e − τ ˆ p ξ µ | ξ (cid:105) = (cid:90) d p ξ (cid:104) ξ (cid:48) | e − τp ξ µ | p (2) ξ (cid:105) (cid:104) p (2) ξ | ξ (cid:105) (D6a)= 14 π (cid:126) (cid:90) d p ξ e − τp ξ µ + ipξ (cid:126) ( ξ (cid:48) − ξ ) + 14 π (cid:126) (cid:90) d p ξ e − τp ξ µ − ipξ (cid:126) ( ξ (cid:48) − ξ ) − π (cid:126) (cid:90) d p ξ e − τp ξ µ + ipξ (cid:126) ( ξ (cid:48) + ξ ) − π (cid:126) (cid:90) d p ξ e − τp ξ µ − ipξ (cid:126) ( ξ (cid:48) + ξ ) (D6b)= (cid:114) µ π (cid:126) τ (cid:104) e − µ (cid:126) τ ( ξ (cid:48) − ξ ) − e − µ (cid:126) τ ( ξ (cid:48) + ξ ) (cid:105) . (D6c) G. M. Torrie and J. P. Valleau, J. Comp. Phys. , 187 (1977). S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen, andP. A. Kollman, J. Comp. Chem. , 1011 (1992). B. Roux, Comp. Phys. Comm. , 275 (1995). E. Carter, G. Ciccotti, J. T. Hynes, and R. Kapral, Chem. Phys.Lett. , 472 (1989). M. Sprik and G. Ciccotti, J. Chem. Phys. , 7737 (1998). G. Ciccotti, R. Kapral, and E. Vanden-Eijnden, ChemPhysChem , 1809 (2005). A. Laio and M. Parrinello, Proc. Natl. Acad. Sci. U.S.A. , 12562(2002). G. Bussi, A. Laio, and M. Parrinello, Phys. Rev. Lett. , 090601(2006). T. M¨ulders, P. Kr¨uger, W. Swegat, and J. Schlitter, J. Chem.Phys. , 4869 (1996). W. K. den Otter and W. J. Briels, Mol. Phys. , 773 (2000). M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie,A. Michaelides, M. A. Morales, and T. E. Markland, Chem.Rev. , 7529 (2016). T. E. Markland and M. Ceriotti, Nat. Rev. Chem. , 0109 (2018). C. L. Vaillant, S. C. Althorpe, and D. J. Wales, J. Chem. TheoryComput. , 33 (2018). K. P. Bishop and P.-N. Roy, J. Chem. Phys. , 102303 (2018). E. M´endez and D. Laria, J. Chem. Phys. , 054302 (2020). K. Hinsen and B. Roux, J. Chem. Phys. , 3567 (1997). D. T. Major, M. Garcia-Viloca, and J. Gao, J. Chem. TheoryComput. , 236 (2006). D. T. Major and J. Gao, J. Chem. Theory Comput. , 949 (2007). B. Walker and A. Michaelides, J. Chem. Phys. , 174306 (2010). A. Vardi-Kilshtain, D. T. Major, A. Kohen, H. Engel, andD. Doron, J. Chem. Theory Comput. , 4786 (2012). J. Cao and G. A. Voth, J. Chem. Phys. , 5093 (1994). J. Cao and G. A. Voth, J. Chem. Phys. , 6168 (1994). N. Blinov and P.-N. Roy, J. Chem. Phys. , 3759 (2004). M. Ceriotti, M. Parrinello, T. E. Markland, and D. E. Manolopou-los, J. Chem. Phys. , 124104 (2010). B. Leaf, Found. Phys. , 581 (1980). R. Wilcox, J. Math. Phys. , 962 (1967). D. Iouchtchenko, “SnaggedNecklace,” https://github.com/0/SnaggedNecklace.jl , accessed 2021-01-03. Wolfram Research Inc., “Hermite polynomials: Summation(formula 05.01.23.0013),” http://functions.wolfram.com/05.01.23.0013.01 , accessed 2020-01-20. D. Thirumalai, E. J. Bruskin, and B. J. Berne, J. Chem. Phys. , 5063 (1983). R. L. Liboff, I. Nebenzahl, and H. H. Fleischmann, Am. J. Phys. , 976 (1973). G. Paz, J. Phys. A35