How to remove the spurious resonances from ring polymer molecular dynamics
HHow to remove the spurious resonances from ring polymer molecular dynamics
Mariana Rossi, Michele Ceriotti, and David E. Manolopoulos Physical and Theoretical Chemistry Laboratory, University of Oxford,South Parks Road, Oxford OX1 3QZ, United Kingdom Laboratory of Computational Science and Modeling, IMX,´Ecole Polytechnique F´ed´erale de Lausanne, 1015 Lausanne, Switzerland
Two of the most successful methods that are presently available for simulating the quantumdynamics of condensed phase systems are centroid molecular dynamics (CMD) and ring polymermolecular dynamics (RPMD). Despite their conceptual differences, practical implementations ofthese methods differ in just two respects: the choice of the Parrinello-Rahman mass matrix andwhether or not a thermostat is applied to the internal modes of the ring polymer during the dynamics.Here we explore a method which is halfway between the two approximations: we keep the pathintegral bead masses equal to the physical particle masses but attach a Langevin thermostat to theinternal modes of the ring polymer during the dynamics. We justify this by showing analytically thatthe inclusion of an internal mode thermostat does not affect any of the desirable features of RPMD:thermostatted RPMD (TRPMD) is equally valid with respect to everything that has actually beenproven about the method as RPMD itself. In particular, because of the choice of bead masses, theresulting method is still optimum in the short-time limit, and the transition state approximationto its reaction rate theory remains closely related to the semiclassical instanton approximation inthe deep quantum tunneling regime. In effect, there is a continuous family of methods with theseproperties, parameterised by the strength of the Langevin friction. Here we explore numericallyhow the approximation to quantum dynamics depends on this friction, with a particular emphasison vibrational spectroscopy. We find that a broad range of frictions approaching optimal dampinggive similar results, and that these results are immune to both the resonance problem of RPMDand the curvature problem of CMD.
I. INTRODUCTION
Centroid molecular dynamics (CMD) and ring poly-mer molecular dynamics (RPMD) are two closely re-lated approximate techniques for including quantum me-chanical zero point energy and tunnelling effects in molec-ular dynamics simulations of condensed phase systems.Both are based on the imaginary time path integral for-mulation of quantum statistical mechanics, but in differ-ent ways: RPMD is classical molecular dynamics in theextended phase space of the imaginary time path integral(or ring polymer ), whereas CMD is classical moleculardynamics on the potential of mean force generated bythe thermal fluctuations of this ring polymer around itscentroid.Both of these techniques are now well established andhave been used in a wide variety of applications. Theyare especially useful for calculating diffusion coefficientsand other transport properties of liquids, and for in-cluding zero-point energy and tunnelling effects in thecalculation of chemical reaction rates.
In these con-texts, both methods appear to be reasonably reliable, onthe basis of comparisons with full quantum mechanicalresults in cases where these are available, andexact quantum mechanical moment constraints in caseswhere they are not.
One explanation for whyRPMD works well for these problems is that diffusioncoefficients and reaction rate coefficients are determinedby the zero-frequency components of the relevant (veloc-ity and reactive flux) autocorrelation spectra, and aretherefore comparatively insensitive to the artificial high- frequency dynamics of the ring polymer internal modes. There are however some important situations in whichRPMD and CMD are known to give unphysical results,the classic example being in the simulation of vibrationalspectroscopy. Here RPMD fails because of the so-called resonance problem : when the frequency of a phys-ical vibration (such as an O-H stretch in liquid water)comes into resonance with an excited internal mode ofthe ring polymer in another degree of freedom (such asa librational mode of the liquid), the resulting interac-tion causes a spurious splitting of the physical peak inthe calculated spectrum. In a simulation at room tem-perature, the first internal excitation of the ring polymeroccurs at a wavenumber of kT / (cid:126) c ∼ − , and soany spectral feature beyond this wavenumber could inprinciple be affected by this problem.The resonance problem does not occur in CMD be-cause the high-frequency vibrations of the internal modesof the ring polymer are averaged over in this method tocalculate the centroid potential of mean force: they donot appear as dynamical variables. However, CMD canalso give unphysical vibrational spectra, for a differentreason known as the curvature problem . This prob-lem is again inherently multidimensional and is expectedto arise whenever a high frequency stretching coordinateis combined with an angular coordinate, such as a rota-tional or a torsional degree of freedom. At low tempera-tures, the ring polymer spreads out around the angularcoordinate, resulting in a softening of the centroid poten-tial of mean force along the stretch.
Because of this,the stretching peak in the simulated spectrum is artifi-cially red-shifted and broadened. In the extreme case of a r X i v : . [ phy s i c s . c h e m - ph ] J un a combination with a free rotation, the frequency of thestretching peak can even tend to zero as T → In this paper, we shall show that both of these prob-lems can be avoided by adopting a method that is halfwaybetween RPMD and CMD. Despite their conceptual dif-ferences, practical implementations of these two methodsare very similar. In RPMD, the masses of the ring poly-mer beads are chosen to be the physical particle masses,and the dynamics that is used to calculate correlationfunctions is microcanonical. In the adiabatic implemen-tation of CMD, the internal modes of the ring polymerare given much smaller masses so that they vibrate atfrequencies well beyond the physical spectral range, anda thermostat is attached to these internal modes to gen-erate the canonical centroid potential of mean force. Weshall show here that, in the context of vibrational spec-troscopy, attaching a thermostat to the internal modes ofthe ring polymer during the dynamics is a good idea, be-cause it damps out the spurious resonances from RPMD.However, adiabatically separating the internal vibrationsof the ring polymer from the centroid motion is a notsuch a good idea, because it is the origin of the curvatureproblem in CMD.The reason why it is legitimate to attach a thermo-stat to the internal modes of the ring polymer duringthe dynamics is discussed in Sec. II: provided the ther-mostat only acts on the internal modes and not on thecentroid, it does not have any impact on any of the estab-lished properties of RPMD. Sec. III goes on to presentsome numerical examples to illustrate how attaching athermostat to the internal modes while keeping the beadmasses equal to the physical masses cures both the res-onance problem of RPMD and the curvature problemof CMD, for a variety of model and more realistic vi-brational problems. It also happens to be the case thatthermostatted RPMD (TRPMD) provides the simplestcomputational scheme, because it avoids both the non-ergodicity of the microcanonical dynamics in RPMD andthe need to use a small integration time step in CMD.Finally, to emphasise that it is not restricted to the calcu-lation of vibrational spectra, we show that TRPMD givesa diffusion coefficient and orientational relaxation timesfor liquid water that are very similar to those obtainedfrom RPMD and CMD. Sec. IV contains some concludingremarks. II. SOME PROPERTIES OF THERMOSTATTEDRPMDA. Standard ring polymer molecular dynamics
Although RPMD can also be used to approximatecross correlation functions, we shall confine our atten-tion here to Kubo-transformed autocorrelation functionsof the form ˜ c AA ( t ) = 1 βZ (cid:90) β dλ tr (cid:104) e − ( β − λ ) ˆ H ˆ A e − λ ˆ H ˆ A ( t ) (cid:105) , (1) where Z = tr (cid:104) e − β ˆ H (cid:105) , (2)and ˆ A ( t ) = e + i ˆ Ht/ (cid:126) ˆ A e − i ˆ Ht/ (cid:126) . (3)Correlation functions of this form arise in a wide varietyof contexts, including the theory of diffusion (in whichˆ A is the operator for the velocity of a tagged particlein a liquid), the theory of reaction rates (in which ˆ A isthe operator for the flux of particles through a transi-tion state dividing surface), and, most relevantly to thepresent study, the theory of vibrational spectroscopy (inwhich ˆ A is the dipole moment operator of the entire sys-tem).Since the multi-dimensional generalisation isstraightforward, and it does not add anything otherthan more indices, we shall present all of the followingequations for a model one-dimensional system with aHamiltonian of the formˆ H = ˆ p m + V (ˆ q ) . (4)Furthermore, since most physically interesting correla-tion functions involving non-local operators can be ob-tained by differentiating correlation functions of local op-erators (for example, the velocity autocorrelation func-tion of a tagged particle in a liquid is minus the secondtime derivative of its position autocorrelation function),we shall restrict our attention to the case where ˆ A is alocal Hermitian operator A (ˆ q ).The n -bead RPMD approximation to the correlationfunction in Eq. (1) is then simply ˜ c AA ( t ) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) A n ( q ) A n ( q t ) , (5)where N n = (2 π (cid:126) ) n Z n with Z n = 1(2 π (cid:126) ) n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) . (6)Here H n ( p , q ) is the classical Hamiltonian of a harmonicring polymer with a potential of V ( q ) acting on eachbead, H n ( p , q ) = n (cid:88) j =1 (cid:34) p j m + 12 mω n ( q j − q j − ) + V ( q j ) (cid:35) , (7)with ω n = 1 / ( β n (cid:126) ), β n = β/n , and q ≡ q n . The coor-dinates q t = q t ( p , q ) in Eq. (5) are obtained from theclassical dynamics generated by this Hamiltonian,˙ p = − ∂H n ( p , q ) ∂ q , (8)˙ q = + ∂H n ( p , q ) ∂ p , (9)and the functions A n ( q ) and A n ( q t ) are averaged overthe beads of the ring polymer necklace at times 0 and t , A n ( q ) = 1 n n (cid:88) j =1 A ( q j ) . (10)The approximation in Eq. (5) is clearly ad hoc , and noone has yet been able to derive it from first principles.However, it is consistent with the path integral descrip-tion of static equilibrium properties, it gives a correlationfunction with the same time-reversal and detailed bal-ance symmetries as the quantum correlation function (which in the present context combine to make ˜ c AA ( t ) areal and even function of t ), it is exact for problems withcontinuous spectra in the high-temperature limit wherethe ring polymer shrinks to a single classical bead, it is asaccurate as possible in the short-time limit for a methodbased on the classical evolution of imaginary-time pathintegral variables, and it is exact at all times for theposition autocorrelation function of a simple harmonicoscillator. Furthermore, Richardson and Althorpe have recentlyestablished a close connection between RPMD transitionstate theory (TST) and the semiclassical instanton ap-proximation in the deep quantum tunnelling regime, andHele and Althorpe have gone on to show that RPMDTST is in fact the only known quantum mechanical ver-sion of TST with positive definite statistics that can bederived from the t → + limit of a reactive flux-side cor-relation function. These important findings, which legit-imise the use of RPMD for calculating chemical reactionrates more so than in any other context, are both clearlyassociated with the accuracy of the RPMD approxima-tion in the short-time limit.The key observation that motivated the present studyis that, as we shall show next, all of the above featuresof RPMD continue to hold when the Hamiltonian dy-namics of the ring polymer necklace is augmented with aLangevin thermostat, provided the thermostat only actson the internal modes of the ring polymer and not onits centroid. In effect, RPMD is just one of a familyof dynamical approximations that share the same cor-rect symmetries and limiting cases, the family being pa-rameterised by the strength of the Langevin friction. Soby exploring the whole family, one might hope to find amethod that alleviates some of the known problems ofRPMD (such as the resonance problem mentioned in theIntroduction), without sacrificing any of its establishedfeatures. B. The effect of an internal mode thermostat
The Hamiltonian equations of motion in Eqs. (8)and (9) can be written out more explicitly as˙ p = F ( q ) = f ( q ) − Kq , (11) ˙ q = 1 m p . (12)Here K is a real, symmetric, positive semi-definite n × n spring constant matrix with elements K jj (cid:48) = mω n (cid:0) δ jj (cid:48) − δ jj (cid:48) − − δ jj (cid:48) +1 (cid:1) , (13)all indices being understood to be translated by a multi-ple of n to lie between 1 and n , and f ( q ) is a force vectorwith elements f j ( q ) = − V (cid:48) ( q j ).An important feature of the spring constant matrix isthat it does not have any impact on the centroid¯ p = e T p , ¯ q = e T q (14)of the ring polymer, but only on the ring polymer internalmodes. To see this, it suffices to note that e T K = T and K e = , (15)where e is a vector with elements e j = 1 /n and is avector of zeroes. This feature was used repeatedly byBraams and Manolopoulos to investigate the short-timelimit of RPMD, and we shall use it again in this contextbelow.Using the same notation, suppose that we now at-tach a Langevin thermostat to the dynamics by replacingEq. (11) with˙ p = F ( q ) − γ p + (cid:114) m γ β n ξ ( t ) , (16)where γ is a real, symmetric, positive semi-definite n × n friction matrix and ξ ( t ) is a vector of uncorrelated normaldeviates with unit variance and zero mean ( (cid:104) ξ j ( t ) (cid:105) = 0and (cid:104) ξ j (0) ξ j (cid:48) ( t ) (cid:105) = δ jj (cid:48) δ ( t )). Then the thermostat willalso be detached from the centroid (i.e., only act directlyon the ring polymer internal modes) if e is an eigenvectorof γ (and therefore also of √ γ ) with eigenvalue zero, e T γ = T and γ e = . (17)So let us consider whether any of the established prop-erties of RPMD are affected by replacing Eq. (11) withEq. (16) subject to this constraint on γ .First of all (and perhaps most importantly), the con-sistency with imaginary time path integral expressionsfor static equilibrium properties is maintained. One caneither calculate these properties by averaging over shortRPMD trajectories with initial conditions sampled fromthe Boltzmann distribution e − β n H n ( p , q ) /N n , or by at-taching a thermostat to the dynamics as in Eq. (16). Infact, for the calculation of static equilibrium properties,attaching a thermostat to the dynamics is widely recog-nised as the better way to proceed.The symmetry properties of RPMD are also unaffectedby replacing Eq. (11) with Eq. (16), at least in the casethat we are considering here: ˜ c AA ( t ) is still a real andeven function of t . This follows because the Langevindynamics is a stationary process – a process that con-serves the equilibrium distribution and is independent ofthe origin of time – despite the appearance of t in thenoise term on the right-hand side of Eq. (16). (The noiseis uncorrelated from one instant of time to the next, andit has the same form at all instants of time.)In more detail, the autocorrelation function that re-sults from the Langevin dynamics is˜ c AA ( t ) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) A n ( q ) ¯ A n ( q t ) , (18)where the overbar denotes an average over the stochasticprocess that takes the integration variables p and q to q t ( p , q ) in time t > c AA ( t ) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) ¯ A n ( q t ) ¯ A n ( q t + t ) , (19)or˜ c AA ( t ) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) ¯ A n ( q t − t ) ¯ A n ( q t ) , (20)for any time origin t > t , and comparing these two equa-tions we see that ˜ c AA ( t ) = ˜ c AA ( − t ). Since A (ˆ q ) is Her-mitian, ˜ c AA ( t ) is also real, and therefore a real and evenfunction of t .So far, we have not made any use of the constraint inEq. (17), which detaches the thermostat from the cen-troid of the ring polymer. However, this constraint isclearly needed for the dynamics to give the correct resultin the high-temperature limit, and also for the positionautocorrelation function of a harmonic oscillator.In the high temperature limit, where a single bead suf-fices to converge the autocorrelation function, the con-straint makes γ zero, so the thermostat has no effect.One therefore recovers standard RPMD in this limit,which gives the correct (classical) result for systems withcontinuous spectra. As for the harmonic limit, whenEq. (17) is satisfied the only term in Eq. (16) that couplesthe centroid to the internal modes of the ring polymer isthe force vector f ( q ): in general, e T f ( q ) is not just afunction of the centroid coordinate ¯ q = e T q . However,when the potential V ( q ) is harmonic, V ( q ) = mω q ,one has e T f ( q ) = − mω ¯ q , and the coupling of the cen-troid to the internal modes disappears. With the con-straint in Eq. (17), thermostatted RPMD will thereforegive the same position (¯ q ) autocorrelation function for aharmonic oscillator as RPMD.Finally, let us consider the effect of the thermostat onthe short-time limit of the correlation function ˜ c AA ( t ).Braams and Manolopoulos have given a detailed analy-sis of the short-time limit of the standard RPMD approx-imation in Eq. (5), in which they argued that it has anerror of O ( t ) for the position autocorrelation functionand of O ( t ) for a more general autocorrelation function of a non-linear operator A (ˆ q ). However, this argumentassumed the odd time derivatives of the RPMD auto-correlation function to be continuous at t = 0, as theyare in quantum mechanics. By considering the specialcase of a simple harmonic oscillator, for which the au-tocorrelation function can be worked out analytically, Jang et al. have shown that the third derivative of theRPMD approximation to ˜ c q q ( t ) becomes discontinuousat t = 0 in the limit as n → ∞ (the correlation functioncontains a term proportional to | t | in this limit ). Thisimplies that the error in the RPMD approximation to thisnon-linear autocorrelation function is O ( t ) rather than O ( t ). Moreover a simple argument shows that the er-ror in the RPMD position autocorrelation function willtherefore be O ( t ) rather than O ( t ). This follows be-cause ˜ c (4) qq ( t ) = ˜ c ff ( t ), where f ( q ) = − V (cid:48) ( q ) is the force.In an anharmonic potential, f ( q ) is a non-linear functionof q (for example it could contain a component propor-tional to q ), so the RPMD approximation to ˜ c ff ( t ) hasan error of O ( t ), and the approximation to ˜ c qq ( t ) anerror of O ( t ).To see how thermostatting the internal modes of thering polymer affects these results, we shall now com-pare the first few initial time derivatives of the Langevin-thermostatted correlation function in Eq. (18) with thoseof the standard RPMD correlation function in Eq. (5).The k -th (forward) initial time derivative of ˜ c AA ( t ) isgiven by Eq. (18) as˜ c ( k ) AA (0) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) A n ( q ) ¯ A ( k ) n ( q ) , (21)where ¯ A ( k ) n ( q ) is the k -th time derivative of the noise-averaged ¯ A n ( q t ) at t = 0. This is given by ¯ A ( k ) n ( q ) = (cid:0) D † (cid:1) k A n ( q ) , (22)where D † is the adjoint of the Fokker-Planck opera-tor associated with the Langevin dynamics in Eqs. (12)and (16), D † = p m · ∂∂ q + F ( q ) · ∂∂ p − p · γ · ∂∂ p + mβ n ∂∂ p · γ · ∂∂ p . (23)The time derivatives of the standard RPMD autocorre-lation are recovered by setting γ = in this operator,which converts it into the adjoint of the Liouville operatorassociated with the Hamiltonian dynamics in Eqs. (11)and (12). The goal is therefore to find the smallest valueof k for which ˜ c ( k ) AA (0) involves any reference to γ , whichsuffices to establish that the thermostatted and standardRPMD autocorrelation functions differ by O ( t k ).In the special case of the position autocorrelation func-tion, for which A (ˆ q ) = ˆ q and A n ( q ) = e T q , a straight-forward but lengthy calculation that combines Eqs. (21)to (23) with the simplifications in Eqs. (15) and (17) gives˜ c (0) qq (0) = (cid:10) e T q q T e (cid:11) , (24 a )˜ c (1) qq (0) = 0 , (24 b )˜ c (2) qq (0) = − β n m (cid:10) e T e (cid:11) ≡ − βm , (24 c )˜ c (3) qq (0) = 0 , (24 d )˜ c (4) qq (0) = + 1 β n m (cid:10) e T H ( q ) e (cid:11) , (24 e )˜ c (5) qq (0) = 0 , (24 f )˜ c (6) qq (0) = − β n m (cid:10) e T H ( q ) e (cid:11) , (24 g )˜ c (7) qq (0) = + 1 β n m (cid:10) e T H ( q ) γ H ( q ) e (cid:11) . (24 h )Here the angular brackets denote an equilibrium average (cid:104) F ( q ) (cid:105) = 1 N n (cid:90) d p (cid:90) d q e − β n H n ( p , q ) F ( q ) , (25)and H ( q ) is a diagonal matrix with diagonal elements H jj ( q ) = V (cid:48)(cid:48) ( q j ).When the potential V ( q ) is harmonic, e T H ( q ) γ H ( q ) e is proportional to e T γ e , which vanishes by virtue ofEq. (17). All subsequent occurrences of γ in the higherderivatives ˜ c ( k ) qq (0) also vanish for the same reason, sothe thermostat has no effect on the position autocorre-lation function. But when the potential V ( q ) is anhar-monic, there is no such simplification, and it follows fromEqs. (24) that the difference between the thermostattedand standard RPMD autocorrelation functions is O ( t ).Since this is exactly the same order as the error in theRPMD position autocorrelation function (see above), theshort-time accuracy of the approximation is unaffectedby attaching a thermostat to the internal modes of thering polymer during the dynamics.By contrast, changing the Parrinello-Rahman massmatrix so as to adiabatically separate the internal modesfrom the centroid mode, as is done in the adiabatic imple-mentation of CMD, certainly does have an effect on theshort-time accuracy. The leading error is then O ( t ) evenbefore a thermostat is attached to the internal modes, and the leading error in the fully converged CMD posi-tion autocorrelation function in which the centroid movesclassically on its potential of mean force is only O ( t ). In the more general case of the autocorrelation func-tion of a non-linear operator A (ˆ q ), RPMD is considerablyless accurate in the short-time limit, and the same istrue of thermostatted RPMD. In this case one finds fromEqs. (21) to (23) that the first time derivative ˜ c ( k ) AA (0) toinvolve γ is˜ c (3) AA (0) = + 1 β n m (cid:10) e T A (cid:48) ( q ) γ A (cid:48) ( q ) e (cid:11) , (26) where A (cid:48) ( q ) is a diagonal matrix with diagonal elements A (cid:48) ( q j ). If A (ˆ q ) were a linear operator, these diagonal el-ements would all be the same, and ˜ c (3) AA (0) would be zeroby virtue of Eq. (17). But for a non-linear operator, theright-hand side of Eq. (26) is non-zero, so TRPMD dif-fers from RPMD by a term of O ( t ). Again, since this isthe same order as the error in the RPMD autocorrelationfunction, nothing is lost by attaching a thermostat tothe internal modes during the dynamics. (CMD is not re-ally applicable to autocorrelation functions of non-linearoperators, but it is nevertheless often used for this pur-pose in the form of the “CMD with classical operators”approximation. The leading error is then O (1) – again3 orders less accurate than RPMD.)A closely related question is whether attaching aLangevin thermostat to the internal modes of the ringpolymer will affect any of the results of Richardson andAlthorpe or Hele and Althorpe concerning RPMDTST. The RPMD TST rate is obtained from the t → + limit of a flux-side correlation function, which is (minus)the first time derivative of a side-side correlation func-tion. The question is therefore whether the thermostatwill affect the short-time limit of this side-side correla-tion function to O ( t ). This can be investigated usingEqs. (21) to (23) by replacing A n ( q ) in Eq. (10) with A n ( q ) = h [ f ( q )], where h ( x ) is the Heaviside step func-tion and f ( q ) = 0 is a transition state dividing surfacein the ring polymer configuration space. Since thisgives ¯ A (1) n ( q ) = δ [ f ( q )] ∂f ( q ) /∂ q · p /m , which does notinvolve the friction matrix γ , it is clear that attaching athermostat to the dynamics will not have any effect onRPMD TST. C. Specification of the friction matrix
So far, our discussion of TRPMD has been fairly gen-eral: we have not specified the friction matrix γ inEq. (16) beyond saying that it is real, symmetric, andpositive semi-definite and that it satisfies the constraintin Eq. (17). One could go further by considering a gen-eralised Langevin equation thermostat in the form of aLangevin thermostat in an extended momentum space, which has the potential to provide an even milder pertur-bation to the Hamiltonian ring polymer dynamics. Butthat would lengthen the discussion considerably, and wehave not yet explored this option in any detail. So in-stead, to make things more concrete, we shall now sim-ply describe the particular family of friction matrices wehave used for illustrative purposes in the calculations pre-sented in Sec. III.These friction matrices are based on the path inte-gral Langevin equation (PILE) thermostat introduced inRef. . The basic idea is to transform to the normal moderepresentation by diagonalising the harmonic spring con-stant matrix K in Eq. (13), and then to construct a di-agonal friction matrix ˜ γ c in this representation. Thishas the advantage that the normal mode friction matrixcan be chosen to give critical damping (optimum canoni-cal sampling of the harmonic oscillator Hamiltonian) foreach excited (non-centroid) normal mode of the free ringpolymer. To make the connection with the notation wehave used above, ˜ γ c can then be transformed back intothe ring-polymer bead representation to give a frictionmatrix for use in Eq. (16).The relevant equations are simply C T KC = m ˜ ω , (27)˜ γ c = 2 ˜ ω , (28) γ c = C ˜ γ c C T , (29)where C is an orthogonal transformation matrix and˜ ω is a diagonal matrix of free ring polymer vibra-tional frequencies: ˜ ω kk (cid:48) = 2 ω n sin( kπ/n ) δ kk (cid:48) for k, k (cid:48) =0 , , . . . , n − ω = 0, the friction matrix γ c doesnot have any effect on the centroid ( k = 0) mode ofthe ring polymer, consistent with the constraint on γ in Eq. (17). (In Ref. , a separate Langevin thermostatwas applied to the centroid to facilitate the canonicalsampling of static equilibrium properties, but to do sohere would invalidate most of the results we have estab-lished for the autocorrelation functions obtained fromthermostatted RPMD.) Note also that, while criticaldamping of the free ring polymer internal modes is cer-tainly a good idea for the calculation of static equilibriumproperties, there is no a priori reason to suppose that itwill be the best thing to do in thermostatted RPMD. Wehave therefore explored a range of dampings by definingthe family of friction matrices γ = λ γ c , (30)where λ > λ = 1 gives critical damping, λ < λ = 0 recovers standard RPMD.The upshot is that λ = 1 / In retrospect, λ = 1 / λ = 1 for the calculation of staticequilibrium properties. We shall use this ‘optimal’ valueof λ throughout the following section unless we specifyotherwise. III. A SOLUTION TO THE RESONANCE ANDCURVATURE PROBLEMS
In the previous section we have shown that adding astochastic thermostat to the internal degrees of freedom σ ( ω ) wavenumber (cm -1 ) CMDTRPMD RPMDTRPMD(a)(b) (c)(d)
FIG. 1. Absorption cross sections for the harmonic OHmolecule: (a) and (b) CMD and TRPMD methods at 100,200, and 300 K; (c) and (d) RPMD and TRPMD methods at109, 350, and 436 K. The dotted grey line indicates the correctharmonic vibrational frequency of the model, ω OH = 3715 . − . of the ring polymer without changing its mass matrix pre-serves all of the established properties of RPMD. Here wewill demonstrate that, for several model and more real-istic systems, the use of a PILE thermostat with optimaldamping produces a method that is immune to both theresonance problem of RPMD and the curvature problemof CMD. A. Model molecules: harmonic OH and CH We shall start by considering two of the modelmolecules that were originally used to highlight the cur-vature and resonance problems. These are simplisticapproximations to the OH and CH molecules with har-monic interaction potentials V = (cid:88) bonds k b r − r e ) + (cid:88) angles k a θ − θ e ) , (31)in which the values of k b , k a , r e , and θ e are given inRef. .For each of these model molecules, we have calculatedthe dipole absorption cross section σ ( ω ) ∝ βω ˜ I ( ω ),where ˜ I ( ω ) is the Fourier transform of the Kubo-transformed autocorrelation function of the moleculardipole moment. RPMD, CMD, and TRPMD all yieldapproximations to this autocorrelation function. The re-sulting spectra are shown in Fig. 1 for the harmonic OHmolecule and in Fig. 2 for CH . For the sake of compar-ison, these simulations were performed with translationsand rotations removed using the procedure described in σ ( ω ) wavenumber (cm -1 ) CMDTRPMD RPMDTRPMD(a)(b) (c)(d)
FIG. 2. Absorption cross sections for the harmonic CH molecule: (a) and (b) CMD and TRPMD methods at 100,200, and 300 K; (c) and (d) RPMD and TRPMD methods at136, 273, and 400 K. Ref. , and at the same temperatures and with the samenumber of beads as in that study.In both figures we show in panels (a) and (c) the CMDand RPMD spectra, which are fully consistent with theresults reported in Ref. . The curvature problem ofCMD is apparent for both molecules. The high-frequencystretching peaks shift (quite dramatically) to lower fre-quencies and broaden significantly with decreasing tem-perature. The shift is more pronounced for the OHmolecule than for the CH molecule, because the pres-ence of a harmonic bending term in the potential formethane reduces the spread of the ring polymer in thedirections orthogonal to the stretch, thereby limiting thesoftening of the centroid potential of mean force alongthe stretching coordinate. For RPMD, one clearly seesthat the peaks in the spectrum are split by resonancesat temperatures where the free ring polymer frequenciescoincide with the physical frequencies (109 K and 436 Kfor OH and 136 K and 273 K for CH ). This resonanceproblem is just as severe for both molecules, because it iscaused by the coupling between the physical vibrations ofthe stretching coordinate and comparatively low ampli-tude vibrations of the internal modes of the ring polymerin the directions orthogonal to the stretch.In panels (b) and (d) we show the optimally-damped( λ = 1 /
2) TRPMD results for both molecules at the sametemperatures used for the CMD and RPMD calculations.In all cases the resonant peak splitting disappears and thecurvature problem is absent. This shows that the curva-ture problem originates from the choice of the Parrinello-Rahman mass matrix used in CMD, rather than from thethermostatting of the internal modes. As we will discussmore in detail in the next subsection, the peaks in theTRPMD spectrum are found to broaden as the temper-
RPMDTRPMD3200 3400 3600 3800 4000 wavenumber (cm -1 ) σ ( ω )
600 800 1000 (a) λ = 0.05 λ = 0.15 λ = 0.5 λ = 1.5 λ = 5(b)(c)(d)(e) x10x10x10x10x10 FIG. 3. Absorption cross sections for the harmonic OHmolecule at 109 K. The various plots compare the RPMDspectrum to the TRPMD spectrum with different dampingparameters λ . On the left side of the figure, the cross sectionshave been multiplied by a factor 10 for better visualisation. ature is lowered. This is because the internal vibrationsof the ring polymer are damped by the Langevin frictionbut not eliminated. With optimal damping, they do notappear as sharp resonances, but they do lead to an in-direct damping of the physical vibrations that becomesmore pronounced as the temperature is lowered and thespacing between the internal frequencies of the ring poly-mer decreases.To demonstrate that optimal damping is indeed thebest choice for the TRPMD method, we have investigatedthe sensitivity of the results to the choice of the dampingparameter. Figure 3 shows the TRPMD spectrum of themodel OH molecule at 109 K calculated with differentvalues of λ , in both the underdamped and overdampedregimes. The underdamped spectrum with λ = 0 .
15 givesa peak very close to the correct OH stretching frequency,but it still shows a spurious peak at about 900 cm − that corresponds to one of the internal excitations of thering polymer in a direction orthogonal to the stretch.The strongly overdamped TRPMD spectrum with λ = 5shows no trace of this spurious peak, but the peak in theOH stretching region is now significantly red shifted (al-though not nearly by so much as it is at this temperaturein CMD). Overall, optimal damping appears to be a rea-sonable compromise that washes out the spurious peakwithout dramatically red-shifting the physical vibration.Although the choice of λ does affect the precise positionof the stretching peak, there is a fairly large range ofvalues approaching λ = 1 / σ ( ω ) wavenumber (cm -1 ) CMDTRPMD RPMDTRPMD(a)(b) (c)(d)
FIG. 4. Absorption cross sections for the anharmonic OHmolecule: (a) and (b) CMD and TRPMD methods at 100,200, and 300 K ; (c) and (d) RPMD and TRPMD methodsat 109, 350, and 436 K. The dotted grey line indicates thecorrect fundamental transition frequency ω ← = 3568 cm − . position and masking the resonance artefact. B. Anharmonic OH
As a second test case, we have considered an OHmolecule with an anharmonic (Morse) interaction poten-tial V = D e (cid:104) − e − α ( r − r e ) (cid:105) , (32)where D e = hc ω e / ω e x e and α = (cid:112) µ OH hc ω e x e / (cid:126) with ω e = 3737 .
76 cm − , ω e x e = 84 .
881 cm − , and r e = 0 . Since the vibrational energy levels ofthis potential are known exactly, this test case providesa convenient way to separate the curvature problem ofCMD from the correct anharmonic red shift of the vi-brational fundamental band, and also to study how theTRPMD method performs for a more realistic potential.Figure 4 is analogous to Fig. 1, but now for the anhar-monic OH molecule. For the CMD and RPMD methods(panels (a) and (c)) we see very similar curvature and res-onance problems to those observed in the harmonic case.The position of the CMD stretching peak is noticeablyred shifted with respect to the correct anharmonic resulteven at 300 K, and the RPMD stretching peak is plaguedby resonances at all three temperatures considered in thefigure.The optimally-damped TRPMD results shown in pan-els (b) and (d) of Fig. 4 are again unaffected by the cur-vature and resonance problems. The TRPMD stretchingpeak is within 35 cm − of the correct fundamental tran-sition frequency, whereas the harmonic vibrational fre- F W H M ( c m - ) TRPMDCMDClassical temperature (K) w a v e nu m b e r ( c m - ) FIG. 5. Full width at half maximum FWHM (top) and peakposition (bottom) as a function of temperature for the CMDand TRPMD spectra shown in Fig. 4, plus additional datafrom a classical simulation and points at 1000 K for all threemethods. The dotted line in the bottom panel indicates thecorrect position of the peak at 3568 cm − . quency of the molecule (which would be obtained froma classical simulation at low temperature) is blue shiftedfrom the fundamental frequency by 170 cm − . The po-sition of the TRPMD peak is also independent of thesimulation temperature, although it does broaden as thetemperature is lowered, and more so than in the case ofthe harmonic OH molecule considered in Fig. 1.In order to quantify the extent of this broadening, wehave fitted Lorentzian curves to the CMD and TRPMDpeaks in Fig. 4. The resulting peak positions and fullwidths at half maximum (FWHM) are shown as a func-tion of temperature in Fig. 5, along with the results ofclassical simulations and additional calculations at 1000K. The peaks proved to be almost perfectly Lorentzian inall cases except at 1000 K, where they were slightly asym-metric. From the top panel of Fig. 5 it is clear that, whilethe TRPMD vibrational peak broadens with decreasingtemperature, the broadening is less severe than in thecase of CMD. The peak position, shown in the bottompanel, is also stable in the TRPMD and classical simula-tions, but not in CMD. The classical simulation producesa peak with a FWHM of only 11cm − at 100 K, but thewidth increases with temperature, and becomes compa-rable to the widths of the CMD and TRPMD peaks at1000 K. ExperimentClassical σ ( ω ) RPMDCMD
800 1200 1600 2000 wavenumber (cm -1 ) TRPMD (a) σ ( ω ) wavenumber (cm -1 ) (b) FIG. 6. (a) Top panel: experimental infrared multi-photon dissociation (IRMPD) spectrum of H O +2 in the 800-2000 cm − region at 100 K from Ref. . The other panels show simulated spectra at this temperature obtained from the Fourier transformof the dipole autocorrelation function using the CCSD(T) parameterised potential energy and dipole moment surfaces of Ref. .From top to bottom: classical molecular dynamics, RPMD, CMD, and optimally-damped TRPMD. (b) Top panel: experimentalIRMPD spectrum of H O +2 in the OH stretch region at 300 K from Ref. . The other panels are the same as in (a), but at 300K. C. The Zundel cation
Next, to illustrate the limitations of using any methodlike RPMD, CMD, or TRPMD to simulate gas phasevibrational spectra, we have considered a significantlymore complex test case: the Zundel cation H O +2 de-scribed by the CCSD(T) potential energy and dipolemoment surfaces of Huang, Braams and Bowman. The infrared spectroscopy of this cation has been stud-ied extensively in the past both theoretically andexperimentally, and it contains features that oneshould not even hope to capture using a method thatneglects real time quantum interference effects. Amongother things, the low temperature spectrum in the sharedproton region contains a doublet feature that has beeninterpreted as a Fermi resonance involving a fourth-ordercoupling between the proton transfer mode, the O-Ostretching mode, and the H-O-H wagging mode.
Thisfeature can only be captured by high levels of theorythat include anharmonic effects and a fully quantum me-chanical treatment of the nuclear dynamics, such as themulti-configurational time-dependent Hartree (MCTDH)method at 0 K. With this caveat in mind, let us first considerthe frequency range corresponding to the OH stretch(Fig. 6(b)). The present results were obtained from 2006 ps trajectories at 300 K, using 16 beads and withoutremoving rotations. Experiments at this temperature show two distinct peaks at about 3600 and 3700 cm − .As in the case of the anharmonic OH discussed above, the classical spectrum is blue shifted from these peaks by al-most 150 cm − . The RPMD spectrum is closer to theexperiment, but it has an incorrect profile and a broadtail towards lower frequencies. In a multi-dimensionalproblem such as this, there can be resonances with theinternal modes of the ring polymer associated with a vari-ety of different molecular vibrational modes, which makesit hard to assess whether and how the RPMD spectrumis contaminated by resonances. The CMD spectrum isred shifted with respect to the experiment, and consistsof a single broad peak. Given the high accuracy of thepotential energy surface, we believe that the red shiftcan be attributed to the curvature problem, and we havechecked this by decreasing the temperature: the red shiftthen becomes even more pronounced in the same way asin Figs. 1 and 4. Finally, the optimally-damped TRPMDspectrum shows neither resonances nor any curvature-related red shift. The position of the absorption band isin reasonable agreement with the experiment, althoughthe significant broadening that was also observed for theanharmonic OH molecule washes out the doublet struc-ture and yields a single peak centred at about 3650 cm − .We have also run simulations at 100 K using 64 beadsto compare RPMD, CMD and TRPMD with an experi-ment performed in the intermediate frequency region ofthe spectrum (Fig. 6(a)). All three methods generatevery broad peaks in this region, and the RPMD spectrumdisplays oscillations in its bending peak at around 1750cm − that can probably be attributed to a resonance ef-fect. Interestingly, the CMD bending peak is red-shifted0relative to experiment, suggesting that at such low tem-peratures the curvature problem can also affect lower fre-quency modes than OH stretching vibrations. The bend-ing peak from the optimally-damped TRPMD calcula-tion is closest to the experiment and shows no signs ofresonances. As discussed above, none of the methods in-vestigated here can reproduce the experimental doubletat around 1000 cm − , which is due to a subtle Fermi res-onance effect. The TRPMD calculation simply pro-duces a single broad absorption band that spans therange of all three experimental peaks between 800 and1400 cm − , and the RPMD and CMD calculations donot fare any better. D. Liquid water
Having highlighted a situation in which none of thepresent methods is really satisfactory, it is only fair toend with an application for which methods like RPMDand CMD are designed: the simulation of liquid water.This is a fundamentally different problem because the ex-act quantum mechanical dipole autocorrelation functionof a liquid decays to zero at long times, resulting in a con-tinuous spectrum. The classical approximation to thisspectrum (to which both RPMD and CMD are pinnedat high temperatures) becomes exact in the limit where β (cid:126) ω (cid:28)
1, which is not the case for the discrete spectrumof a small gas phase molecule. Methods like RPMD andCMD (and indeed even classical MD) are therefore on farfirmer ground for the simulation of liquids than they arefor isolated molecules.With this in mind, we have used classical MD, RPMD,CMD and optimally-damped TRPMD to calculate thedipole absorption spectrum of the q-TIP4P/F watermodel. This water model was chosen because it re-produces a wide variety of structural, thermodynamicand dynamical properties of liquid water in path in-tegral simulations, including the vibrational band fre-quencies which were optimised in the parameterisationof the model using a partially-adiabatic approximationto CMD. Because it is a simple point charge model ratherthan a polarisable model, it is missing the induced dipolecontribution to the change in dipole moment along theOH stretch, which is known to make a significant contri-bution to the intensity of the OH stretching band. How-ever, the frequency of the stretching band in the model isroughly correct, and the same is true of the intramolec-ular bending and intermolecular librational bands. In this case, since the volume V of the system is welldefined, it is possible to calculate the absolute signal n ( ω ) α ( ω ) = πβω cV (cid:15) ˜ I ( ω ) , (33)where n ( ω ) is the frequency-dependent refractive indexof the liquid, α ( ω ) is its Beer-Lambert absorption co-efficient, and ˜ I ( ω ) is the Fourier transform of its Kubo-transformed dipole autocorrelation function. The present calculations were performed with the i-PI path integralcode, using LAMMPS as the backend to calculate en-ergies and forces. We used 216 water molecules in a pe-riodically replicated simulation box with the experimen-tal density at 300 K. The classical results were averagedover 8 independent 100 ps trajectories initiated from anequilibrated simulation, with a time step of 0.25 fs. Theoptimally-damped TRPMD calculations were performedin exactly the same way, but with n = 32 ring polymerbeads. In the CMD calculations, the masses of the inter-nal modes of the free ring polymer were adjusted to bringthem to a frequency of 16000 cm − , and we decreasedthe integration time step to 0.025 fs to correctly inte-grate the rapid oscillations of these modes. RPMD cal-culations are typically done by averaging over hundredsof short trajectories, with the momenta resampled fromthe Boltzmann distribution between each one in order toovercome the non-ergodicity of the microcanonical ringpolymer dynamics. However, since RPMD is the λ → λ = 0 . in Fig. 7. All four simulations give very similar spectrain the low frequency librational region, where they arein good agreement with the experiment. In the interme-diate frequency bending region, the peaks of the threepath integral methods are again seen to be in agreement,but they are red shifted from the classical bending peakby around 35 cm − . Given the good agreement betweenthe three path integral methods, the absence of any evi-dent curvature problem in the CMD bending peak, andthe absence of any apparent resonance problem in theRPMD bending peak, we conclude that this is the cor-rect anharmonic red shift of the bending band for theq-TIP4P/F water model.The differences between the path integral methods oc-cur in the high frequency OH stretching region. Here theRPMD calculation shows clear indications of the reso-nance problem, although the resonances are not as pro-nounced as those seen previously using the TTM3-F water model, which has a non-linear dipole momentoperator that accentuates the resonance coupling. Theresonances are washed out in the TRPMD spectrum,which exhibits a single peak that is red shifted from theclassical peak by some 75 cm − . The CMD peak is redshifted by a further 80 cm − . This peak is clearly inthe best agreement with the experiment, but this is sim-ply because the q-TIP4P/F potential was parameterisedto agree with this experiment in a (partially adiabatic)CMD calculation. The issue of the curvature problem in CMD simula-tions of liquid water has been the subject of some de-1
Experiment
Classical n ( ω ) α ( ω ) ( c m - ) RPMD
CMD wavenumber (cm -1 ) TRPMD (a) n ( ω ) α ( ω ) ( c m - ) wavenumber (cm -1 ) (b) FIG. 7. (a) Top panel: experimental infrared absorption spectrum of room temperature liquid water from Ref. between 0 and2000 cm − . The other panels show simulated spectra at this temperature obtained from the Fourier transform of the dipoleautocorrelation function using the q-TIP4P/F water model. From top to bottom: classical molecular dynamics, RPMD,CMD, and optimally-damped TRPMD. (b) Same as in (a), but in the region between 3000 and 4000 cm − . bate. Paesani and Voth have argued, on the basis ofpartially adiabatic CMD simulations with a weak adia-batic decoupling parameter of a solvated HDO moleculein H O and D O, and comparisons of the OD and OHstretching peaks of the solute with those given by mixedquantum-classical calculations, that the curvature prob-lem is unlikely to play a significant role in simulations ofliquids like q-TIP4P/F with anharmonic stretching po-tentials. On the other hand, Ivanov et al. have ar-gued that the curvature problem is indeed present in thestretching bands of CMD simulations of neat liquid HDO,H O, and D O with both harmonic and anharmonic OHstretching potentials.In the present simulations, a comparison with the posi-tion of the TRPMD OH stretching peak, which we wouldnot expect on the basis of the other vibrational problemswe have considered above to be shifted very far fromthe correct position, suggests that the CMD peak doesshow some evidence of the curvature problem. This isnot to say that the TPRMD OH stretching frequency isexactly right, and indeed it may well not be: recall the ∼
35 cm − blue shift that this method gives from thecorrect fundamental transition frequency of the anhar-monic OH molecule (Fig. 4). In view of this uncertainty,which prevents us from drawing any more definite con-clusions here, we believe it would be worthwhile in a fu-ture study to compare the present CMD and TRPMD re-sults for q-TIP4P/F liquid water with those of a more ac-curate (and expensive) quantum mechanical calculation,for example using the recently-developed local monomerapproximation. Finally, Table I compares the diffusion coefficients and various rotational correlation times τ ηl obtained from thepresent classical, RPMD, CMD, and optimally-dampedTRPMD simulations. One sees that all three path inte-gral methods are in qualitative agreement in predictinga 12% increase in the diffusion coefficient and a ∼ . Indeed the differences between the path integralmethods for these zero-frequency dynamical propertiesare mostly within the statistical error bars of our simu-lations (obtained from 8 independent 100 ps trajectoriesfor each method). TABLE I. Dynamical properties of liquid water at 300 K ob-tained from classical, RPMD, CMD and TRPMD simulationsof the q-TIP4P/F water model. D is the diffusion coefficientand τ ηl is the l th order orientational relaxation time for molec-ular axis η . Classical RPMD CMD TRPMD D H O (˚A ps − ) 0.194(5) 0.218(3) 0.219(2) 0.217(2) τ HH1 (ps) 6.1(2) 5.0(1) 5.1(1) 5.1(1) τ OH1 (ps) 5.8(2) 4.8(1) 4.9(1) 4.8(1) τ µ (ps) 5.2(2) 4.4(1) 4.5(2) 4.4(1) τ HH2 (ps) 2.55(8) 2.15(7) 2.17(5) 2.08(6) τ OH2 (ps) 2.23(7) 1.85(5) 1.90(4) 1.81(4) τ µ (ps) 1.87(8) 1.42(4) 1.49(4) 1.45(5) IV. CONCLUDING REMARKS
In this paper, we have shown that it is possible to applyan internal mode thermostat to RPMD without alteringany of its established properties (see Sec. II), and thatthe resulting method – which is halfway between CMDand RPMD – is significantly better than either of thesealternatives for vibrational spectroscopy (see Sec. III).The results are still not perfect, and we would certainlynot recommend using any of the methods we have consid-ered here to simulate the spectroscopy of small gas phasemolecules. However, for more complex systems, rang-ing from large biomolecules to liquids, optimally-dampedTRPMD might at least be expected to provide a reason-able prediction of the anharmonic shifts in vibrationalfundamental bands, at a modest computational cost.Since TRPMD is no less valid (or ad hoc !) than stan-dard RPMD, and we have shown here that it worksbetter for vibrational spectroscopy, it is natural to askwhether it might also offer an improvement for other ap-plications, such as the calculation of diffusion coefficientsand chemical reaction rates. Does the TRPMD veloc-ity autocorrelation function of a quantum liquid suchas para-hydrogen satisfy imaginary time moment con-straints more accurately than the RPMD velocity auto- correlation function? Their transition state theory limitsare the same, but are the TRPMD transmission coeffi-cients of some simple chemical reactions (for which theexact rates can be computed for comparison) more ac-curate than the RPMD transmission coefficients? Theseare interesting questions for future work.Finally, we should stress that the PILE thermostat wehave considered here is by no means unique. By movingto a coloured noise (generalised Langevin equation) ther-mostat, in the form of a Langevin thermostat in an ex-tended momentum space, it may well be possible to findsomething better. There is clearly a tremendous amountof leeway in the construction of extended phase spaceapproximations to quantum dynamics, and it is almostcertain that we have yet to find the best one. Methodslike CMD and RPMD have simply scratched the surfaceof what might be possible.
ACKNOWLEDGMENTS
We are grateful to David Wilkins for providing us witha code for calculating rotational correlation times and toTom Markland for extensive comments on a preliminarydraft of this manuscript. This work was supported bythe Deutsche Forschungsgemeinschaft (DFG), the Wolf-son Foundation, and the Royal Society. J. Cao and G. A. Voth, J. Chem. Phys. , 5106 (1994). S. Jang and G. A. Voth, J. Chem. Phys. , 2371 (1999). I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. ,3368 (2004). B. J. Braams and D. E. Manolopoulos, J. Chem. Phys. , 124105 (2006). R. P. Feynman and A. R. Hibbs,
Quantum Mechanics andPath Integrals (McGraw-Hill, New York, 1965). D. Chandler and P. G. Wolynes, J. Chem. Phys. , 4078(1981). G. A. Voth, Adv. Chem. Phys. , 135 (1996). S. Habershon, D. E. Manolopoulos, T. E. Markland andT. F. Miller III, Ann. Rev. Phys. Chem. , 387 (2013). M. Pavese and G. A. Voth, Chem. Phys. Lett. , 231(1996). J. Lobaugh and G. A. Voth, J. Chem. Phys. , 2400(1997). Y. Yonetani and K. Kinugawa, J. Chem. Phys. , 9651(2003). T. D. Hone and G. A. Voth, J. Chem. Phys. , 6412(2004). T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. , 184503 (2005). T. F. Miller III and D. E. Manolopoulos, J. Chem. Phys. , 154504 (2005). T. D. Hone, P. J. Rossky and G. A. Voth, J. Chem. Phys. , 154103 (2006). T. E. Markland, S. Habershon and D. E. Manolopoulos, J.Chem. Phys. , 194506 (2008). T. F. Miller III, J. Chem. Phys. , 194502 (2008). S. Habershon, T. E. Markland and D. E. Manolopoulos, J.Chem. Phys. , 024501 (2009). A. R. Menzeleev and T. F. Miller III, J. Chem. Phys. ,034106 (2010). T. E. Markland, J. A. Morrone, B. J. Berne, K. Miyazaki,E. Rabani and D. R. Reichman, Nat. Phys. , 134 (2011). G. A. Voth, J. Phys. Chem. , 8365 (1993). J. Cao and G. A. Voth, J. Chem. Phys. , 6856 (1996). J. Cao and G. A. Voth, J. Chem. Phys. , 1769 (1997). S. Jang and G. A. Voth, J. Chem. Phys. , 8747 (2000). E. Geva, Q. Shi and G. A. Voth, J. Chem. Phys. , 9209(2001). I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. ,084106 (2005). I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. ,034102 (2005). R. Collepardo-Guevara, I. R. Craig and D. E. Manolopou-los, J. Chem. Phys. , 144502 (2008). R. Collepardo-Guevara, Y. V. Suleimanov andD. E. Manolopoulos, J. Chem. Phys. , 174713(2009). Y. V. Suleimanov, R. Collepardo-Guevara andD. E. Manolopoulos, J. Chem. Phys. , 044131(2011). A. R. Menzeleev, N. Ananth and T. F. Miller III, J. Chem.Phys. , 074106 (2011). N. Boekelheide, R. Salomon-Ferrer and T. F. Miller III,Proc. Nat. Acad. Sci. USA , 16159 (2011). R. Perez de Tudela, F. J. Aoiz, Y. V. Suleimanov andD. E. Manolopoulos, J. Phys. Chem. Lett. , 493 (2012). J. W. Allen, W. H. Green, Y. Li, H. Guo andY. V. Suleimanov, J. Chem. Phys. , 221103 (2013). Y. Li, Y. V. Suleimanov and H. Guo, J. Phys. Chem. Lett. , 700 (2014). A. R. Menzeleev, F. Bell and T. F. Miller III, J. Chem.Phys. , 064103 (2014). B. J. Braams, T. F. Miller III and D. E. Manolopoulos,Chem. Phys. Lett. , 179 (2006). S. Habershon, B. J. Braams and D. E. Manolopoulos, J.Chem. Phys. , 174108 (2007). A. Perez, M. E. Tuckerman and M. H. Muser, J. Chem.Phys. , 184105 (2009). M. Shiga and A. Nakayama, Chem. Phys. Lett. , 175(2008). S. Habershon, G. S. Fanourgakis and D. E. Manolopoulos,J. Chem. Phys. , 074501 (2008). A. Witt, S. D. Ivanov, M. Shiga, H. Forbert and D. Marx,J. Chem. Phys. , 194510 (2009). S. D. Ivanov, A. Witt, M. Shiga and D. Marx, J. Chem.Phys. , 031101 (2010). J. Cao and G. A. Voth, J. Chem. Phys. , 6168 (1994). R. Kubo, J. Phys. Soc. Jpn. , 570 (1957). J. O. Richardson and S. C. Althorpe, J. Chem. Phys. ,214106 (2009). T. J. H. Hele and S. C. Althorpe, J. Chem. Phys. ,084108 (2013). S. C. Althorpe and T. J. H. Hele, J. Chem. Phys. ,084115 (2013). T. J. H. Hele and S. C. Althorpe, J. Chem. Phys. ,084116 (2013). A. Horikoshi and K. Kinugawa, J. Chem. Phys. ,174104 (2005). S. Jang, A. V. Sinitskiy and G. A. Voth, J. Chem. Phys. , 154103 (2014). R. Zwanzig,
Nonequilibrium Statistical Mechanics (OxfordUniversity Press, New York, 2001), pp. 42 and 43. M. Parrinello and A. Rahman, J. Chem. Phys. , 860(1984). J. Cao and G. A. Voth, J. Chem. Phys. , 6157 (1994). M. Ceriotti, G. Bussi and M. Parrinello, J. Chem. TheoryComput. , 1170 (2010) M. Ceriotti, M. Parrinello, T. E. Markland andD. E. Manolopoulos, J. Chem. Phys. , 124104 (2010). K. P. Huber and G. Herzberg,
Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules (Van Nostrand Reinhold, New York, 1979), p. 508. X. Huang, B. J. Braams and J. M. Bowman, J. Chem.Phys. , 044308 (2005). O. Vendrell, F. Gatti and H.-D. Meyer, J. Chem. Phys. , 184303 (2007). G. Niedner-Schatteburg, Angew. Chem. Intl. Ed. , 1008(2008). X. Huang, S. Habershon and J. M. Bowman, Chem. Phys.Lett. , 253 (2008). M. Kaledin, A. L. Kaledin, J. M. Bowman, J. Ding andK. D. Jordan, J. Phys. Chem. A , 7671 (2009). M. Baer, D. Marx and G. Mathias, Angew. Chem. Intl.Ed. , 7346 (2010). F. Agostini, R. Vuilleumier and G. Ciccotti, J. Chem.Phys. , 084302 (2011). L. I. Yeh, M. Okumura, J. D. Meyers, J. M. Price andY. T. Lee, J. Chem. Phys. , 7319 (1989). K. R. Asmis, N. L. Pivonka, G. Santambrogio,M. Br¨ummer, C. Kaposta, D. M. Neumark and L. W¨oste,Science , 1375 (2003). T. D. Fridgen, T. B. McMahon, L. MacAleese, J. Lemaireand P. Maitre, J. Phys. Chem. A , 9008 (2004). N. I. Hammer, E. G. Diken, J. R. Roscoli, M. A. Johnson,E. M. Myshakin, K. D. Jordan, A. B. McCoy, X. Huang,J. M. Bowman and S. Carter, J. Chem. Phys. , 244301(2005). T. L. Guasco, M. A. Johnson and A. B. McCoy, J. Phys.Chem. A , 5847 (2011). M. Ceriotti, J. More and D. E. Manolopoulos, Comp. Phys.Comm. , 1019 (2014) S. Plimpton, J. Comp. Phys. , 1 (1995) R. W. Hall and B. J. Berne, J. Chem. Phys. , 3641(1984). J. E. Bertie and Z. Lan, Appl. Spectrosc. , 1047 (1996). G. S. Fanourgakis and S. S. Xantheas, J. Chem. Phys. ,074506 (2008). F. Paesani and G. A. Voth, J. Chem. Phys. , 014105(2010). Y. Wang and J. M. Bowman, J. Chem. Phys. , 154510(2011). H. Liu, Y. Wang and J. M. Bowman, J. Phys. Chem. Lett.3