Progress and challenges in ab initio simulations of quantum nuclei in weakly bonded systems

Atomistic simulations based on the first-principles of quantum mechanics are reaching unprecedented length scales. This progress is due to the growth in computational power allied with the development of new methodologies that allow the treatment of electrons and nuclei as quantum particles. In the realm of materials science, where the quest for desirable emergent properties relies increasingly on soft weakly bonded materials, such methods have become indispensable. In this Perspective, an overview of simulation methods that are applicable for large system sizes and that can capture the quantum nature of electrons and nuclei in the adiabatic approximation is given. In addition, the remaining challenges are discussed, especially regarding the inclusion of nuclear quantum effects (NQEs) beyond a harmonic or perturbative treatment, the impact of NQEs on electronic properties of weakly bonded systems, and how different first-principles potential energy surfaces can change the impact of NQEs on the atomic structure and dynamics of weakly bonded systems


I. INTRODUCTION
Quantum mechanical simulations with atomistic resolution are currently reaching length and time scales that had been almost unimaginable only a few decades ago.][6][7] The fundamental theoretical basis for the majority of the atomistic simulations that can be found in the literature to date is the celebrated Born-Oppenheimer (BO) approximation.Within this approximation, the quantum dynamics of electrons and nuclei can be treated separately, and one is free to employ further approximations either on the electronic or on the nuclear problem, as shown in Fig. 1(a).
The need for a quantum mechanical treatment of electrons has long been accepted. 8Despite this fact, classical-like functional forms (or force fields) that do not explicitly treat the electronic degrees of freedom individually but as an integral part of an effective atom can still offer the best compromise between cost and accuracy for simulations involving system sizes beyond a few thousand atoms or time scales reaching beyond the nanosecond regime. 9These force fields can often deliver good qualitative insights, but they typically lack broad transferability and quantum mechanical accuracy.Treating the electronic degrees of freedom with "first-principles" methods restores the transferability, while improving the accuracy, of the overall approach.
1][12][13] However, as experimental measurements become more precise (e.g., single molecule or time resolved ultrafast experiments) and theoretical methods that can treat systems on a full quantum mechanical picture become more computationally feasible, examples spanning all classes of materials start to emerge, where the quantum nature of the nuclei provokes A nuclear configuration provides the commonly called "external potential" (V ext ) that the electrons are subject to, while the ground state expectation value of the electronic Hamiltonian provides the Born-Oppenheimer potential (V BO ), which enters as the potential energy of the nuclear degrees of freedom.Under this approximation, it is possible to treat the nuclei as quantum or classical particles (b).In turn, in both cases, one can assume that the forces acting on the nuclei are approximately harmonic, or not, being described by the gradient with respect to nuclear positions q of the Born-Oppenheimer potential.
a qualitative change in several quantities, 14 even in the absence of electronically non-adiabatic transitions.
In some situations, the quantum nature of the nuclei has been acknowledged: The addition of zero-point energy (ZPE) corrections within the harmonic approximation is a standard practice, for example, when analyzing relative and reaction energies in the gas-phase and in some solid-state problems.However, treating nuclei as quantum particles still presents a challenge in a number of ways, including (i) going beyond the harmonic approximation, (ii) addressing nuclear dynamics, and (iii) addressing soft materials and liquids.For the most part, these calculations would require quantum nuclei with anharmonic forces [see Fig. 1(b)].Taking these effects into account can lead to a proper treatment of electron-phonon coupling, which has several implications for superconductivity, metal-insulator transitions, charge transfer and transport, and others.
This Perspective aims to provide an overview of the methodology that can be employed to treat both electrons and nuclei as quantum particles in large-scale atomistic simulations, with a focus on the field of ab initio path integral molecular dynamics (PIMD).It does not cover any techniques going beyond the adiabatic approximation for electronic and nuclear degrees of freedom.This Perspective includes considerations to estimate whether nuclear quantum effects (NQEs) are relevant for a specific problem and how the shape of the potential energy surface can influence the degree to which NQEs lead to the observed structures and dynamics of the system.Situations where the use of perturbative approaches to include NQE may not be sufficient are shown, and the impact of these nuclear fluctuations on the electronic structure of weakly bonded systems is discussed.

