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

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.⁸ 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.⁹ 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,¹⁴ 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.

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¹⁵ 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,²⁴ and range-separated functionals can, moreover, approximate electronic screening.²⁵ 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²⁶ or thousands of force evaluations for smaller systems.²⁷ 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³⁵ 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.²⁴ 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⁴¹ 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.

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.⁵⁴ 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.

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,⁶⁰ 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.⁶¹ In particular, the ring-polymer instanton method,⁶² 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

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⁷⁹ 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)⁶⁵ can be obtained from it.^75,80 This theory has already opened the door to the improvement of some of these approximate methods⁸¹ and offers potential avenues for further developments.

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_BT), 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.⁷⁷ 

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),⁸² 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.⁸³ 

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.

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²⟩ 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.⁸⁹ 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,⁹⁰ 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 $Δi=∑sS(Fs,iDFT−Fs,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⁸⁵ 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.⁹³ 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²¹ 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.⁷¹ 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³⁷ (DFT-BLYP+vdW), which are known to be necessary for the description of H-bonded systems.³¹ 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_BT 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_BT 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⁻¹, is broadened and red-shifted when approximating NQEs with TRPMD. The TRPMD result matches experiment very well,²³ 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.⁷⁷ 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.⁶³ 

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.

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.