Nuclear Quantum Effects in liquid water at near classical computational cost using the adaptive Quantum Thermal Bath
Nastasia Mauger, Thomas Plé, Louis Lagardère, Sara Bonella, Etienne Mangaud, Jean-Philip Piquemal, Simon Huppert
NNuclear Quantum Effects in liquid water at near classical computational cost usingthe adaptive Quantum Thermal Bath
Nastasia Mauger , Thomas Pl´e , Louis Lagard`ere*, Sara Bonella, ´Etienne Mangaud, Jean-Philip Piquemal*, and Simon Huppert* Sorbonne Universit´e, LCT, UMR 7616 CNRS, F-75005, Paris, France CNRS, Sorbonne Universit´e, Institut des NanoSciences de Paris,UMR 7588, 4 Place Jussieu, F-75005 Paris, France CECAM Centre Europ´een de Calcul Atomique et Mol´eculaire,´Ecole Polytechnique F´ed´erale de Lausanne, Batochimie, Avenue Forel 2, 1015 Lausanne, Switzerland
These authors contributed equally to this work CNRS, Sorbonne Universit´e, Institut des NanoSciences de Paris,UMR 7588, 4 Place Jussieu, F-75005 Paris, France (Dated: February 8, 2021)We demonstrate the accuracy and efficiency of a recently introduced approach to account fornuclear quantum effects (NQE) in molecular simulations: the adaptive Quantum Thermal Bath(adQTB). In this method, zero point energy is introduced through a generalized Langevin ther-mostat designed to precisely enforce the quantum fluctuation-dissipation theorem. We propose arefined adQTB algorithm with improved accuracy and we report adQTB simulations of liquid wa-ter. Through extensive comparison with reference path integral calculations, we demonstrate thatit provides excellent accuracy for a broad range of structural and thermodynamic observables aswell as infrared vibrational spectra. The adQTB has a computational cost comparable to classicalmolecular dynamics, enabling simulations of up to millions of degrees of freedom.
Nuclear quantum effects play a major role in a widerange of physical and chemical processes where lightatoms, and especially hydrogen, are involved[1–4]. Inparticular, a few studies point to their importance in bi-ological systems[5–7], where hydrogen-bonding is ubiqui-tous, but realistic atomic-scale simulations in that arearemain scarce. For such large and complex systems,the most common approach has been to include NQEsimplicitly, by fitting analytical potential energy surfacemodels in order to recover experimental thermodynamicproperties when performing simulations with classicalnuclei[8, 9]. This strategy potentially limits transferabil-ity and its ability to make predictions outside the fittingdata set is questionable. Furthermore, the recent devel-opments of new generation polarizable force fields[8, 10–13] and machine learning (ML) potentials[8, 14–17] haveopened perspectives for atomistic simulations of con-densed matter systems. These approaches enable highfidelity modeling of the Born-Oppenheimer (BO) energy,and reproduce advanced quantum chemical calculationsat a fraction of their computational cost. When reach-ing such precision on the BO energy, it becomes crucialto account for NQEs explicitly to accurately reproduceexperimental observation and take full advantage of thehigh accuracy achieved[18–21].The conceptual and computational complexity of themethods that account for NQEs explicitly has hinderedtheir spread to a broad community. Reliable results canbe obtained in the imaginary-time path integrals (PI)framework [22, 23], by simulating multiple classical repli-cas of the system (also called beads). PI provides a nu-merically exact reference for static properties (approx- imations have also been derived for dynamical observ-ables, as discussed below), but their numerical cost in-creases linearly with the number of replicas and can be-come very large compared to classical molecular dynam-ics (MD). Several solutions have been proposed to miti-gate this cost, such as multiple timestepping in real andimaginary time[24–27]. However, this method is based ona decomposition of the energy as a sum of cheap high-frequency and expensive low-frequency term, that is notalways feasible (in particular in ML approaches). Otherdevelopments, such as high-order PI [28, 29] or PI pertur-bation theory[30, 31] allow decreasing the number of nec-essary replicas, but the computational overhead remainsimportant - typically increasing the simulation load byan order of magnitude for hydrogen-bonded systems atroom temperature.Recently, a different approach was introduced for theexplicit treatment of NQEs with the Quantum Ther-mal Bath (QTB) [32, 33] and the related quantumthermostat[34, 35], relying on generalized Langevin ther-mostats to approximate the zero-point motion of thenuclei. Although elegant and inexpensive, these meth-ods suffer from zero-point energy (ZPE) leakage fromhigh to low frequency modes which can lead to mas-sive errors[36, 37]. One possible workaround is to com-bine the generalized thermostat approach with pathintegrals[38, 39]. Even though the number of requiredreplicas is reduced compared to standard PIMD simula-tions, the computational cost remains significant (at least6 beads are needed for water at ambient conditions[40]).In this letter, we focus on an alternative approach, theadaptive QTB (adQTB) that completely avoids resorting a r X i v : . [ phy s i c s . c h e m - ph ] F e b to PI.In adQTB, the ZPE leakage is compensated di-rectly, using a quantitative criterion derived from thefluctuation-dissipation theorem (FDT). The method wassuccessfully tested on model systems[41], but its applica-bility to more realistic problems remained to be demon-strated. In the following, we report the main theoreti-cal aspects of the QTB and adQTB methodologies andintroduce two refinements to the adQTB algorithm, im-proving its efficiency and accuracy and broadening therange of its possible applications, in particular enablingreliable constant pressure simulations. We then applythe method to liquid water. Careful comparison with PIreferences for structural and thermodynamic propertiesas well as infrared absorption spectra (IRS) shows that,contrary to standard QTB which is plagued by massiveZPE leakage, adQTB is able to capture NQEs with aremarkable accuracy, while its computational overheadremains limited to less than 25% compared to classicalMD, allowing to scale up the system size to over a millionatoms.In (ad)QTB simulations, each nuclear degree of free-dom follows a Langevin equation [32]: m i d r i dt ( t ) = − ∂V∂r i − m i γ dr i dt ( t ) + R i ( t ) i = 1 , ..., N (1)where V ( r , ..., r i , ..., r N ) is the interatomic potential ( i denotes both the atom number and the direction x , y or z ). Eq. (1) comprises a dissipative force (with friction co-efficient γ ) balanced by a random force R i ( t ) that injectsenergy in the system. In classical Langevin dynamics, R i ( t ) is a white noise, whose amplitude is proportionalto temperature. In QTB, the random force is coloredwith the following correlation spectrum: C R i R j ( ω ) = 2 m i γ i ( ω ) θ ( ω, T ) δ ji (2)where γ i ( ω ) is the random force amplitude and θ ( ω, T ) = (cid:126) ω (cid:20)
12 + 1 e (cid:126) ω/k B T − (cid:21) (3)corresponds to the average thermal energy in a quan-tum harmonic oscillator at frequency ω and tempera-ture T . Therefore, the aim of the QTB is to accountfor ZPE contributions in an otherwise classical dynam-ics by thermalizing each vibrational mode with an effec-tive energy θ ( ω, T ) instead of the classical thermal en-ergy k B T . However, in the original formulation of theQTB (where γ i ( ω ) = γ, ∀ ω ), the ZPE provided to high-frequency modes leaks towards low frequencies, whichleads to an incorrect energy distribution and can dra-matically alter the results. In adQTB, this leakage isquantified precisely using a general result of linear re-sponse theory: the quantum fluctuation-dissipation the-orem (FDT)[42]. For each degree of freedom i , we define the deviation from the FDT as:∆ F DT,i ( ω ) = Re [ C v i R i ( ω )] − m i γ i ( ω ) C v i v i ( ω ) (4) v i = dr i dt denotes the velocity, while C v i v i ( ω ) and C v i R i ( ω ) are respectively its autocorrelation and itscross-correlation spectrum with the random force R i .The FDT characterizes the frequency-dependent distri-bution of energy in a quantum system at thermal equi-librium. It implies that ∆ F DT,i ( ω ) should be zero forany ω , a condition violated in standard QTB, due toZPE leakage. In adQTB, ∆ F DT,i ( ω ) is estimated at reg-ular intervals and the coefficients γ i ( ω ) are adjusted onthe fly via a first-order dynamics to correct for this vio-lation: a negative ∆ F DT,i ( ω ) reveals an excess of energyat frequency ω , so γ i ( ω ) is reduced, and conversely forpositive deviations. The adQTB results are producedonce the γ i ( ω ) are adapted and ∆ F DT,i ( ω ) vanishes onaverage.Here, we introduce two refinements with respect toRef. 41, both of which are presented in full detail in Sup-plementary Materials. First, to improve the adaptationefficiency, the coefficients γ i ( ω ) are adjusted accordingto the mean FDT deviation, averaged over all equiva-lent degrees of freedom (i.e. over the 3 directions andover all same-type atoms). Second, we account for thefact that, due to the spectral broadening induced by thefriction force, the QTB (and adQTB) tends to slightlyunderestimate the average potential energy and to over-estimate the kinetic energy. This error (unrelated to ZPEleakage) can be predicted and quantified for a harmonicoscillator[43, 44]. We use this harmonic reference andthe deconvolution procedure of Ref. 45 to correct for thisinaccuracy: we slightly modify θ ( ω, T ) to compensate forthe effect of γ on the potential energy, while the kineticenergy is corrected a posteriori . The kinetic energy cor-rection is significant (more than 10%) and essential toenable reliable isobaric simulations, as its neglect causeslarge errors on the pressure estimation.The role of NQEs in liquid water has been extensivelyinvestigated both experimentally and theoretically[46–50]. It also represents a major challenge for the adQTB,as massive ZPE leakage takes place from the high-frequency intramolecular vibrations (O-H stretching andH-O-H bending modes) toward the slow intermolecularmotion[36, 51]. Moreover, net NQEs on the structuralproperties of water are relatively weak due to the compe-tition between two opposite trends: the stretching ZPEstrengthens hydrogen bonding, while the bending ZPEweakens it[40, 52]. The ability of the adQTB to cap-ture this subtle balance is an important indication of itsrobustness that opens perspectives for its broader appli-cation.Interatomic interactions are modeled by the q-TIP4P/F potential [53] which was included in a localversion of the Tinker-HP massively parallel package[54], g O - O (r) r (Å) 0510150.5 1 1.5 2 2.5 g O - H (r) r (Å)ClassicalQTBadQTBPIMD 01234 1 1.5 2 2.5 3 3.5 g H - H (r) r (Å) (a) (b) (c) FIG. 1: Radial distribution functions at 300 K and constant volume corresponding to the density ρ = 0 .
995 g.cm − .where we also implemented PIMD and (ad)QTB. Simu-lations are performed with 1000 water molecules. PIMDsimulations are essentially converged with 32 beads (thenumber typically reported in the literature) and requireshort timesteps, we used a 0.2 fs timestep for all methodsand checked that increasing it to 1 fs had only a limitedeffect on the accuracy of the adQTB results. In classicalLangevin MD and PIMD simulations, static averages areindependent of the parameter γ , and we use γ = 1 ps − in both cases to limit its effect on dynamical properties.On the other hand, adQTB requires relatively large fric-tion coefficients γ , to prevent vanishing of γ i ( ω ) duringadaptation (which whould results in incorrect compen-sation of the ZPE leakage[41]). We use γ = 20 ps − forall QTB and adQTB simulation (the influence of theseparameters and the scalabilty of the algorithm for largesystems is assessed in Supplementary Material).In Figure 1, the QTB and adQTB Radial Distribu-tion Functions (RDFs) are compared with their classicaland PIMD counterparts. The most salient NQE for thisobservable is the strong broadening of the intramolecu-lar peaks caused by ZPE in the O-H and H-H RDFs.This effect is very well captured by the adQTB simula-tions, while it is slightly underestimated by the standardQTB due to ZPE leakage. Apart from this, the clas-sical and quantum RDFs are very similar, due to theaforementioned competition of NQEs. In standard QTBsimulations, the leakage of the intramolecular ZPE desta-bilizes the hydrogen bond network completely and theintermolecular peaks are excessively broadened, but theadQTB procedure efficiently suppresses the leakage andthe corresponding curves almost superimpose with thePIMD reference.This analysis is further confirmed by Table I, reportingthe average of the different the q-TIP4P/F energy terms.Intermolecular interactions (labeled Coulomb and VdW)are only slightly affected by NQEs and their classicaland PIMD values are close. In standard QTB, the to-tal intermolecular energy is overestimated by more than1 kcal.mol − due to ZPE leakage, but this is well cor-rected in adQTB, where accurate values are recovered.The adQTB is remarkably precise for intramolecular en- E k AB BS VdW Coul. r OH (˚A) θ HOH (deg)Classical 2.69 0.41 1.18 2.20 -14.00 0.96 104.8QTB 8.39 1.23 5.81 1.72 -12.38 0.98 104.6adQTB 8.60 1.17 6.37 2.11 -13.76 0.98 104.7PIMD 8.41 1.17 6.26 2.15 -13.87 0.98 104.7
TABLE I: Observables at 300 K. The kinetic energy( E k ), the Angular Bending (AB), Bond Stretching(BS), Vand der Waals (VdW) and Coulomb (Coul.)energy terms of the q-TIP4P/F potential[53] arereported in kcal.mol − per water molecule (thestandard error is inferior to 0.01 kcal.mol − ), togetherwith the average oxygen-hydrogen distance r OH andmolecular angle θ HOH .ergies (labeled AB and BS) and for the kinetic energy(that comprise large amounts of ZPE). It also capturesthe elongation of the OH distance induced by NQEs,while the molecular angle is essentially unaffected. Thedielectric constant computed from the adQTB simula-tions at 300 K is 57, in good agreement with our PIMDestimation of 58 and with the value in Ref. 53, given therelatively large statistical uncertainties.Although PIMD provides a numerically exact referencefor static quantum properties, the computation of dy-namical observables, such as infrared absorption spectra(IRS), represents a much steeper theoretical challenge,subject of intense research[55–61]. There is no refer-ence method to compute IRS exactly while accounting forNQEs in large systems, but various approximations havebeen developed[62–65]. Recently, Benson et al. com-pared different state-of-the-art approximate methods forIRS calculation in liquid water and ice[66]. They showthat the Linearized semiclassical initial value represen-tation (LSC-IVR) method[64] - where time-correlationfunctions are computed from short classical trajectoriesinitialized from an approximate sampling of the Wignerdistribution - provides the most accurate IRS withintheir broad set of approaches, while the PI-based ther- n ( ω ) α ( ω ) ( a . u ) ω (cm -1 )ClassicaladQTBTRPMD 0 0.2 0.4 0.6 0.8 1 1.2 6000 7000 8000 FIG. 2: Infrared absorption spectra at 300 K (in arbitrary units).mostated ring-polymer MD (TRPMD)[55] is presented asthe cheapest available approach yielding reliable results.QTB has formerly been used with some success as anempirical method to compute approximate IRS[33, 67].Although not formally derivable from first principles ex-cept for the harmonic oscillator case, the use of QTB andadQTB for IRS calculations can be justified qualitativelyby noting that the short-time dynamics is only little af-fected by the thermostat and thus essentially classical.Therefore, much like LSC-IVR, the QTB combines clas-sical dynamics with approximate quantum initial valuesampling. Furthermore, the deconvolution procedure ofRef. 45 efficiently eliminates the main effect of the ther-mostat: the broadening of the spectral peaks.Figure 2 compares IRS computed in adQTB to thoseobtained in classical MD and TRPMD (for which amild Langevin thermostat with γ = 1 ps − was ap-plied). Compared to TRPMD, the low-frequency ab-sorption band computed with adQTB is slightly moreintense, and the bending peak (around 1500 cm − ) isa little blue-shifted and broadened. The OH stretchingpeak at 3500 cm − is sharper in adQTB than in TRPMDand its overtone at 7000 cm − has a much larger inten-sity. These two discrepancies are in favor of the adQTBapproach since TRPMD has been shown to cause a spu-rious broadening of the spectral features and to stronglyunderestimate anharmonic resonances[66]. Overall, theadQTB IRS are very similar to the LSC-IVR results re-ported in Ref. 66. This should be further confirmed bystudies on different systems but it is extremely promisinggiven the almost classical computational cost of adQTB.The dynamical properties related to slow molecularmotions, on the other hand, cannot be quantitatively as-sessed in our present adQTB implementation, due to theneed for relatively large friction coefficients. The diffu-sion coefficient D (cid:39) . s − is underestimated byalmost a factor 3 with respect to its RPMD value[53] (asimilar decrease of D is observed in classical LangevinMD using γ = 20 ps − ). The deconvolution procedure isof no help here, since D corresponds to the zero-frequency D en s i t y ( g . c m - ) Temperature (K)ClassicalQTBadQTBPIMD 9 9.5 10 10.5 11 11.5 12 12.5 240 260 280 300 320 340 360 380(b) E n t ha l p y o f V apo r i z a t i on ( kc a l . m o l - ) Temperature (K)ClassicalQTBadQTBPIMD
FIG. 3: Density of liquid water (a) and Enthalpy ofvaporization (b) at constant pressure P = 1 atm.component of the vibration spectrum, and the deconvo-lution does not provide reliable results in that spectralregion[45]. Improved diffusion estimates might be ob-tained in future works by decreasing γ selectively at lowfrequencies using a generalized friction force, or by ap-propriately redesigning the adQTB algorithm, for exam-ple using the recently introduced fast-forward Langevinmethod[68].We now explore the use of adQTB to perform fixedpressure simulations using a Langevin piston barostat[69,70]. Pressure is a challenging quantity to evaluate inthe (ad)QTB framework: its estimator is a difference be-tween two large terms that almost cancel (a potentialand a kinetic term, of the order of 10 atm each). There-fore, even small inaccuracies on either of these contribu-tions can result in non-negligible errors (see Supplemen-tary materials). The results obtained for the density asa function of temperature at P = 1 atm are shown inFigure 3.Because of the competition between NQEs, the clas-sical and PIMD results are very similar, both show-ing a characteristic bell shape with a maximum around280 K. NQEs are only responsible for a small decrease ofthe density in the intermediate temperature range (270-330 K). The standard QTB completely fails to capturethis temperature-dependence. It decreases monotonouslyand strongly overestimates the variations of the density.Compensating the leakage in adQTB allows recoveringthe overall bell shape and a good agreement with thePIMD reference. In the intermediate temperature range(most relevant for biological systems), adQTB is very ac-curate. The curvature of the density curve is only slightlyunderestimated, leading to small errors of the order of0 .
005 g.cm − in the low-temperature and in the high-temperature limits.These results show that adQTB can be a useful and in-expensive tool for constant pressure simulations of phys-ical and chemical properties. As an illustration, wepresent on Figure 3.b the enthalpy of vaporization ∆ H vap computed from the same isobaric simulations. The clas-sical ∆ H vap is systematically overestimated compared tothe corresponding PIMD values[12, 71]. When NQEsare included with the standard QTB, ∆ H vap decreasesmarkedly, and becomes even underestimated, but this isdue to ZPE leakage and the adQTB recovers an almostperfect agreement with the PIMD reference.Finally, we discuss the computational overhead of theadQTB simulations with respect to classical LangevinMD. A first additional cost comes from the generationof the colored random forces and the adaptation of the γ i ( ω ) coefficients. It represents approximately 20% ofthe total simulation time and the scalability tests pro-vided in Supplementary Materials show that, even forsystems over one million atoms, it remains inferior to25% in our present implementation - that will be furtheraccelerated using Graphics Processing Units (GPUs)[72].The q-TIP4P/F water model is particularly inexpensive,and we expect this overhead to become negligible in com-parison to atomic force calculations with more realisticmodels. A second additional cost comes from the adap-tation procedure that requires time for the γ i ( ω ) to con-verge. This necessary time can vary from one systemto another. In our liquid water simulations, we show inSupplementary Materials that with an appropriate choiceof adaptation parameters, the γ i ( ω ) coefficients can con-verge in about 10 ps. The minimum adaptation timeis thus small compared with the several ns required toreach statistical convergence on some of the physical ob- servables, as the density and the dielectric constant.The adQTB renews the original promise of the QTBmethod to provide approximate quantum simulations atan almost classical cost, but with a much improved re-liability. It is a promising alternative to PI methodsto account for NQEs explicitely in the calculation ofstatic properties as well as vibrational spectra. Combinedwith accurate ML potentials or polarizable force fields, itshould provide a powerful tool with broad applications,in particular for the large-scale simulations required inbiophysics and biochemistry. ACKNOWLEDGEMENTS
This work was made possible thanks to funding fromthe European Research Council (ERC) under the Eu-ropean Union’s Horizon 2020 research and innovationprogram (grant agreement No 810367), project EMC2.Computations have been performed at CINES on the Oc-cigen machine on grant no A0070707671. The authorsare grateful to Fabio Finocchi and Philippe Depondt formany interesting discussions.* [email protected],* [email protected],* [email protected] [1] M. Benoit, D. Marx, and M. Parrinello, Nature , 258(1998).[2] S. Miura, M. E. Tuckerman, and M. L. Klein, J. Chem.Phys. , 5290 (1998).[3] M. Rossi, P. Gasparotto, and M. Ceriotti, Phys. Rev.Lett. , 115702 (2016).[4] L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Na-ture Physics , 1 (2020).[5] P. K. Agarwal, S. R. Billeter, P. R. Rajagopalan, S. J.Benkovic, and S. Hammes-Schiffer, Proc. Natl. Acad.Sci. USA , 2794 (2002).[6] A. P´erez, M. E. Tuckerman, H. P. Hjalmarson, and O. A.Von Lilienfeld, Journal of the American Chemical Society , 11510 (2010).[7] L. Wang, S. D. Fried, S. G. Boxer, and T. E. Markland,Proc. Natl. Acad. Sci. USA , 18454 (2014).[8] G. A. Cisneros, K. T. Wikfeldt, L. Ojam¨ae, J. Lu, Y. Xu,H. Torabifard, A. P. Bart´ok, G. Cs´anyi, V. Molinero, andF. Paesani, Chemical Reviews , 7501 (2016).[9] A. V. Onufriev and S. Izadi, WIREs ComputationalMolecular Science , e1347 (2018).[10] C. Liu, J.-P. Piquemal, and P. Ren, Journal of ChemicalTheory and Computation , 4122 (2019).[11] C. Liu, J.-P. Piquemal, and P. Ren, The Journal of Phys-ical Chemistry Letters , 419 (2020).[12] S. K. Reddy, S. C. Straight, P. Bajaj, C. Huy Pham,M. Riera, D. R. Moberg, M. A. Morales, C. Knight,A. W. G¨otz, and F. Paesani, J. Chem. Phys. , 194504(2016). [13] J. Melcr and J.-P. Piquemal, Front. Mol. Biosci. , 143(2019).[14] T. Morawietz and J. Behler, The Journal of PhysicalChemistry A , 7356 (2013).[15] L. Zhang, J. Han, H. Wang, R. Car, and W. E, Phys.Rev. Lett. , 143001 (2018).[16] J. S. Smith, O. Isayev, and A. E. Roitberg, ChemicalScience , 3192 (2017).[17] A. Singraber, J. Behler, and C. Dellago, Journal ofChemical Theory and Computation , 1827 (2019).[18] G. S. Fanourgakis and S. S. Xantheas, J. Chem. Phys. , 074506 (2008).[19] F. Paesani, S. Yoo, H. J. Bakker, and S. S. Xantheas, TheJournal of Physical Chemistry Letters , 2316 (2010).[20] L. Pereyaslavets, I. Kurnikov, G. Kamath, O. Butin,A. Illarionov, I. Leontyev, M. Olevanov, M. Levitt, R. D.Kornberg, and B. Fain, PNAS , 8878 (2018).[21] B. Cheng, E. A. Engel, J. Behler, C. Dellago, and M. Ce-riotti, Proceedings of the National Academy of Sciences , 1110 (2019).[22] R. P. Feynman, A. R. Hibbs, and D. F. Styer,Quantum mechanics and path integrals (Courier Corpo-ration, 2010).[23] D. Chandler and P. G. Wolynes, J. Chem. Phys. , 4078(1981).[24] T. E. Markland and D. E. Manolopoulos, J. Chem. Phys. , 024105 (2008).[25] X. Cheng, J. D. Herr, and R. P. Steele, J. Chem. Th.Comput. , 1627 (2016).[26] V. Kapil, J. VandeVondele, and M. Ceriotti, J. Chem.Phys. , 054111 (2016).[27] O. Marsalek and T. E. Markland, J. Chem. Phys. ,054112 (2016).[28] A. P´erez and M. E. Tuckerman, J. Chem. Phys. ,064104 (2011).[29] V. Kapil, J. Behler, and M. Ceriotti, J. Chem. Phys. , 234103 (2016).[30] I. Poltavsky and A. Tkatchenko, Chem. Sci. , 1368(2016).[31] I. Poltavsky, V. Kapil, M. Ceriotti, K. S. Kim, andA. Tkatchenko, . Chem. Th. Comput. , 1128 (2020).[32] H. Dammak, Y. Chalopin, M. Laroche, M. Hayoun, andJ.-J. Greffet, Phys. Rev. Lett. , 190601 (2009).[33] Y. Bronstein, P. Depondt, F. Finocchi, and A. M. Saitta,Phys. Rev. B , 214101 (2014).[34] M. Ceriotti, G. Bussi, and M. Parrinello, Phys. Rev.Lett. , 030603 (2009).[35] M. Ceriotti, G. Bussi, and M. Parrinello, J. Chem. The-ory Comput. , 1170 (2010).[36] J. Hern´andez-Rojas, F. Calvo, and E. G. Noya, J. Chem.Theory Comput. , 861 (2015).[37] F. Brieuc, Y. Bronstein, H. Dammak, P. Depondt,F. Finocchi, and M. Hayoun, J. Chem. Theory Com-put. , 5688 (2016).[38] F. Brieuc, H. Dammak, and M. Hayoun, J. Chem. The-ory Comput. , 1351 (2016).[39] M. Ceriotti, D. E. Manolopoulos, and M. Parrinello, J.Chem. Phys. , 084104 (2011).[40] M. Ceriotti, W. Fang, P. G. Kusalik, R. H. McKenzie,A. Michaelides, M. A. Morales, and T. E. Markland,Chemical Reviews , 7529 (2016).[41] E. Mangaud, S. Huppert, T. Pl´e, P. Depondt, S. Bonella,and F. Finocchi, J. Chem. Theory Comput. , 2863(2019). [42] R. Kubo, Rep. Prog. Phys. , 255 (1966).[43] J.-L. Barrat and D. Rodney, J. Stat. Phys. , 679(2011).[44] M. Basire, D. Borgis, and R. Vuilleumier, Phys. Chem.Chem. Phys. , 12591 (2013).[45] M. Rossi, V. Kapil, and M. Ceriotti, J. Chem. Phys. , 102301 (2018), https://doi.org/10.1063/1.4990536.[46] J. A. Morrone and R. Car, Phys. Rev. Lett. , 017801(2008).[47] F. Paesani and G. A. Voth, The Journal of Physi-cal Chemistry B , 5702 (2009), pMID: 19385690,https://doi.org/10.1021/jp810590c.[48] G. Reiter, J. C. Li, J. Mayers, T. Abdul-Redah, andP. Platzman, Brazilian Journal of Physics , 142 (2004).[49] G. Hura, J. M. Sorenson, R. M. Glaeser, andT. Head-Gordon, J. Chem. Phys. , 9140 (2000),https://doi.org/10.1063/1.1319614.[50] G. Romanelli, M. Ceriotti, D. E. Manolopoulos, C. Pan-talei, R. Senesi, and C. Andreani, The Journal of Phys-ical Chemistry Letters , 3251 (2013).[51] S. Habershon and D. E. Manolopoulos, J. Chem. Phys. , 244518 (2009).[52] X.-Z. Li, B. Walker, and A. Michaelides, Proc. Natl.Acad. Sci. USA , 6369 (2011).[53] S. Habershon, T. E. Markland, and D. E. Manolopoulos,J. Chem. Phys. , 024501 (2009).[54] L. Lagard`ere, L.-H. Jolly, F. Lipparini, F. Aviat,B. Stamm, Z. F. Jing, M. Harger, H. Torabifard, G. A.Cisneros, M. J. Schnieders, N. Gresh, Y. Maday, P. Y.Ren, J. W. Ponder, and J.-P. Piquemal, Chem. Sci. ,956 (2018).[55] M. Rossi, M. Ceriotti, and D. E. Manolopoulos, J. Chem.Phys. , 234116 (2014).[56] T. J. Hele, M. J. Willatt, A. Muolo, and S. C. Althorpe,J. Chem. Phys. , 134103 (2015).[57] J. Beutier, R. Vuilleumier, S. Bonella, and G. Ciccotti,Mol. Phys. , 2894 (2015).[58] M. Ceotto, G. Di Liberto, and R. Conte, Phys. Rev.Lett. , 010401 (2017).[59] M. Basire, F. Mouhat, G. Fraux, A. Bordage, J.-L. Haze-mann, M. Louvel, R. Spezia, S. Bonella, and R. Vuilleu-mier, J. Chem. Phys. , 134102 (2017).[60] G. Trenins, M. J. Willatt, and S. C. Althorpe, J. Chem.Phys. , 054109 (2019).[61] T. Pl´e, S. Huppert, F. Finocchi, P. Depondt, andS. Bonella, J. Chem. Phys. , 114114 (2019),https://doi.org/10.1063/1.5099246.[62] J. Cao and G. A. Voth, J. Chem. Phys. , 10070 (1993).[63] J. Cao and G. A. Voth, J. Chem. Phys. , 5093 (1994).[64] W. H. Miller, The Journal of Physical Chemistry A ,2942 (2001).[65] I. R. Craig and D. E. Manolopoulos, J. Chem. Phys. ,3368 (2004).[66] R. L. Benson, G. Trenins, and S. C. Althorpe, FaradayDiscussions , 350 (2019).[67] Y. Bronstein, P. Depondt, L. E. Bove, R. Gaal, A. M.Saitta, and F. Finocchi, Phys. Rev. B , 024104 (2016).[68] M. Hijazi, D. M. Wilkins, and M. Ceri-otti, J. Chem. Phys. , 184109 (2018),https://doi.org/10.1063/1.5029833.[69] S. E. Feller, Y. Zhang, R. W. Pastor, and B. R. Brooks,J. Chem. Phys. , 4613 (1995).[70] M. Ceriotti, J. More, and D. E. Manolopoulos, ComputerPhysics Communications , 1019 (2014). [71] B. Guillot and Y. Guissani, J. Chem Phys.108