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.

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

FIG. 1.

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

FIG. 1.

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

Close modal

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.

Under the Born–Oppenheimer approximation, the electronic structure problem is solved at a fixed configuration of the nuclei, which defines the external potential Vext, and the electrons create the Born–Oppenheimer (BO) potential VBO 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 Physics15 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 = δExc[ρ]/δρ, where ρ is the electronic density. Even though all known approximations to νxc are not exact, they yield, in many cases, good electronic properties for molecules, liquids, and solids. In addition, they can provide an accurate potential energy surface, from which structural and dynamical properties of the nuclei can be calculated,18–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 atoms26 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 functionals35 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 methods41 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.

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,
VBOqVBO(q0)+12i,j3NqiVBO(q)qiqj0qj,
(1)
where qi are nuclear coordinates, N is the number of atoms in the system, and 0 denotes a global or local minimum of the BO potential. The second term on the right of Eq. (1) involves the second derivative of the potential with respect to nuclear positions (the Hessian) and can be written in diagonal form by a coordinate transformation. The solution of such a Schrödinger equation for the nuclei yields the nuclear vibrational modes and their respective energies. These solutions lead to several analytical results for diverse quantities, i.e., the evaluation of ZPE, free energies, vibrational spectra, and reaction rates.42,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.
FIG. 2.

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

FIG. 2.

(a) Example of a two-dimensional potential energy surface represented by a quartic potential along q1 and q2 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.

Close modal

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 potential48–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.

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 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
Z=limPZP=limPi=13Nm(i)2π2βP3P2i=13Ndq1dqP×ExpβPi=13Nj=1Pm(i)ωP22(qj+1(i)qj(i))2+j=1PVBO(qj(1),,qj(3N)),
(2)
where P is a convergence parameter that corresponds to the discretization of imaginary time, or the commonly called “beads” of a ring polymer, βP = 1/(PkBT), ωP = PkBT/, and qj(i) 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,
HP=j=1Pi=13N(pj(i))22m(i)+12m(i)ωP2(qj(i)qj1(i))2+j=1PVBOqj(1),qj(2),,qj(3N).
(3)
When P tends to infinity, sampling configurations with the equations of motion generated by HP gives access to the quantum Boltzmann distribution of particles in the system and hence guarantees the correct quantum statistics in the limit of distinguishable particles. The masses m′ entering the kinetic energy term in Eq. (3) can assume any value, as they only serve to assist with sampling.

Unlike ab initio molecular dynamics with classical nuclei, where the dynamics algorithm is both a tool for statistical sampling and a means of real time propagation of nuclear coordinates, in PIMD, the time evolution of the equations of motion serves only as a sampling tool. It is possible to formulate real-time path integral methods to access the quantum real-time evolution of nuclei. However, even with notable progress being achieved in the past years,60 the applicability of these methods remains rather limited due to the sign problem arising from the oscillatory nature of the real-time propagator. Semiclassical methods, based on different semiclassical approximations of the real-time propagator, can overcome the sign problem.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

Attempts to approximate quantum dynamics with techniques based on imaginary-time path integrals65–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¯AB in RPMD/TRPMD is given by
c¯AB(t)limP1(2π)PZPdq1dqPdp1dpP×AP(q(0))BP(q(t))eβPHP,
(4)
where AP(q)=j=1PA(qj)/P and q (p) is the vector of all 3N Cartesian positions (momenta) of the system. The time evolution of the position and momenta in TRPMD is given by
qj(i)t=pj(i)/m(i),
(5)
pj(i)t=iVBO(qj(1)qj(3N))jm(i)ωP2(2δj,jδj,j1δj,j+1)qj(i)jγjjpj(i)+2m(i)βPjγ12jjξj(i),
(6)
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 Kjj=ωP2(2δj,jδj,j1δj,j+1) and its transformation to the internal modes of the free ring polymer K̃=CTKC=ω̃2, where the transformation matrices C are given in Ref. 74 and ω̃kk=2ωPsin(kπ/P)δkk. The original scheme proposes γ̃=ω̃ (γ=Cγ̃CT), 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 the 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 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 systems19,21 and molecules with many degrees of freedom23,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 dynamics79 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 methods81 and offers potential avenues for further developments.

NQEs on equilibrium and response properties are more pronounced at lower temperatures and for lighter particles. A quantity that illustrates this particularly well is the thermal de Broglie wavelength, given by
λ=2π2mkBT12.
(7)
Even though this expression assumes non-interacting free particles, it tends to be a good estimator of the importance of NQEs. Classical Boltzmann statistics typically break down if λ is of the same order of, or larger than, the typical interparticle spacing. As an example, for a proton at room temperature, λH = 1.00 Å, while for a heavy nucleus such as Cu, λCu = 0.12 Å. However, at a much lower temperature such as 10 K, these values are λH = 5.50 Å and λCu = 0.69 Å.