II. QUANTUM MECHANICS FOR ELECTRONS IN LARGE SCALE SIMULATIONS
Under the Born-Oppenheimer approximation, the electronic structure problem is solved at a fixed configuration of the nuclei, which defines the external potential Vext, and the electrons create the Born-Oppenheimer (BO) potential V BO that defines the potential energy surface where the nuclei move.Accurately capturing potential energy surfaces, polarization, electronic density rearrangement, etc., is essential for a reliable assessment of the impact of NQEs on diverse properties of materials.
The BO potential can be obtained with an electronic structure method of choice, for example, using wave-function based methods or density-functional theory (DFT).Such methods, which consider electrons as explicit quantum particles, are very accessible and tractable nowadays.A recent special issue of the Journal of Chemical Physics 15 highlights the advances in many electronic structure packages, showing how algorithms to solve electronic structure and dynamics problems have evolved together with computer architectures in the past few decades.Because of these algorithmic developments, very accurate methodologies, such as coupled cluster theory, have become available to treat larger systems including periodic solids. 16,17ithin the realm of DFT, one can obtain very reliable results by employing good approximations to the exchange-correlation potential, νxc = δExc[ρ]/δρ, where ρ is the electronic density.Even though all known approximations to νxc are not exact, they yield, in many cases, good electronic properties for molecules, liquids, and solids.In addition, they can provide an accurate potential energy surface, from which structural and dynamical properties of the nuclei can be calculated, [18][19][20][21][22][23] attesting to their capacity to deliver reliable results for diverse nuclear configurations.In particular, approximations that take into account fractions of exact exchange in this term mitigate the self-interaction error, 24 and range-separated functionals can, moreover, approximate electronic screening. 25This family of so-called "hybrid" functionals is more computationally demanding, but recent implementations of these methods are very efficient and allow for the treatment of either thousands of atoms 26 or thousands of force evaluations for smaller systems. 27Tailor-made and local hybrid functionals capable of delivering good accuracy for interfaces have also been recently developed. 28,291][32][33][34] These interactions are highly non-local and thus absent from standard DFT functionals due to their inherent spatial locality.7][38][39][40] With these functionals and corrections, dispersion interactions can be accurately captured in a wide range of systems spanning from condensed phase solids and surfaces to soft molecular matter.So far, no correction or functional in this area was shown to perform equally well across several different material classes, which still presents a challenge especially to new material architectures.
There are still a number of open challenges to be tackled within electronic structure theory, in general, and within DFT, in particular. 24The lack of an exact density functional means that a very large number of exchange-correlation approximations exist and it is often difficult to evaluate whether a particular one can yield predictive results for a given (large-scale) problem.Employing approximations that can capture the relevant physics of the problem at hand guarantees a certain degree of reliability.However, without the availability of benchmarks from high-level quantum-chemistry methods or appropriate experiments, it is usually not possible to ascertain the accuracy of DFT approximations for certain properties of complex systems.Nevertheless, the quality of the current DFT approximations and the possibility of applying many-body perturbation methods 41 or quantum-chemistry wave-function methods to larger systems have led to a situation where a proper theoretical treatment of the nuclear problem becomes necessary to unravel further physical phenomena that can be assessed in modern experiments.

