{'fraction_non_alphanumeric': 0.026411148131034314, 'fraction_numerical': 0.008949633989153433, 'mean_word_length': 5.1226920786182255, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
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 (NQE) beyond a harmonic or perturbative treatment, the impact of NQE on electronic properties of weakly-bonded systems, and how different first-principles potential energy surfaces can change the impact of NQE 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. The reasons for this impressive progress can be attributed to a number of conspiring factors including the increase of the availability of high performance computers, improved algorithms in many community software packages [1][2][3] , and the sympatico relationship of machine learning with quantum mechanical methods [4][5][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 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 sketched in Fig. 1(a).
The need for a quantum mechanical treatment of electrons has long been well accepted 8 . Despite this fact, classical-like functional forms (or, force fields) that do not explicitly treat the electronic degrees of freedom individually but rather 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 9 . These 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 "firstprinciples" methods restores the transferability, while improving the accuracy, of the overall approach.
On the other hand, the assumption that the nuclei behave classically in most situations has been left essentially unchallenged for decades in many areas of physical chemistry and materials science, with few notable exceptions [10][11][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, 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, as shown in (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.
where the quantum nature of the nuclei provokes 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 corrections within the harmonic approximation is a standard practice, for example, when analysing relative and reaction energies in the gas-phase and in some solid-state problems. However, treating nuclei as quantum particles still represents 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 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. It does not cover any techniques going beyond the adiabatic approximation for electronic and nuclear degrees of freedom. The perspective includes considerations to estimate whether NQE are relevant for a specific problem, and how the shape of the potential energy surface can influence the degree to which NQE 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 V ext , 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 NQE on diverse properties of materials.
The BO potentia canl 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 last decades. Because of these algorithmic developments, very accurate methodologies, like coupled cluster theory, have become available to treat larger systems including periodic solids 16,17 .
Within the realm of DFT, one can obtain very reliable results by employing good approximations to the exchangecorrelation potential, ν xc = δ E xc [ρ]/δ ρ, where ρ is the electronic density. Even though all known approximations to ν xc are not exact, they yield, in many cases, reliable electronic properties for molecules, liquids and solids. In addition, they can also provide an accurate potential energy surface, from which structural and dynamical properties of the nuclei can be calculated [18][19][20][21][22][23] , which attests 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 25 . This family of so-called "hybrid" functionals is more computationally demanding, but recent implementations of these methods are very efficient, and allow the treatment of either thousands of atoms 26 or thousands of force evaluations for smaller systems 27 . Tailor-made and local hybrid functionals capable of delivering good accuracy for interfaces have also been recently developed 28,29 .
For weakly bonded systems, capturing long range dispersion interactions (usually simply called van der Waals interactions) is of extreme importance [30][31][32][33][34] . These interactions are highly non-local and thus absent from standard DFT functionals, due to their inherent spatial locality. They can be accounted by non-local functionals 35 , or by corrections that are added to semilocal (and hybrid) functionals [36][37][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 for 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 24 . The lack of an exact density functional means that a very large amount of exchange-correlation approximations exist and it is often difficult to evaluate whether a particular one can yield predictive results for a given (largescale) 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 has 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 q i 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 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, as for example the evaluation of zero point energy, free energies, vibrational spectra, reaction rates, and others 42,43 . In 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 44 . However, what this picture completely disregards is the coupling between different vibrational modes, as well as anharmonic terms pertaining to any given mode. In Figure 2(a,b) a pictorial representation of a harmonic approximation to a considerably anharmonic 2dimensional 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. Keeping with a quantum picture of nuclei, one route is to propose harmonic potentials that effectively capture anharmonic terms [45][46][47][48] . Another 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 54 . However, 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 are 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-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 classical Hamiltonian and thus behave as classical particles. Instead, the method of ab initio path integral molecular dynamics 10,12 , which was pioneered in the 80s and has gained broad popularity in the past decade, allows 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. 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" 58 . Under 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/(Pk B T ), ω P = Pk B T /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 When P tends to infinity, sampling configurations with the equations of motion generated by H P 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 path integral molecular dynamics 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 progresses being achieved in the last 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 realtime propagator, can overcome the sign problem 61 . In particular, the ring-polymer instanton method 62 , which uses trajectories in imaginary time to describe nuclear tunnelling along well-defined reaction coordinates, has been successfully applied to quite high-dimensional systems, including molecules adsorbed on metallic surfaces 63,64 .
Attempts 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. Nevertheless, these methods have met considerable success in the simulation of vibrational spectroscopy, diffusion properties, reaction rates, and others 21,23,[68][69][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 functionsc AB in RPMD/TRPMD is given bȳ . . dq P dp 1 . . . dp P × A P (q(0))B P (q(t))e −β P H P where A P (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 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 K j j = ω 2 P (2δ j, j − δ j, j −1 − δ j, j +1 ) 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 quantum-mechanical limits that RPMD autocorrelation functions were shown to obey 67 . Because the formalism imposes only a few restrictions on the friction coefficients, other white-noise friction parameters inspired by harmonicwell problems have been proposed 75 , and the extension of the formalism to generalized Langevin equation (GLE) thermostats has also been explored 76 . The approach is applicable, in particular, to vibrational spectra of condensed phase systems 19,21 and molecules with many degrees of freedom 23,77 at moderate temperatures. It is free from the spurious artefacts of RPMD and centroid molecular dynamics on vibrational spectra 78 , albeit introducing an unphysical broadening to the lineshape of the peaks 67 . This broadening could be partially mitigated with GLE thermostats 76 . At a higher simulation cost, path integral Liouville dynamics can deliver excellent results for the vibrational spectra of small molecules 68 .
The 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 like RPMD, TRPMD and CMD can be obtained from it 75,80 . This 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?
NQE 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 NQE. 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 like Cu, λ Cu = 0.12 Å. However, at a much lower temperature like 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 depend 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 dynamical properties of materials like 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 for these properties is the ratio between the harmonic estimate of the zero point energy of vibrational states and their thermal occupation,hω/(k B T ), which, if much larger than 1 hints at the (possibly large) impact of NQE. Further, 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 77 .
Beyond effects related to the zero-point-energy, 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 T c , below which nuclear tunneling will be pronounced. Assuming a parabolic barrier, T c =hω † /(2πk B ) 82 , where ω † is the absolute value of the imaginary frequency at the transition state geometry. Deviations from a parabolic barrier, however, leads to cases where this crossover temperature ceases to be a good estimate for tunneling contributions 83 .
One could continue listing a few of these "rules-of-thumb", like the estimation of the Debye temperature for heat conduction studies, and others. However, what is clear is that these 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". estimations are based on predictions that employ harmonic potentials, as results can be obtained analytically. Indeed, NQE are relatively easy to account for when most degrees of freedom behave harmonically or only exhibit small deviations from harmonic behavior. This fact underlies a number of longstanding successes in the field solid-state physics 51,84 . Challenges 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 the previous section 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. 3a, 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 increasing temperature. Furthermore, while in the classical case the anharmonic contributions to these fluctuations are very 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. 3b, 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). 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).
For static properties of materials, it is indeed often found that even when anharmonic effects are pronounced, they can be well-captured 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 like 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 89 . Therefore, 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,86 .
However, there are situations where the harmonic approximation for NQE 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 like hydrogen play a prominent role. This is exemplified in Figure 4a for one water molecule adsorbed on a 2×2 Pd(111) surface with 4 layers (FHI-aims program light settings 90 , 6×6×1 k-points, 100 Å vacuum). The anharmonic contributions to the forces of this system was 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 PBE functional including the vdW interactions of Ref. 91 . Subsequently, the root mean square 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 Fig. 4(a) and (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 NQE, in particular for the modes that bend the hydrogen atoms towards the surface.
These anharmonic NQE 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. Another striking example of these effects can be seen in a related system of water molecules adsorbed on the steps of a Pt(221) slab. The work function change ∆φ of the interface was calculated with quantum and classical nuclei 85 by performing ab initio path integral molecular dynamics using an acceleration scheme that reduces the necessary amount of PIMD replicas in weakly bonded inter- faces, proposed in Ref. 85 . As reported in that work, a significant difference in ∆φ is introduced when considering quantum nuclei, as reproduced in Fig 4(b). The origin of this effect was explained to be the increase in magnitude of the molecular dipole component pointing towards the surface upon the inclusion of NQEs. In contrast, for a largely non-polar molecule, namely cyclohexane adsorbed on Rh(111), NQE were shown to impact the ∆φ both in experiment and theory 34,92 . However, in this case the origin of the effect does not lie on anharmonic NQE on the internal degrees of freedom of the adsorbate. As depicted in Fig. 4(c), 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 adsorbate and surface and the pushback effect. Therefore, a variety of consequences of the interplay between anharmonic NQE and electronic structure can emerge at interfaces.
The potential energy surface itself plays a decisive role in how NQE 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 93 . When analysing the impact of NQE in dynamical properties of water with PIMDbased 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 NQE on this quantity is small. However, this is not because NQE in water are negligible. Instead, NQE 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 NQE 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 NQE on certain properties. For liquid water, excellent agreement with experiment for diffusion coefficients and vibrational properties have 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, like biomolecules 94,95 .
Further examples where changing the ab initio potential can provoke substantial changes of the impact of NQE on dynamical properties of different systems is shown in Fig. 5. In Ref. 71 the diffusion of proton 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 71 . A depiction of the location of the proton in one simulation snap-  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. This data was 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. This data was presented in Ref. 23 . shot is shown on the top of Fig. 5a. Using this technique, it was observed that the diffusion coefficient of a proton in a water wire at 300 K is almost doubled when considering NQE on a potential energy surface that does not contain long-range van der Waals interactions (DFT-BLYP). However, the diffusion remains basically unchanged when NQE 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 31 . The underlying reason is depicted in the free energy plots in Fig. 5(a), as also discussed in Ref. 71 . The oxygen-oxygen distance d OO in such wires is smaller when including vdW interactions, such that the barrier at the proton transfer coordinate ν is similar to k B T at room temperature, even without including NQE. 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 k B T at room temperature for classical nuclei and likely becomes the rate limiting step. However, when considering quantum nuclei, zero-point energy 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. 5b. The in-ternal 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 NQE. In contrast, the PBE functional with many-body dispersion interactions 40,96 (MBD) 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 classical-nuclei MD, the NH-stretch peak is already considerably red-shifted (almost centered at the position that B3LYP-TRPMD predicts). When including NQE, 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 77 . This 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 appropriate methodology. It is also coupling between vibrational modes of surface and adsorbate that results in a surprising tunneling rate enhancement of intramolecular hydrogen transfer when porphycene is adsorbed on metallic surfaces 63 .

VI. FINAL REMARKS
Including NQE 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 quantum-mechanical description. They are also relevant for understanding matter and chemistry in extreme conditions as found, e.g., 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 NQE, 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 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 role -possibly also in non-equilibrium situations. The open questions in this area are considerably more numerous, but will hopefully also profit from the advances achieved so far for the adiabatic regime.

ACKNOWLEDGEMENTS
I would like to acknowledge the members of my research group, who have greatly inspired me throughout the past years. In the spirit of this special issue, I would also like to thank the numerous mentors, collaborators, and especially female mentors, collaborators, and colleagues who work in my area. All of them have given me strength and inspiration to continue on this path. I also thank Michele Ceriotti and Aaron Kelly for reading earlier versions of this perspective and for very useful discussions.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.