One telltale sign of the importance of NQE is the appearance of geometric isotope effects. Classically, the mass of an atom cannot change static structural properties of materials, but within quantum mechanics, it can. However, the magnitude of these changes depends on the character of the potential energy surface. If the local potential energy is approximately harmonic, the effect will be negligibly small. However, if it is not, the effect can be large.

Nevertheless, even classically, the atomic mass changes the dynamical properties of materials, such as diffusion and transport coefficients or vibrational spectra. In a harmonic potential, quantum and classical mechanics would predict the same mass scalings for vibrational frequencies, but the occupation of the vibrational states would differ. A useful quantity to analyze these properties is the ratio between the harmonic estimate of the ZPE of vibrational states and their thermal occupation, ℏω/(kBT), which, if much larger than 1, hints at the (possibly large) impact of NQEs. Furthermore, if the potential is strongly anharmonic, the dependence of the vibrational energies on the mass of the atoms can be very different in quantum and classical mechanics.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 Tc, below which nuclear tunneling will be pronounced. Assuming a parabolic barrier, Tc = ℏω/(2πkB),82 where ω is the absolute value of the imaginary frequency at the transition state geometry. Deviations from a parabolic barrier, however, lead to cases where this crossover temperature ceases to be a good estimate for tunneling contributions.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.

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 ⟨q2⟩ 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).

FIG. 3.

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

FIG. 3.

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