III. QUANTUM MECHANICS FOR NUCLEI IN LARGE SCALE SIMULATIONS
Regarding the nuclei, the simplest and most straightforward way to account for their quantum nature on first-principles potential energy surfaces is the so-called harmonic approximation.A quantum-mechanical Hamiltonian for the nuclei is proposed, in which the potential is considered to be of a quadratic form with respect to the nuclear coordinates.This is achieved by a truncation after the second order of a Taylor expansion of the potential around a minimum, where qi are nuclear coordinates, N is the number of atoms in the system, and 0 denotes a global or local minimum of the BO potential.
The second term on the right of Eq. ( 1) involves the second derivative of the potential with respect to nuclear positions (the Hessian) and can be written in diagonal form by a coordinate transformation.The solution of such a Schrödinger equation for the nuclei yields the nuclear vibrational modes and their respective energies.These solutions lead to several analytical results for diverse quantities, i.e., the evaluation of ZPE, free energies, vibrational spectra, and reaction rates. 42,43In particular, for periodic solids, such a formulation presents the advantage of easily accounting for the phonon dispersion within the Brillouin zone by considering the periodicity of atomic displacements and their phases in reciprocal space. 44However, what this picture completely disregards is the coupling between different vibrational modes, as well as anharmonic terms pertaining to any given mode.In Figs.2(a) and 2(b), a pictorial representation of a harmonic approximation to a considerably anharmonic two-dimensional double-well potential is shown.In this potential, it is straightforward to visualize how the coupling between two coordinates is neglected with this approximation, and the anharmonic terms around the chosen minimum are also absent.
There are different routes to include anharmonic contributions in the nuclear motion.6][47][48] Another route is to consider higher-order terms in a perturbative Taylor series expansion of the potential [48][49][50][51][52][53] or, for certain quantities, simply include anharmonic corrections from fitted potential forms to vibrational modes that are particularly anharmonic. 54However, the cost of the more rigorous theories tends to have a very steep scaling with the number of degrees of freedom, making them feasible only for small systems or for larger systems if a small subset of essential degrees of freedom is sufficient.In general, such theories in their most widely implemented forms are only applicable for soft matter and liquids where approximations can be found.The challenges in these cases involve the existence of multiple quasi-degenerate local minima, the lack of a well-defined reference configuration, and situations where anharmonic terms cannot be treated as just a small perturbation.
In contrast, a complete anharmonic picture can be obtained by statistically sampling the full BO potential.This can be achieved through ab initio molecular dynamics, [55][56][57] but with a caveat: evolving the nuclei in time (even if only for sampling purposes) supposes they follow Newtonian equations of motion generated by a The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcpclassical Hamiltonian and thus behave as classical particles.Instead, the method of ab initio PIMD, 10,12 which was pioneered in the 1980s and has gained broad popularity in the past decade, allows one to achieve a quantum mechanical description of both nuclei and electrons within the Born-Oppenheimer approximation.This method delivers quantum statistics for the nuclear degrees of freedom by sampling the potential of a special classical system, making use of the "quantum-classical isomorphism." 58Under the assumption of distinguishable particles and the BO approximation, 59 the quantum partition function in the position representation within the (ab initio) path integral formalism can be written as where P is a convergence parameter that corresponds to the discretization of imaginary time, or the commonly called "beads" of a ring polymer, βP = 1/(PkBT), ωP = PkBT/ ̵ h, and q (i) j represents the position degree of freedom i in the jth bead of the ring polymer.In practice, these imaginary-time slices, or beads, represent replicas of the full physical system of interest, which are connected to each other through the harmonic spring terms.The reintroduction of a kinetic energy term in the exponential allows the formulation of a Hamiltonian that can be sampled through molecular dynamics in the extended phase space of the ring polymer, 10 namely, j , . . ., q (3N) j ). ( When P tends to infinity, sampling configurations with the equations of motion generated by HP gives access to the quantum Boltzmann distribution of particles in the system and hence guarantees the correct quantum statistics in the limit of distinguishable particles.The masses m ′ entering the kinetic energy term in Eq. ( 3) can assume any value, as they only serve to assist with sampling.Unlike ab initio molecular dynamics with classical nuclei, where the dynamics algorithm is both a tool for statistical sampling and a means of real time propagation of nuclear coordinates, in PIMD, the time evolution of the equations of motion serves only as a sampling tool.It is possible to formulate real-time path integral methods to access the quantum real-time evolution of nuclei.However, even with notable progress being achieved in the past years, 60 the applicability of these methods remains rather limited due to the sign problem arising from the oscillatory nature of the real-time propagator.Semiclassical methods, based on different semiclassical approximations of the real-time propagator, can overcome the sign problem. 61In particular, the ring-polymer instanton method, 62 which uses trajectories in imaginary time to describe nuclear tunneling along well-defined reaction coordinates, has been successfully applied to quite high-dimensional systems, including molecules adsorbed on metallic surfaces. 63,64ttempts to approximate quantum dynamics with techniques based on imaginary-time path integrals [65][66][67][68] can also circumvent the sign problem and thus dramatically extend their applicability realm-but these approximations to quantum dynamics are often not mathematically justified.9][70][71][72][73] Thermostatted ring polymer molecular dynamics (TRPMD) 67 is a very efficient approximation in this area, which makes it useful when working with costly ab initio potentials.This approximation boils down to the dynamics in the extended phase space of the ring polymer (with m ′(i) = m (i) ), as proposed in ring polymer molecular dynamics (RPMD), 66 but with thermostats attached only to the internal modes of the ring polymer.The approximation to Kubo-transformed correlation functions cAB in RPMD/TRPMD is given by where AP(q) = ∑ P j=1 A(q j )/P and q (p) is the vector of all 3N Cartesian positions (momenta) of the system.The time evolution of the position and momenta in TRPMD is given by where ξ k (t) is a Gaussian-distributed random number with ⟨ξ k (t)⟩ = 0 and ⟨ξ k (0)ξ k (t)⟩ = δ(t) and γ is a friction matrix determined in the following way.We define the matrix K with elements ) and its transformation to the internal modes of the free ring polymer K = C T KC = ω2 , where the transformation matrices C are given in Ref. 74 and ωkk ′ = 2ωP sin(kπ/P)δ kk ′ .The original scheme proposes γ = ω (γ = CγC T ), which ensures optimal damping of the potential.In addition, γ00 = 0 for the centroid mode, which ensures that the formalism maintains all quantummechanical limits that the RPMD autocorrelation functions were shown to obey. 67Because the formalism imposes only a few restrictions on the friction coefficients, other white-noise friction parameters inspired by harmonic-well problems have been proposed, 75 and the extension of the formalism to generalized Langevin equation (GLE) thermostats has also been explored. 76The approach is applicable particularly to vibrational spectra of condensed phase systems 19,21 and molecules with many degrees of freedom 23,77 at The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcpmoderate temperatures.It is free from the spurious artifacts of RPMD and centroid molecular dynamics on vibrational spectra, 78 albeit introducing an unphysical broadening to the line shape of the peaks. 67This broadening could be partially mitigated with GLE thermostats. 76At a higher simulation cost, path integral Liouville dynamics can deliver excellent results for the vibrational spectra of small molecules. 68he successes of these approximate methods rely on their exact treatment of quantum Boltzmann statistics (i.e., the quantum statistics of distinguishable particles) for the initial conditions and their conservation of the quantum Boltzmann distribution along different types of time evolution based on classical trajectories.Their drawbacks stem from their lack of any quantum-mechanical phase information and the fact that their largely ad hoc nature makes it difficult to construct systematic improvements.An interesting recent development on this front has been the derivation of an approximate quantum dynamics technique called Matsubara dynamics 79 in which quantum statistics and classical trajectories can be combined for real time properties through a rigorous route starting from first-principles quantum dynamics.Even though this method also suffers from the sign problem (rendering it impractical), it is possible to see how methods such as RPMD, TRPMD, and centroid molecular dynamics (CMD) 65 can be obtained from it. 75,80This theory has already opened the door to the improvement of some of these approximate methods 81 and offers potential avenues for further developments.

