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. The reasons for this impressive progress can be attributed to a number of contributing factors, including the increase in the availability of high performance computers, improved algorithms in many community software packages,^{1–3} and the simpatico relationship between machine learning and quantum mechanical methods.^{4–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.^{8} Despite 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.^{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 “first-principles” 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–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 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 *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 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,17}

Within the realm of DFT, one can obtain very reliable results by employing good approximations to the exchange–correlation potential, *ν*_{xc} = *δE*_{xc}[*ρ*]/*δρ*, 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–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.^{25} This 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.^{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–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 (also hybrid) functionals.^{36–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.^{24} The 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

*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 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,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 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. In keeping with a quantum picture of nuclei, one route is to propose harmonic potentials that effectively capture anharmonic terms.^{45–48} Another route is to consider higher-order terms in a perturbative Taylor series expansion of the potential^{48–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 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.

*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*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.”

^{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

*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*/

*ℏ*, and $qj(i)$ represents the position degree of freedom

*i*in the

*j*th 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,

*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 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.^{61} In 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,64}

^{65–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–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 $c\xafAB$ in RPMD/TRPMD is given by

**(**

*q***) is the vector of all 3**

*p**N*Cartesian positions (momenta) of the system. The time evolution of the position and momenta in TRPMD is given by

*ξ*

_{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**

*γ***with elements $Kjj\u2032=\omega P2(2\delta j,j\u2032\u2212\delta j,j\u2032\u22121\u2212\delta j,j\u2032+1)$ and its transformation to the internal modes of the free ring polymer $K\u0303=CTKC=\omega \u03032$, where the transformation matrices**

*K***are given in Ref. 74 and $\omega \u0303kk\u2032=2\omega P\u2061sin(k\pi /P)\delta kk\u2032$. The original scheme proposes $\gamma \u0303=\omega \u0303$ ($\gamma =C\gamma \u0303CT$), which ensures optimal damping of the potential. In addition, $\gamma \u030300=0$ for the centroid mode, which ensures that the formalism maintains all quantum-mechanical limits that the RPMD autocorrelation functions were shown to obey.**

*C*^{67}Because 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.

^{76}The approach is applicable particularly 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 artifacts of RPMD and centroid molecular dynamics on vibrational spectra,

^{78}albeit introducing an unphysical broadening to the line shape 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 such as RPMD, TRPMD, and centroid molecular dynamics (CMD)^{65} 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?

*λ*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, *ℏω*/(*k*_{B}*T*), 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.^{77}

Beyond 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 *T*_{c}, below which nuclear tunneling will be pronounced. Assuming a parabolic barrier, *T*_{c} = *ℏω*^{†}/(2*πk*_{B}),^{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.^{83}

One 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,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 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).

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 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.^{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 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 k-points, 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 $\Delta i=\u2211sS(Fs,iDFT\u2212Fs,ih)2/S$ was calculated, where $Fs,iDFT/h$ 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 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,92} However, 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.^{93} When 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,95}

Further 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.^{71} A depiction of the location of the proton in one simulation snapshot is shown on the top of Fig. 5(a). 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 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.^{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 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 *k*_{B}*T* 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 classical-nuclei 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.^{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 the appropriate methodology. It is also coupling between vibrational modes of the 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 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 quantum-mechanical 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 role—possibly 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.

## ACKNOWLEDGMENTS

I acknowledge the members of my research group, who have greatly inspired me throughout the past years. In the spirit of this special issue, I thank the numerous mentors and collaborators and especially also thank 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.

## REFERENCES

*ab initio*accuracy to 100 million atoms with machine learning

*F*center in molten KCl

^{+}

_{5}

*Ab initio*path integral molecular dynamics: Basic ideas

_{n}-LysH

^{+}polyalanine peptides (n = 5,10,15) in Vacuo: Helical or not?

*ab initio*liquid water: The interplay of nuclear and electronic quantum effects

*ab initio*parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu

*Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra*

*Statistical Mechanics*

*Solid State Physics*

_{3}NH

_{3})

*Quantum Mechanics and Path Integrals*

*Ab Initio*Molecular Dynamics: Basic Theory and Advanced Methods

^{+}and OH

^{−}diffusion along confined water wires

*Ab initio*molecular simulations with numeric atom-centered orbitals