Close modal

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 Δi=sS(Fs,iDFTFs,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.

FIG. 4.

(a) Estimation of anharmonic contributions to forces according to the quantity ε defined in the text for all Cartesian degrees of freedom of a water molecule adsorbed on a Pd(111) surface with four layers, ordered as x, y, and z for each atom. The layers are labeled with different colors, and L4 is in contact with the molecule. Dark blue circles indicate ε calculated from classical Boltzmann statistical sampling, while green squares were calculated from quantum Boltzmann statistical sampling for the degrees of freedom. (b) Same quantities as in (a) but on the normal mode representation, showing only the modes located mostly on the adsorbate, pictured at the bottom. (c) Change in the work function (DFT-PBE+vdWsurf) of water molecules adsorbed on the steps of a Pt(221) surface, as shown in Ref. 85. The molecular dipole component in the direction of the surface is larger if quantum statistics for the nuclei are considered, which causes a different work function change. (d) Change in the work function (DFT-PBE+vdWsurf) of cyclohexane adsorbed on Rh(111), as shown in Ref. 34. The molecule lies slightly closer to the surface if quantum statistics for the nuclei are considered, which causes a difference in work function change.

FIG. 4.

(a) Estimation of anharmonic contributions to forces according to the quantity ε defined in the text for all Cartesian degrees of freedom of a water molecule adsorbed on a Pd(111) surface with four layers, ordered as x, y, and z for each atom. The layers are labeled with different colors, and L4 is in contact with the molecule. Dark blue circles indicate ε calculated from classical Boltzmann statistical sampling, while green squares were calculated from quantum Boltzmann statistical sampling for the degrees of freedom. (b) Same quantities as in (a) but on the normal mode representation, showing only the modes located mostly on the adsorbate, pictured at the bottom. (c) Change in the work function (DFT-PBE+vdWsurf) of water molecules adsorbed on the steps of a Pt(221) surface, as shown in Ref. 85. The molecular dipole component in the direction of the surface is larger if quantum statistics for the nuclei are considered, which causes a different work function change. (d) Change in the work function (DFT-PBE+vdWsurf) of cyclohexane adsorbed on Rh(111), as shown in Ref. 34. The molecule lies slightly closer to the surface if quantum statistics for the nuclei are considered, which causes a difference in work function change.

Close modal

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 nuclei85 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 reported21 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 interactions37 (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 dOO in such wires is smaller when including vdW interactions such that the barrier at the proton transfer coordinate ν is similar to kBT at room temperature, even without including NQEs. The hydrogen transfer is, thus, not the rate limiting step for proton diffusion at this temperature, which lies instead on the rearrangement of the water molecules in the wire. When disregarding vdW interactions, the hydrogen transfer barrier is higher than kBT at room temperature for classical nuclei and likely becomes the rate limiting step. However, when considering quantum nuclei, ZPE makes the effective barrier much smaller, and this ceases to be the case.

FIG. 5.

(a) Analysis of the dynamics of diffusion of a proton in a water wire confined by a cylindrical potential, as presented in Ref. 71. The centers of the Boys orbitals of each water monomer as calculated from DFT-BLYP(+vdW), shown in black at the top, were used as descriptors of a machine-learning clustering algorithm to identify the position of the proton at any point in time (green translucent sphere at the top). Free energy surfaces projected on the oxygen–oxygen distances dOO and the proton transfer coordinate ν = dOH1dOH2, as calculated from MD and PIMD simulations with the BLYP and the BLYP+vdW functionals. These data were presented in Ref. 71. (b) Infrared spectrum of the porphycene molecule as calculated from MD and TRPMD calculations (T = 300 K) on potential energy surfaces calculated with the B3LYP+vdW functional and the PBE+MBD functional. Overlaid snapshots of the B3LYP+vdW TRPMD simulation are shown in the upper panel. These data were presented in Ref. 23.

FIG. 5.

(a) Analysis of the dynamics of diffusion of a proton in a water wire confined by a cylindrical potential, as presented in Ref. 71. The centers of the Boys orbitals of each water monomer as calculated from DFT-BLYP(+vdW), shown in black at the top, were used as descriptors of a machine-learning clustering algorithm to identify the position of the proton at any point in time (green translucent sphere at the top). Free energy surfaces projected on the oxygen–oxygen distances dOO and the proton transfer coordinate ν = dOH1dOH2, as calculated from MD and PIMD simulations with the BLYP and the BLYP+vdW functionals. These data were presented in Ref. 71. (b) Infrared spectrum of the porphycene molecule as calculated from MD and TRPMD calculations (T = 300 K) on potential energy surfaces calculated with the B3LYP+vdW functional and the PBE+MBD functional. Overlaid snapshots of the B3LYP+vdW TRPMD simulation are shown in the upper panel. These data were presented in Ref. 23.

Close modal

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) interactions40,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 

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.

1.
V. W.-z.
Yu
,
C.
Campos
,
W.
Dawson
,
A.
García
,
V.
Havu
,
B.
Hourahine
,
W. P.
Huhn
,
M.
Jacquelin
,
W.
Jia
,
M.
Keçeli
,
R.
Laasner
,
Y.
Li
,
L.
Lin
,
J.
Lu
,
J.
Moussa
,
J. E.
Roman
,
Á.
Vázquez-Mayagoitia
,
C.
Yang
, and
V.
Blum
, “
ELSI—An open infrastructure for electronic structure solvers
,”
Comput. Phys. Commun.
256
,
107459
(
2020
).
2.
M. J. T.
Oliveira
,
N.
Papior
,
Y.
Pouillon
,
V.
Blum
,
E.
Artacho
,
D.
Caliste
,
F.
Corsetti
,
S.
de Gironcoli
,
A. M.
Elena
,
A.
García
,
V. M.
García-Suárez
,
L.
Genovese
,
W. P.
Huhn
,
G.
Huhs
,
S.
Kokott
,
E.
Küçükbenli
,
A. H.
Larsen
,
A.
Lazzaro
,
I. V.
Lebedeva
,
Y.
Li
,
D.
López-Durán
,
P.
López-Tarifa
,
M.
Lüders
,
M. A. L.
Marques
,
J.
Minar
,
S.
Mohr
,
A. A.
Mostofi
,
A.
O’Cais
,
M. C.
Payne
,
T.
Ruh
,
D. G. A.
Smith
,
J. M.
Soler
,
D. A.
Strubbe
,
N.
Tancogne-Dejean
,
D.
Tildesley
,
M.
Torrent
, and
V. W.-z.
Yu
, “
The CECAM electronic structure library and the modular software development paradigm
,”
J. Chem. Phys.
153
,
024117
(
2020
).
3.
V.
Kapil
,
M.
Rossi
,
O.
Marsalek
,
R.
Petraglia
,
Y.
Litman
,
T.
Spura
,
B.
Cheng
,
A.
Cuzzocrea
,
R. H.
Meißner
,
D. M.
Wilkins
,
B. A.
Helfrecht
,
P.
Juda
,
S. P.
Bienvenue
,
W.
Fang
,
J.
Kessler
,
I.
Poltavsky
,
S.
Vandenbrande
,
J.
Wieme
,
C.
Corminboeuf
,
T. D.
Kühne
,
D. E.
Manolopoulos
,
T. E.
Markland
,
J. O.
Richardson
,
A.
Tkatchenko
,
G. A.
Tribello
,
V.
Van Speybroeck
, and
M.
Ceriotti
, “
i-PI 2.0: A universal force engine for advanced molecular simulations
,”
Comput. Phys. Commun.
236
,
214
223
(
2019
).
4.
R.
Jinnouchi
,
J.
Lahnsteiner
,
F.
Karsai
,
G.
Kresse
, and
M.
Bokdam
, “
Phase transitions of hybrid perovskites simulated by machine-learning force fields trained on the fly with Bayesian inference
,”
Phys. Rev. Lett.
122
,
225701
(
2019
).
5.
W.
Jia
,
H.
Wang
,
M.
Chen
,
D.
Lu
,
J.
Liu
,
L.
Lin
,
R.
Car
,
W.
E
, and
L.
Zhang
, “
Pushing the limit of molecular dynamics with ab initio accuracy to 100 million atoms with machine learning
,” arXiv:2005.00223 (
2020
).
6.
S. J.
Buxton
and
S.
Habershon
, “
Accelerated path-integral simulations using ring-polymer interpolation
,”
J. Chem. Phys.
147
,
224107
(
2017
).
7.
S.
Chmiela
,
H. E.
Sauceda
,
I.
Poltavsky
,
K.-R.
Müller
, and
A.
Tkatchenko
, “
sGDML: Constructing accurate and data efficient molecular force fields using machine learning
,”
Comput. Phys. Commun.
240
,
38
45
(
2019
).
8.
N.
Bohr
, “
I. On the constitution of atoms and molecules
,”
London, Edinburgh Dublin Philos. Mag. J. Sci.
26
,
1
25
(
1913
).
9.
D. E.
Shaw
,
P.
Maragakis
,
K.
Lindorff-Larsen
,
S.
Piana
,
R. O.
Dror
,
M. P.
Eastwood
,
J. A.
Bank
,
J. M.
Jumper
,
J. K.
Salmon
,
Y.
Shan
, and
W.
Wriggers
, “
Atomic-level characterization of the structural dynamics of proteins
,”
Science
330
,
341
346
(
2010
).
10.
M.
Parrinello
and
A.
Rahman
, “
Study of an F center in molten KCl
,”
J. Chem. Phys.
80
,
860
867
(
1984
).
11.
D.
Marx
and
M.
Parrinello
, “
Structural quantum effects and three-centre two-electron bonding in CH+5
,”
Nature
375
,
216
218
(
1995
).
12.
D.
Marx
and
M.
Parrinello
, “
Ab initio path integral molecular dynamics: Basic ideas
,”
J. Chem. Phys.
104
,
4077
4082
(
1995
).
13.
J. A.
Poulsen
,
G.
Nyman
, and
P. J.
Rossky
, “
Practical evaluation of condensed phase quantum correlation functions: A Feynman–Kleinert variational linearized path integral method
,”
J. Chem. Phys.
119
,
12179
(
2003
).
14.
T. E.
Markland
and
M.
Ceriotti
, “
Nuclear quantum effects enter the mainstream
,”
Nat. Rev. Chem.
2
,
0109
(
2018
).
15.
C. D.
Sherrill
,
D. E.
Manolopoulos
,
T. J.
Martínez
, and
A.
Michaelides
, “
Electronic structure software
,”
J. Chem. Phys.
153
,
070401
(
2020
).
16.
C.
Riplinger
and
F.
Neese
, “
An efficient and near linear scaling pair natural orbital based local coupled cluster method
,”
J. Chem. Phys.
138
,
034106
(
2013
).
17.
I. Y.
Zhang
and
A.
Grüneis
, “
Coupled cluster theory in materials science
,”
Front. Mater.
6
,
123
(
2019
).
18.
M.
Rossi
,
V.
Blum
,
P.
Kupser
,
G.
von Helden
,
F.
Bierau
,
K.
Pagel
,
G.
Meijer
, and
M.
Scheffler
, “
Secondary structure of Ac-Alan-LysH+ polyalanine peptides (n = 5,10,15) in Vacuo: Helical or not?
,”
J. Phys. Chem. Lett.
1
,
3465
3470
(
2010
).
19.
M.
Rossi
,
H.
Liu
,
F.
Paesani
,
J.
Bowman
, and
M.
Ceriotti
, “
Communication: On the consistency of approximate quantum dynamics simulation methods for vibrational spectra in the condensed phase
,”
J. Chem. Phys.
141
,
181101
(
2014
).
20.
N.
Raimbault
,
V.
Athavale
, and
M.
Rossi
, “
Anharmonic effects in the low-frequency vibrational modes of aspirin and paracetamol crystals
,”
Phys. Rev. Mater.
3
,
053605
(
2019
).
21.
O.
Marsalek
and
T. E.
Markland
, “
Quantum dynamics and spectroscopy of ab initio liquid water: The interplay of nuclear and electronic quantum effects
,”
J. Phys. Chem. Lett.
8
,
1545
1551
(
2017
).
22.
X.-Z.
Li
,
B.
Walker
, and
A.
Michaelides
, “
Quantum nature of the hydrogen bond
,”
Proc. Natl. Acad. Sci. U. S. A.
108
,
6369
6373
(
2011
).
23.
Y.
Litman
,
J. O.
Richardson
,
T.
Kumagai
, and
M.
Rossi
, “
Elucidating the nuclear quantum dynamics of intramolecular double hydrogen transfer in porphycene
,”
J. Am. Chem. Soc.
141
,
2526
2534
(
2019
).
24.
A. J.
Cohen
,
P.
Mori-Sánchez
, and
W.
Yang
, “
Challenges for density functional theory
,”
Chem. Rev.
112
,
289
320
(
2011
).
25.
A. V.
Krukau
,
G. E.
Scuseria
,
J. P.
Perdew
, and
A.
Savin
, “
Hybrid functionals with local range separation
,”
J. Chem. Phys.
129
,
124103
(
2008
).
26.
S. M.
Janke
,
M.
Rossi
,
S. V.
Levchenko
,
S.
Kokott
,
M.
Scheffler
, and
V.
Blum
, “
Pentacene and tetracene molecules and films on H/Si(111): Level alignment from hybrid density functional theory
,”
Electron. Struct.
2
,
035002
(
2020
).
27.
B.
Cheng
,
J.
Behler
, and
M.
Ceriotti
, “
Nuclear quantum effects in water at the triple point: Using theory as a link between experiments
,”
J. Phys. Chem. Lett.
7
,
2210
2215
(
2016
).
28.
P.
Borlido
,
M. A. L.
Marques
, and
S.
Botti
, “
Local hybrid density functional for interfaces
,”
J. Chem. Theory Comput.
14
,
939
947
(
2018
).
29.
H.
Zheng
,
M.
Govoni
, and
G.
Galli
, “
Dielectric-dependent hybrid functionals for heterogeneous materials
,”
Phys. Rev. Mater.
3
,
073803
(
2019
).
30.
E. R.
Johnson
,
S.
Keinan
,
P.
Mori-Sánchez
,
J.
Contreras-García
,
A. J.
Cohen
, and
W.
Yang
, “
Revealing noncovalent interactions
,”
J. Am. Chem. Phys. Soc.
132
,
6498
6506
(
2010
).
31.
N.
Marom
,
A.
Tkatchenko
,
M.
Rossi
,
V. V.
Gobre
,
O.
Hod
,
M.
Scheffler
, and
L.
Kronik
, “
Dispersion interactions with density-functional theory: Benchmarking semiempirical and interatomic pairwise corrected density functionals
,”
J. Chem. Theory Comput.
7
,
3944
3951
(
2011
).
32.
M.
Rossi
,
A.
Tkatchenko
,
S. B.
Rempe
, and
S.
Varma
, “
Role of methyl-induced polarization in ion binding
,”
Proc. Natl. Acad. Sci. U. S. A.
110
,
12978
12983
(
2013
).
33.
A.
Tkatchenko
,
M.
Rossi
,
V.
Blum
,
J.
Ireta
, and
M.
Scheffler
, “
Unraveling the stability of polypeptide helices: Critical role of van der Waals interactions
,”
Phys. Rev. Lett.
106
,
118102
(
2011
).
34.
K.
Fidanyan
,
I.
Hamada
, and
M.
Rossi
, “
Quantum nuclei at weakly bonded interfaces: The case of cyclohexane on Rh(111)
,”
Adv. Theory Simul.
2000241
(published online
2021
).
35.
K.
Berland
,
V. R.
Cooper
,
K.
Lee
,
E.
Schröder
,
T.
Thonhauser
,
P.
Hyldgaard
, and
B. I.
Lundqvist
, “
van der Waals forces in density functional theory: A review of the vdW-DF method
,”
Rep. Prog. Phys.
78
,
066501
(
2015
).
36.
A. D.
Becke
and
E. R.
Johnson
, “
A density-functional model of the dispersion interaction
,”
J. Chem. Phys.
123
,
154101
(
2005
).
37.
A.
Tkatchenko
, “
Accurate molecular van der Waals interactions from ground-state electron density and free-atom reference data
,”
Phys. Rev. Lett.
102
,
073005
(
2009
).
38.
S.
Grimme
,
J.
Antony
,
S.
Ehrlich
, and
H.
Krieg
, “
A consistent and accurate ab initio parametrization of density functional dispersion correction (DFT-D) for the 94 elements H-Pu
,”
J. Chem. Phys.
132
,
154104
(
2010
).
39.
T.
Gould
,
S.
Lebègue
,
J. G.
Ángyán
, and
T.
Bučko
, “
A fractionally ionic approach to polarizability and van der Waals many-body dispersion calculations
,”
J. Chem. Theory Comput.
12
,
5920
5930
(
2016
).
40.
A.
Tkatchenko
,
R. A.
Distasio
,
R.
Car
, and
M.
Scheffler
, “
Accurate and efficient method for many-body van der Waals interactions
,”
Phys. Rev. Lett.
108
,
236402
(
2012
).
41.
L.
Reining
, “
The GW approximation: Content, successes and limitations
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
8
,
e1344
(
2018
).
42.
E. B.
Wilson
,
J. C.
Decius
, and
P. C.
Cross
,
Molecular Vibrations: The Theory of Infrared and Raman Vibrational Spectra
, Dover Books on Chemistry (
Dover Publications
,
New York, NY
,
1980
).
43.
D. A.
McQuarrie
,
Statistical Mechanics
, Harper’s Chemistry Series (
Harper & Row
,
New York, NY
,
1976
).
44.
N. W.
Ashcroft
and
N. D.
Mermin
,
Solid State Physics
(
Holt Rinehart and Winston
,
New York, NY
,
1976
).
45.
O.
Hellman
,
I. A.
Abrikosov
, and
S. I.
Simak
, “
Lattice dynamics of anharmonic solids from first principles
,”
Phys. Rev. B
84
,
180301
(
2011
).
46.
S. E.
Brown
,
I.
Georgescu
, and
V. A.
Mandelshtam
, “
Self-consistent phonons revisited. II. A general and efficient method for computing free energies and vibrational spectra of molecules and clusters
,”
J. Chem. Phys.
138
,
044317
(
2013
).
47.
I.
Errea
,
M.
Calandra
, and
F.
Mauri
, “
Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides
,”
Phys. Rev. B
89
,
064302
(
2014
).
48.
V.
Kapil
,
E.
Engel
,
M.
Rossi
, and
M.
Ceriotti
, “
Assessment of approximate methods for anharmonic free energies
,”
J. Chem. Theory Comput.
15
,
5845
5857
(
2019
).
49.
Q.
Yu
and
J. M.
Bowman
, “
Vibrational second-order perturbation theory (VPT2) using local monomer normal modes
,”
Mol. Phys.
113
,
3964
3971
(
2015
).
50.
C.
König
, “
Tailored multilevel approaches in vibrational structure theory: A route to quantum mechanical vibrational spectra for complex systems
,”
Int. J. Quantum Chem.
121
,
e26375
(
2021
).
51.
B.
Monserrat
,
N. D.
Drummond
, and
R. J.
Needs
, “
Anharmonic vibrational properties in periodic systems: Energy, electron-phonon coupling, and stress
,”
Phys. Rev. B
87
,
144302
(
2013
).
52.
L. S.
Norris
,
M. A.
Ratner
,
A. E.
Roitberg
, and
R. B.
Gerber
, “
Møller–Plesset perturbation theory applied to vibrational problems
,”
J. Chem. Phys.
105
,
11261
11267
(
1998
).
53.
B. R.
Westbrook
,
E. M.
Valencia
,
S. C.
Rushing
,
G. S.
Tschumper
, and
R. C.
Fortenberry
, “
Anharmonic vibrational frequencies of ammonia borane (BH3NH3)
,”
J. Chem. Phys.
154
,
041104
(
2021
).
54.
J.
Hoja
,
H.-Y.
Ko
,
M. A.
Neumann
,
R.
Car
,
R. A.
DiStasio
, and
A.
Tkatchenko
, “
Reliable and practical computational description of molecular crystal polymorphs
,”
Sci. Adv.
5
,
eaau3338
(
2019
).
55.
R.
Car
and
M.
Parrinello
, “
Unified approach for molecular dynamics and density-functional theory
,”
Phys. Rev. Lett.
55
,
2471
2474
(
1985
).
56.
T. D.
Kühne
, “
Second generation Car–Parrinello molecular dynamics
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
4
,
391
406
(
2014
).
57.
A. M. N.
Niklasson
, “
Next generation extended Lagrangian first principles molecular dynamics
,”
J. Chem. Phys.
147
,
054103
(
2017
).
58.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
, International Series in Pure and Applied Physics (
McGraw-Hill
,
New York, NY
,
1965
).
59.
D.
Marx
and
J.
Hutter
,
Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods
(
Cambridge University Press
,
Cambridge
,
2009
).
60.
N.
Makri
, “
Modular path integral methodology for real-time quantum dynamics
,”
J. Chem. Phys.
149
,
214108
(
2018
).
61.
W. H.
Miller
, “
The semiclassical initial value representation: A potentially practical way for adding quantum effects to classical molecular dynamics simulations
,”
J. Phys. Chem. A
105
,
2942
2955
(
2001
).
62.
J. O.
Richardson
, “
Perspective: Ring-polymer instanton theory
,”
J. Chem. Phys.
148
,
200901
(
2018
).
63.
Y.
Litman
and
M.
Rossi
, “
Multidimensional hydrogen tunneling in supported molecular switches: The role of surface interactions
,”
Phys. Rev. Lett.
125
,
216001
(
2020
).
64.
W.
Fang
,
J.
Chen
,
P.
Pedevilla
,
X.-Z.
Li
,
J. O.
Richardson
, and
A.
Michaelides
, “
Origins of fast diffusion of water dimers on surfaces
,”
Nat. Commun.
11
,
1689
(
2020
).
65.
J.
Cao
and
G. A.
Voth
, “
The formulation of quantum statistical mechanics based on the Feynman path centroid density. III. Phase space formalism and analysis of centroid molecular dynamics
,”
J. Chem. Phys.
101
,
6157
6167
(
1994
).
66.
I. R.
Craig
and
D. E.
Manolopoulos
, “
Quantum statistics and classical mechanics: Real time correlation functions from ring polymer molecular dynamics
,”
J. Chem. Phys.
121
,
3368
(
2004
).
67.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
, “
How to remove the spurious resonances from ring polymer molecular dynamics
,”
J. Chem. Phys.
140
,
234116
(
2014
).
68.
J.
Liu
and
Z.
Zhang
, “
Path integral Liouville dynamics: Applications to infrared spectra of OH, water, ammonia, and methane
,”
J. Chem. Phys.
144
,
034307
(
2016
).
69.
G. R.
Medders
and
F.
Paesani
, “
Dissecting the molecular structure of the air/water interface from quantum simulations of the sum-frequency generation spectrum
,”
J. Am. Chem. Soc.
138
,
3912
3919
(
2016
).
70.
S.
Habershon
,
D. E.
Manolopoulos
,
T. E.
Markland
, and
T. F.
Miller
III
, “
Ring-polymer molecular dynamics: Quantum effects in chemical dynamics from classical trajectories in an extended phase space
,”
Annu. Rev. Phys. Chem.
64
,
387
413
(
2013
).
71.
M.
Rossi
,
M.
Ceriotti
, and
D. E.
Manolopoulos
, “
Nuclear quantum effects in H+ and OH diffusion along confined water wires
,”
J. Phys. Chem. Lett.
7
,
3001
3007
(
2016
).
72.
J. R.
Cendagorta
,
A.
Powers
,
T. J. H.
Hele
,
O.
Marsalek
,
Z.
Bačić
, and
M. E.
Tuckerman
, “
Competing quantum effects in the free energy profiles and diffusion rates of hydrogen and deuterium molecules through clathrate hydrates
,”
Phys. Chem. Chem. Phys.
18
,
32169
32177
(
2016
).
73.
K. A.
Jung
,
P. E.
Videla
, and
V. S.
Batista
, “
Ring-polymer, centroid, and mean-field approximations to multi-time Matsubara dynamics
,”
J. Chem. Phys.
153
,
124112
(
2020
).
74.
M.
Ceriotti
,
M.
Parrinello
,
T. E.
Markland
, and
D. E.
Manolopoulos
, “
Efficient stochastic thermostatting of path integral molecular dynamics
,”
J. Chem. Phys.
133
,
124104
(
2010
).
75.
T. J. H.
Hele
, “
On the relation between thermostatted ring-polymer molecular dynamics and exact quantum dynamics
,”
Mol. Phys.
114
,
1461
1471
(
2016
).
76.
M.
Rossi
,
V.
Kapil
, and
M.
Ceriotti
, “
Fine tuning classical and quantum molecular dynamics using a generalized Langevin equation
,”
J. Chem. Phys.
148
,
102301
(
2018
).
77.
Y.
Litman
,
J.
Behler
, and
M.
Rossi
, “
Temperature dependence of the vibrational spectrum of porphycene: A qualitative failure of classical-nuclei molecular dynamics
,”
Faraday Discuss.
221
,
526
546
(
2019
).
78.
A.
Witt
,
S. D.
Ivanov
,
M.
Shiga
,
H.
Forbert
, and
D.
Marx
, “
On the applicability of centroid and ring polymer path integral molecular dynamics for vibrational spectroscopy
,”
J. Chem. Phys.
130
,
194510
(
2009
).
79.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
, “
Boltzmann-conserving classical dynamics in quantum time-correlation functions: ‘Matsubara dynamics
,’”
J. Chem. Phys.
142
,
134103
(
2015
).
80.
T. J. H.
Hele
,
M. J.
Willatt
,
A.
Muolo
, and
S. C.
Althorpe
, “
Communication: Relation of centroid molecular dynamics and ring-polymer molecular dynamics to exact quantum dynamics
,”
J. Chem. Phys.
142
,
191101
(
2015
).
81.
G.
Trenins
,
M. J.
Willatt
, and
S. C.
Althorpe
, “
Path-integral dynamics of water using curvilinear centroids
,”
J. Chem. Phys.
151
,
054109
(
2019
).
82.
M. J.
Gillan
, “
Quantum-classical crossover of the transition rate in the damped double well
,”
J. Phys. C: Solid State Phys.
20
,
3621
3641
(
1987
).
83.
S.
Álvarez-Barcia
,
J. R.
Flores
, and
J.
Kästner
, “
Tunneling above the crossover temperature
,”
J. Phys. Chem. A
118
,
78
82
(
2013
).
84.
F.
Giustino
, “
Electron–phonon interactions from first principles
,”
Rev. Mod. Phys.
89
,
015003
(
2017
).
85.
Y.
Litman
,
D.
Donadio
,
M.
Ceriotti
, and
M.
Rossi
, “
Decisive role of nuclear quantum effects on surface mediated water dissociation at finite temperature
,”
J. Chem. Phys.
148
,
102320
(
2017
).
86.
M.
Rossi
,
P.
Gasparotto
, and
M.
Ceriotti
, “
Anharmonic and quantum fluctuations in molecular crystals: A first-principles study of the stability of paracetamol
,”
Phys. Rev. Lett.
117
,
115702
(
2016
).
87.
V.
Kapil
,
J.
Wieme
,
S.
Vandenbrande
,
A.
Lamaire
,
V.
Van Speybroeck
, and
M.
Ceriotti
, “
Modeling the structural and thermal properties of loaded metal–organic frameworks. An interplay of quantum and anharmonic fluctuations
,”
J. Chem. Theory Comput.
15
,
3237
3249
(
2019
).
88.
R.
Ramírez
,
N.
Neuerburg
,
M.-V.
Fernández-Serra
, and
C. P.
Herrero
, “
Quasi-harmonic approximation of thermodynamic properties of ice Ih, II, and III
,”
J. Chem. Phys.
137
,
044502
(
2012
).
89.
F.
Knoop
,
T. A. R.
Purcell
,
M.
Scheffler
, and
C.
Carbogno
, “
Anharmonicity measure for materials
,”
Phys. Rev. Mater.
4
,
083809
(
2020
).
90.
V.
Blum
,
R.
Gehrke
,
F.
Hanke
,
P.
Havu
,
V.
Havu
,
X.
Ren
,
K.
Reuter
, and
M.
Scheffler
, “
Ab initio molecular simulations with numeric atom-centered orbitals
,”
Comput. Phys. Commun.
180
,
2175
2196
(
2009
).
91.
V. G.
Ruiz
,
W.
Liu
,
E.
Zojer
,
M.
Scheffler
, and
A.
Tkatchenko
, “
Density-functional theory with screened van der Waals interactions for the modeling of hybrid inorganic–organic systems
,”
Phys. Rev. Lett.
108
,
146103
(
2012
).
92.
T.
Koitaya
and
J.
Yoshinobu
, “
The quantum nature of C–H⋯metal interaction: Vibrational spectra and kinetic and geometric isotope effects of adsorbed cyclohexane
,”
Chem. Rec.
14
,
848
856
(
2014
).
93.
S.
Habershon
,
T. E.
Markland
, and
D. E.
Manolopoulos
, “
Competing quantum effects in the dynamics of a flexible water model
,”
J. Chem. Phys.
131
,
024501
(
2009
).
94.
M.
Rossi
,
W.
Fang
, and
A.
Michaelides
, “
Stability of complex biomolecular structures: van der Waals, hydrogen bond cooperativity, and nuclear quantum effects
,”
J. Phys. Chem. Lett.
6
,
4233
4238
(
2015
).
95.
W.
Fang
,
J.
Chen
,
M.
Rossi
,
Y.
Feng
,
X.-Z.
Li
, and
A.
Michaelides
, “
Inverse temperature dependence of nuclear quantum effects in DNA base pairs
,”
J. Phys. Chem. Lett.
7
,
2125
2131
(
2016
).
96.
R. A.
DiStasio
, Jr.
,
V. V.
Gobre
, and
A.
Tkatchenko
, “
Many-body van der Waals interactions in molecules and condensed matter
,”
J. Phys.: Condens. Matter
26
,
213202
(
2014
).
Published open access through an agreement withMax Planck Institute for the Structure and Dynamics of Matter