IV. WHEN ARE NUCLEAR QUANTUM EFFECTS IMPORTANT?
NQEs on equilibrium and response properties are more pronounced at lower temperatures and for lighter particles.A quantity that illustrates this particularly well is the thermal de Broglie wavelength, given by Even though this expression assumes non-interacting free particles, it tends to be a good estimator of the importance of NQEs.Classical Boltzmann statistics typically break down if λ is of the same order of, or larger than, the typical interparticle spacing.As an example, for a proton at room temperature, λ H = 1.00 Å, while for a heavy nucleus such as Cu, λ Cu = 0.12 Å.However, at a much lower temperature such as 10 K, these values are λ H = 5.50 Å and λ Cu = 0.69 Å.One telltale sign of the importance of NQE is the appearance of geometric isotope effects.Classically, the mass of an atom cannot change static structural properties of materials, but within quantum mechanics, it can.However, the magnitude of these changes depends on the character of the potential energy surface.If the local potential energy is approximately harmonic, the effect will be negligibly small.However, if it is not, the effect can be large.
Nevertheless, even classically, the atomic mass changes the dynamical properties of materials, such as diffusion and transport coefficients or vibrational spectra.In a harmonic potential, quantum and classical mechanics would predict the same mass scalings for vibrational frequencies, but the occupation of the vibrational states would differ.A useful quantity to analyze these properties is the ratio between the harmonic estimate of the ZPE of vibrational states and their thermal occupation, ̵ hω/(kBT), which, if much larger than 1, hints at the (possibly large) impact of NQEs.Furthermore, if the potential is strongly anharmonic, the dependence of the vibrational energies on the mass of the atoms can be very different in quantum and classical mechanics. 77eyond effects related to ZPE, tunneling is perhaps the most pronounced NQE in high-dimensional systems at finite temperatures.Especially in chemical reactions, a good estimator for the importance of tunneling is the crossover temperature Tc, below which nuclear tunneling will be pronounced.Assuming a parabolic barrier, Tc = ̵ hω † /(2πkB), 82 where ω † is the absolute value of the imaginary frequency at the transition state geometry.Deviations from a parabolic barrier, however, lead to cases where this crossover temperature ceases to be a good estimate for tunneling contributions. 83ne could continue listing a few of these "rules-of-thumb," such as the estimation of the Debye temperature for heat conduction studies, and others.However, what is clear is that these estimations are based on predictions that employ harmonic potentials, as results can be obtained analytically.Indeed, NQEs are relatively easy to account for when most degrees of freedom behave harmonically or only exhibit small deviations from the harmonic behavior.This fact underlies a number of long-standing successes in the field of solid-state physics. 51,84Challenges arise, however, when the anharmonic nature of potential energy surfaces is essential and especially when systems display complex electronic structural rearrangements that require an ab initio treatment.

V. BEYOND THE HARMONIC APPROXIMATION WITH QUANTUM ELECTRONS AND NUCLEI: EXAMPLES
The considerations discussed in Sec.IV make it clear that there is a sort of "competition" involving anharmonicity and NQE, which is depicted in Fig. 3 using two models as examples.In Fig. 3(a), the anharmonic contributions to the fluctuations of nuclear positions ⟨q 2 ⟩ calculated from classical or quantum mechanics for a model where the anharmonicity is comparable across the whole frequency range are plotted.These anharmonic contributions are similar at lower frequencies and quite different at higher frequencies in this case.The frequency range for which the classical and quantum anharmonic contributions coincide increases with the increase in temperature.Furthermore, while in the classical case, the anharmonic contributions to these fluctuations are small across all frequencies at low temperatures, in the quantum case, these contributions are large for higher frequencies at all temperatures.However, as shown in Fig. 3(b), in a model where low-frequency vibrational modes present a higher degree of anharmonicity than higher frequency ones (which is often observed in real materials, as exemplified below), there is only a finite region at intermediate frequency ranges and intermediate to lower temperatures where quantum anharmonic contributions significantly differ from the classical estimate and are, at the same time, of a significant magnitude (because at higher frequencies, they approach zero for both classical and quantum nuclei in this model).Systems where both anharmonicity and NQE play a significant role at temperatures of interest typically involve light atoms moving on a potential energy surface in which the relevant coordinates are highly anharmonic (e.g., double-well potentials).

The Journal of Chemical Physics
For static properties of materials, it is indeed often found that even when anharmonic effects are pronounced, they can be wellcaptured by classical statistical sampling of the nuclear potential energy surface, and the quantum-nuclear contributions, although important, can be well described within the harmonic approximation.In fact, for situations where it is possible to define a well-suited reference configuration, one can systematically break down the different contributions to a fully anharmonic quantum mechanical quantity into classical and quantum components, as well as harmonic and anharmonic contributions to each.An example of such a breakdown for free energies of weakly bonded molecular crystals has been discussed in Ref. 86.A conclusion from that work, which has also been observed in other materials such as metal-organic frameworks and ice, 87,88 is that the largest contributions to anharmonic motions stem from low energy vibrational modes, which can be treated classically at relevant temperature ranges.This observation resonates with what has been recently reported in a study involving several solid state materials, including inorganic perovskites. 89herefore, for example, in the estimation of free energy differences of molecular crystals, only a small error was observed by treating the anharmonic contributions to this quantity classically and the quantum contributions within the harmonic approximation. 48,86owever, there are situations where the harmonic approximation for NQEs is truly insufficient, and going beyond this approximation unravels important effects also in the electronic structure.One such situation involves weakly bonded interfaces in which light atoms such as hydrogen play a prominent role.This is exemplified in Fig. 4(a) for one water molecule adsorbed on a 2 × 2 Pd(111) surface with four layers (FHI-aims program light settings, 90 6 × 6 × 1 kpoints, 100 Å vacuum).The anharmonic contributions to the forces of this system were estimated following the procedure presented in Ref. 89.In short, the atomic displacements were sampled according to both classical and quantum Boltzmann statistics for the nuclei at T = 100 K, with a harmonic potential defined by the Hessian matrix of the minimum-energy structure, as obtained with DFT employing the Perdew-Burke-Ernzerhof (PBE) functional including the vdW interactions of Ref. 91.Subsequently, the root mean square deviation Δi = √ ∑ S s (F DFT s,i − F h s,i ) 2 /S was calculated, where F DFT/h s,i is the force of degree of freedom i in sample s calculated from the DFT/harmonic potential.A sample size of S = 100 was considered.For ease of comparison, Δi obtained with quantum and classical nuclei were normalized by the standard deviation σi of the force of the corresponding degree of freedom, as sampled from quantum Boltzmann statistics.We define εi = Δi/σi as a measure for comparison.The result of this calculation for all the Cartesian degrees of freedom of this system, as well as for the normal modes of the adsorbate, is presented in Figs.4(a) and 4(b).Anharmonic contributions are larger for the degrees of freedom of the adsorbate and the top/bottom layers of the slab.Interestingly, the z components of the top Pd(111) layer and all the coordinates of the adsorbate show a significantly higher degree of anharmonicity when the nuclei are treated as quantum particles.Analyzing the same quantities in the normal mode representation reveals that the most anharmonic modes belonging to the adsorbate also show a pronounced contribution from anharmonic NQEs, particularly for the modes that bend the hydrogen atoms toward the surface.
These anharmonic NQEs can affect the electronic structure of these interfaces-an effect that can only be captured when both the electronic and nuclear degrees of freedom are treated within quantum mechanics.A striking example of these effects can be seen in a of Chemical Physics related system of water molecules adsorbed on the steps of a Pt(221) slab.The work function change Δϕ in the interface was calculated with quantum and classical nuclei 85 by performing ab initio PIMD using an acceleration scheme that reduces the necessary amount of PIMD replicas in weakly bonded interfaces proposed in Ref. 85.As reported in that work, a significant difference in Δϕ is introduced when considering quantum nuclei, as reproduced in Fig. 4(c).The origin of this effect was explained to be the increase in magnitude of the molecular dipole component pointing toward the surface upon the inclusion of NQEs.Also for a largely non-polar molecule, namely, cyclohexane adsorbed on Rh(111), NQEs were shown to impact the Δϕ in both experiment and theory. 34,92However, in this case, the origin of the effect does not lie on anharmonic NQEs on the internal degrees of freedom of the adsorbate.As depicted in Fig. 4(d), the anharmonicity on the molecule-surface interaction coordinate, which includes an H-Rh bond, causes the adsorbate to lie at different distances from the surface if nuclei are classical or quantum.This changes the effective surface dipole through the modulation of the density overlap of the adsorbate and surface and the pushback effect.Therefore, a variety of consequences of the interplay between anharmonic NQEs and electronic structure can emerge at interfaces.
The potential energy surface itself plays a decisive role in how NQEs affect the structure and dynamics of any given system.One well-known and rather extreme example is that of the competing NQEs in liquid water. 93When analyzing the impact of NQEs in dynamical properties of water with PIMD-based simulations, Habershon, Markland, and Manolopoulos noticed that they were substantially overestimated when calculated on potential energy surfaces where individual water monomers were considered to be rigid, or where the internal degrees of freedom were described by simple harmonic terms.The evaluation of self-diffusion coefficients with fully flexible (and anharmonic) water molecules showed that the impact of NQEs on this quantity is small.However, this is not so because NQEs in water are negligible.Instead, NQEs in directions perpendicular to the hydrogen bond (H-bond) between molecules make the network weaker, contributing to an increase in the diffusion coefficient-but this is compensated by NQEs acting parallel to the H-bonds, which increase the molecular dipole and strengthen intermolecular interactions, causing a decrease in the diffusion coefficient.Ab initio potential energy surfaces do not contain such drastic approximations on any degree of freedom, but they can still considerably change the impact of NQEs on certain properties.For liquid water, excellent agreement with experiment for diffusion coefficients and vibrational properties has been reported 21 when employing TRPMD with the revPBE0 functional, and these simulations also confirmed the competing NQE picture.These competing effects have also been observed in other complex H-bonded materials, such as biomolecules. 94,95urther examples where changing the ab initio potential can provoke substantial changes in the impact of NQE on dynamical properties of different systems are shown in Fig. 5.In Ref. 71 the diffusion of protons and hydroxide ions was studied with ab initio MD and TRPMD.One advantage of using ab initio potentials is that it was possible to define an unbiased descriptor for the position of the proton at any given time in the simulation based on the position of the center of the Boys orbitals and an associated machine-learning clustering algorithm. 71A depiction of the location of the proton in one simulation snapshot is shown on the top of Fig. 5 this technique, it was observed that the diffusion coefficient of a proton in a water wire at 300 K is almost doubled when considering NQEs on a potential energy surface that does not contain long-range van der Waals interactions [DFT with the Becke-Lee-Yang-Parr (BLYP) exchange correlation functional].However, the diffusion remains basically unchanged when NQEs are included on a potential energy surface that contains these interactions 37 (DFT-BLYP+vdW), which are known to be necessary for the description of H-bonded systems. 31The underlying reason is depicted in the free energy plots in Fig. 5(a), as also discussed in Ref. 71.The oxygen-oxygen distance dOO in such wires is smaller when including vdW interactions such that the barrier at the proton transfer coordinate ν is similar to kBT at room temperature, even without including NQEs.The hydrogen transfer is, thus, not the rate limiting step for proton diffusion at this temperature, which lies instead on the rearrangement of the water molecules in the wire.When disregarding vdW interactions, the hydrogen transfer barrier is higher than kBT at room temperature for classical nuclei and likely becomes the rate limiting step.However, when considering quantum nuclei, ZPE makes the effective barrier much smaller, and this ceases to be the case.
A final example is that of the vibrational spectrum of porphycene, presented in Ref. 23 and depicted in Fig. 5(b).The internal hydrogen atoms in this molecule exhibit a strongly anharmonic potential energy profile along the hydrogen-transfer coordinate characterized by a double well.Among several DFT functionals considered in Ref. 23, B3LYP with pairwise vdW corrections (B3LYP+vdW) was the one that best approximated the relative energies of the stationary points in the potential energy surface calculated with coupled-cluster including single, double, and perturbative triple excitations [CCSD(T)].The calculation of the anharmonic infrared (IR) spectrum of porphycene with this functional shows that the NH-stretch signal, located in the region between 2000 and 3000 cm −1 , is broadened and red-shifted when approximating NQEs with TRPMD.The TRPMD result matches experiment very well, 23 and a strengthening of the H-bonds was observed upon including NQEs.In contrast, the PBE functional with many-body dispersion (MBD) interactions 40,96 predicts much lower barriers for the hydrogen transfer in this system compared to CCSD(T).In fact, they are so low that a harmonic estimate of ZPE showed that they would completely "fill" the double well, making the barrier irrelevant, as discussed in Ref. 23.Indeed, when calculating the IR spectrum with this functional from classicalnuclei MD, the NH-stretch peak is already considerably red-shifted (almost centered at the position that B3LYP-TRPMD predicts).When including NQEs, it becomes extremely broad and overly red-shifted, in complete disagreement with experiments.Anharmonic effects in porphycene are so pronounced that intermode coupling leads to a qualitatively different temperature dependence on the position of this same peak if nuclei are considered as quantum or classical particles. 77This phenomenon is bound to be observed in a wide range of high-dimensional anharmonic systems with significant coupling between high and low frequency vibrational modes, which can now be addressed with the appropriate methodology.It is also coupling between vibrational modes of the surface and adsorbate that results in a surprising tunneling rate enhancement of The Journal of Chemical Physics PERSPECTIVE scitation.org/journal/jcpintramolecular hydrogen transfer when porphycene is adsorbed on metallic surfaces. 63

VI. FINAL REMARKS
Including NQEs beyond the harmonic approximation in ab initio simulations is within reach in many applications in chemical physics, condensed matter physics, and materials science.As these areas become more focused on soft organic materials for technological applications, these methods will undoubtedly deliver the much needed and unprecedented fundamental insights into many problems where both electrons and nuclei require a quantummechanical description.They are also relevant for understanding matter and chemistry under extreme conditions as found in interstellar space or the interior of planets.
At the same time, the recent methodological developments discussed here have pushed the boundaries where the real challenges in these simulations lie.In this Perspective, a few methodological challenges were mentioned, namely, the need for practical and more accurate methods that can capture dynamical NQEs and the need of overcoming the barrier of performing these simulations with more accurate electronic structure methods.The former still seems slightly elusive, but the field is quite active and there is hope for new solutions.The latter, on the other hand, can be addressed by taking advantage of the recent developments in machine-learning techniques for atomistic simulations, coupled with the possibility of reusing data stored in large databases.
An important fundamental next step is the development of a practical methodology that is able to capture electron-phonon coupling in non-adiabatic scenarios and where many anharmonic (quantum) nuclear degrees of freedom play an important rolepossibly also in non-equilibrium situations.The open questions in this area are considerably numerous but will hopefully also profit from the advances achieved so far for the adiabatic regime.

FIG. 1 .
FIG. 1. (a)A sketch depicting the conventional model of atoms under the Born-Oppenheimer approximation.A nuclear configuration provides the commonly called "external potential" (V ext ) that the electrons are subject to, while the ground state expectation value of the electronic Hamiltonian provides the Born-Oppenheimer potential (V BO ), which enters as the potential energy of the nuclear degrees of freedom.Under this approximation, it is possible to treat the nuclei as quantum or classical particles (b).In turn, in both cases, one can assume that the forces acting on the nuclei are approximately harmonic, or not, being described by the gradient with respect to nuclear positions q of the Born-Oppenheimer potential.

FIG. 2 .
FIG. 2. (a)Example of a two-dimensional potential energy surface represented by a quartic potential along q 1 and q 2 with quadratic coupling between the coordinates.(b) Corresponding harmonic approximation of the potential energy surface depicted in (a) around its two lowest energy minima in solid colors overlaid on the original surface.The approximation distorts the potential around the minima, completely neglecting the coupling between the two coordinates.

FIG. 3 .
FIG. 3.An example of different regimes of nuclear fluctuations.(a) The anharmonic contribution to the quantum and classical nuclear fluctuations as a function of temperature and oscillation frequency, defined as Δ(T, ω) = 1 − ⟨q 2 (T, ω)⟩ harm./⟨q 2 (T, ω)⟩ anh. .The anharmonic expectation values were calculated based on a quartic polynomial expansion of a one-dimensional Morse potential V(q) = D(1 − e −αq ) 2 , with α 2 = mω 2 /(2D).Parameters for the sketch are m = mp and D = 0.1 Ha.Two viewpoints of the same plot are presented.The region of coincidence of quantum and classical anharmonic contributions on the frequency axis increases with the increase in temperature.Intervals shown in the plot span roughly up to 3500 cm −1 and up to 600 K. (b) The difference ΔΔ(T, ω) = Δ qt (T, ω) − Δ cl (T, ω), where the anharmonic potential of (a) was modified by a scaling term (1 − ω/ωmax) that multiplies the third and fourth orders of the polynomial expansion (potential becomes harmonic approaching ω = ωmax).Only a region at intermediate frequencies and lower/intermediate temperatures exhibits higher degrees of "quantum anharmonicity."

FIG. 4 .
FIG. 4. (a) Estimation of anharmonic contributions to forces according to the quantity ε defined in the text for all Cartesian degrees of freedom of a water molecule adsorbed on a Pd(111) surface with four layers, ordered as x, y, and z for each atom.The layers are labeled with different colors, and L4 is in contact with the molecule.Dark blue circles indicate ε calculated from classical Boltzmann statistical sampling, while green squares were calculated from quantum Boltzmann statistical sampling for the degrees of freedom.(b) Same quantities as in (a) but on the normal mode representation, showing only the modes located mostly on the adsorbate, pictured at the bottom.(c) Change in the work function (DFT-PBE+vdW surf ) of water molecules adsorbed on the steps of a Pt(221) surface, as shown in Ref. 85.The molecular dipole component in the direction of the surface is larger if quantum statistics for the nuclei are considered, which causes a different work function change.(d) Change in the work function (DFT-PBE+vdW surf ) of cyclohexane adsorbed on Rh(111), as shown in Ref. 34.The molecule lies slightly closer to the surface if quantum statistics for the nuclei are considered, which causes a difference in work function change.
FIG. 5. (a) Analysis of the dynamics of diffusion of a proton in a water wire confined by a cylindrical potential, as presented in Ref. 71.The centers of the Boys orbitals of each water monomer as calculated from DFT-BLYP(+vdW), shown in black at the top, were used as descriptors of a machine-learning clustering algorithm to identify the position of the proton at any point in time (green translucent sphere at the top).Free energy surfaces projected on the oxygen-oxygen distances d OO and the proton transfer coordinate ν = d OH1 − d OH2 , as calculated from MD and PIMD simulations with the BLYP and the BLYP+vdW functionals.These data were presented in Ref. 71.(b) Infrared spectrum of the porphycene molecule as calculated from MD and TRPMD calculations (T = 300 K) on potential energy surfaces calculated with the B3LYP+vdW functional and the PBE+MBD functional.Overlaid snapshots of the B3LYP+vdW TRPMD simulation are shown in the upper panel.These data were presented in Ref. 23.