The incorporation of nuclear quantum effects and non-Born–Oppenheimer behavior into quantum chemistry calculations and molecular dynamics simulations is a longstanding challenge. The nuclear–electronic orbital (NEO) approach treats specified nuclei, typically protons, quantum mechanically on the same level as the electrons with wave function and density functional theory methods. This approach inherently includes nuclear delocalization and zero-point energy in molecular energy calculations, geometry optimizations, reaction paths, and dynamics. It can also provide accurate descriptions of excited electronic, vibrational, and vibronic states as well as nuclear tunneling and nonadiabatic dynamics. Nonequilibrium nuclear–electronic dynamics simulations beyond the Born–Oppenheimer approximation can be used to investigate a wide range of excited state processes. This Perspective provides an overview of the foundational NEO methods and enumerates the prospects for using these methods as building blocks for future developments. The conceptual simplicity and computational efficiency of the NEO approach will enhance its accessibility and applicability to diverse chemical and biological systems.
A wide range of chemical and biological processes rely on the quantum mechanical properties of nuclei, such as zero-point energy effects, vibrationally excited states, and nuclear tunneling. Many of these processes also involve a breakdown of the Born–Oppenheimer separation between electrons and nuclei. A prime example is proton-coupled electron transfer (PCET), which is essential for photosynthesis, respiration, DNA synthesis, and many electrocatalytic reactions.1–3 A meaningful theoretical description of PCET reactions often requires the quantum mechanical treatment of the transferring hydrogen nuclei as well as the electrons and a rigorous treatment of hydrogen tunneling and nonadiabatic effects.1 Vibronically nonadiabatic PCET rate constants and kinetic isotope effects can be calculated in terms of nonadiabatic transitions between electron–proton vibronic states.4 Moreover, simulating the nonequilibrium dynamics of PCET reactions typically entails nonadiabatic molecular dynamics methods.
Incorporating nuclear quantum effects and non-Born–Oppenheimer behavior into quantum chemistry calculations and molecular dynamics simulations is challenging. In conventional Born–Oppenheimer approaches, the potential energy surface is generated by solving the electronic Schrödinger equation for fixed nuclear configurations. The nuclei can be propagated either classically or quantum mechanically (i.e., with wavepacket or path integral methods) on the resulting adiabatic potential energy surface. In some cases, the potential energy surface can be generated on-the-fly with ab initio electronic structure methods during the nuclear dynamics. Nonadiabatic dynamics approaches, such as the Ehrenfest,5 surface hopping,6 multiple spawning,7,8 multiconfigurational time-dependent Hartree,9,10 and nonadiabatic ring polymer molecular dynamics11–14 methods, allow nuclei to move on multiple potential energy surfaces. These nonadiabatic dynamics methods invoke different types of approximations and include nuclear quantum effects to varying degrees.15
In contrast to these methods, the nuclear–electronic orbital (NEO) approach16,17 treats specified nuclei quantum mechanically on the same level as the electrons using molecular orbital techniques without invoking the Born–Oppenheimer separation between the electrons and the quantum nuclei. In this approach, the mixed nuclear–electronic time-independent or time-dependent Schrödinger equation is solved with wave function or density functional theory (DFT) methods. Typically, the quantum nuclei are protons, although other nuclei and other types of particles, such as positrons, can also be treated quantum mechanically within this framework. At least two nuclei are usually treated classically to avoid problems with translations or rotations, which alternatively could be removed or mitigated while treating all nuclei quantum mechanically.18
An advantage of the NEO approach is that nuclear quantum effects such as proton delocalization and zero-point energy, as well as non-Born–Oppenheimer effects between electrons and protons, are inherently included during geometry optimizations, reaction paths, and dynamics. The non-Born–Oppenheimer effects between the quantum and classical nuclei, as well as the electrons and classical nuclei, have been found to be negligible for equilibrium molecular properties, as determined by analysis of the NEO diagonal Born–Oppenheimer corrections.19 For nonequilibrium processes, these non-Born–Oppenheimer effects can be incorporated through nonadiabatic dynamics methods, such as Ehrenfest and surface hopping approaches for the classical nuclei moving on the NEO electron–proton vibronic surfaces. A variety of related multicomponent wave function and DFT methods have been explored,20–27 and exact factorization methods28,29 provide an alternative promising solution to these problems. This brief overview does not cover the vast array of non-Born–Oppenheimer or nonadiabatic methods, and more comprehensive reviews are provided elsewhere.15,17
This Perspective briefly reviews the foundational NEO methods and provides a road map for how these methods can be used as building blocks for future developments and applications. Electron–proton correlation is particularly important because of the attractive electron–proton Coulomb interaction.16,17 Two different paths have been explored for including correlation in NEO calculations of ground states: wave function methods, such as coupled cluster (CC),30 and multicomponent DFT.31–33 These paths have been extended to equation-of-motion (EOM)34,35 and time-dependent DFT (TDDFT)36,37 methods for the description of excited electronic, vibrational, and vibronic states. Furthermore, the combination of the NEO real-time TDDFT38 and Ehrenfest dynamics39,40 methods allows the simulation of nonequilibrium nuclear–electronic dynamics beyond the Born–Oppenheimer approximation. Finally, the NEO multistate DFT (NEO-MSDFT) method41 enables the description of hydrogen tunneling systems. The capacity to describe hydrogen tunneling and nonequilibrium nonadiabatic dynamics in a conceptually straightforward and computationally efficient manner opens up many avenues for exploration.
In the NEO framework, the system is divided into electrons, quantum nuclei, and classical nuclei. Herein, the quantum nuclei are assumed to be protons, but extensions to other types of nuclei and even positrons are straightforward. The NEO Hamiltonian, HNEO, includes the kinetic energies of the electrons and quantum protons, the interactions of the electrons and quantum protons with the classical nuclei, and the electron–electron, proton–proton, and electron–proton Coulomb interaction terms. The simplest NEO method is based on the NEO-Hartree–Fock (NEO-HF) wave function,16 which is the product of electronic and protonic Slater determinants,
Here, xe and xp are the collective spatial and spin coordinates of the electrons and quantum protons, respectively. The electronic and protonic Slater determinants, Φe(xe) and Φp(xp), are composed of electronic and protonic spin orbitals, which are expanded in electronic and protonic Gaussian basis sets.42
Minimizing the Hartree–Fock energy with respect to the electronic and protonic spin orbitals leads to two sets of coupled Hartree–Fock–Roothaan equations for the electrons and protons,16
Here, Fe, Ce, Se, and ɛe are the electronic Fock matrix, orbital coefficient matrix, overlap matrix, and orbital energy matrix, respectively, and the protonic matrices are defined analogously. These two equations are strongly coupled because the electronic and protonic Fock matrices each depend on both the electronic and protonic orbital coefficients. As a result, these two equations must be solved self-consistently. Typically each quantum proton is represented by a set of electronic and protonic basis sets with the same center, which is optimized variationally during the self-consistent field (SCF) procedure. Due to the lack of electron–proton correlation, the NEO-HF method produces proton densities that are much too localized, preventing the accurate description of molecular properties such as geometries, energies, and frequencies.16,17 Nevertheless, it serves as the starting point for the correlated NEO methods.
MINIMUM ENERGY PATHS IN THE NEO FRAMEWORK
The implementation of analytical NEO gradients and Hessians43 has enabled the optimization of geometries and transition states, as well as the generation of minimum energy paths (MEPs). The NEO potential energy surface depends on only the classical nuclear coordinates, with the underlying assumption that the basis function centers associated with the quantum protons are optimized variationally for each configuration on this potential energy surface. The MEP is generated by starting at a transition state on the NEO potential energy surface, defined to be a first-order saddle point, and using a steepest descent algorithm to step down to the reactant and product minima, optimizing the quantum proton basis function center positions at each step to remain on the NEO potential energy surface. Figure 1 depicts a hydride transfer reaction between two carbon atoms and the corresponding MEP computed with the conventional electronic HF and NEO-HF methods.43 For this system, the NEO-HF energy barrier is lower than the conventional HF energy barrier because the zero-point energy of the transferring hydrogen nucleus is inherently included in the NEO energy. The analytical NEO Hessian, which depends only on the classical nuclear coordinates, enables the inclusion of the zero-point energy and entropic contributions from the classical nuclei. The entropic contributions from the quantum nuclei can be included with the vibrational analysis methods discussed below,43,44 thereby providing thermochemical properties that include the anharmonic effects of the quantum nuclei.
Analysis of the single imaginary normal mode at the transition state, as well as the geometries along the MEP, provides insight into the dominant nuclear motions driving this hydride transfer reaction. In the conventional electronic HF case, the imaginary normal mode at the transition state is dominated by the motion of the transferring hydrogen, which also contributes significantly to the intrinsic reaction coordinate (IRC) along the MEP. In the NEO-HF case, this hydrogen nucleus is treated quantum mechanically and therefore cannot contribute to the normal mode or the IRC. Instead, the imaginary normal mode at the transition state is dominated by the tetrahedral-to-planar rearrangement about the two central carbon atoms as they switch between sp3 and sp2 hybridizations. The geometries along the MEP confirm that this tetrahedral-to-planar rearrangement represents the dominant contribution to the IRC and thus drives the hydride transfer reaction. Interestingly, this same motion also contributes to the imaginary mode at the transition state for the conventional electronic HF case, even though it is overwhelmed by the transferring hydrogen motion. If the contribution from the transferring hydrogen is eliminated from the conventional HF imaginary normal mode, the resulting vector is nearly identical to the NEO-HF imaginary normal mode after renormalization.
Thus, the NEO MEP provides fundamental insights into the nuclear motions that drive hydrogen transfer reactions, including proton, hydride, and hydrogen atom transfer, as well as PCET. This type of analysis is analogous to analyses of electron transfer reactions that identify the nuclear motions driving electron transfer. In this case, however, the nuclear motions are driving the transfer of a quantum mechanical hydrogen nucleus as well as the rearrangement of electrons. These nuclear motions correspond to the inner-sphere reorganization contributing to the hydrogen transfer reaction. Furthermore, analysis of the reactive electronic intrinsic bond orbitals45,46 and protonic orbitals along the MEP provides insight into whether the electron and proton are transferring synchronously or asynchronously. For the specific reaction depicted in Fig. 1, such an analysis indicates that this hydride transfer reaction entails synchronous electron and proton transfer. These types of mechanistic insights are useful for understanding and tuning chemical reactions.
MULTICOMPONENT COUPLED CLUSTER AND CONFIGURATION INTERACTION METHODS
Correlation effects can be included within the NEO framework with the multicomponent configuration interaction (CI) and coupled cluster (CC) approaches.23,30,47 Analogous to the conventional electronic counterparts,48,49 the CI and CC wave functions can be expressed as
Here, the cluster excitation operator generates all single, double, triple, and further excited determinants among the electrons and protons by acting on the NEO-HF reference state. For example, at the singles and doubles (SD) level, includes electronic and protonic single excitations, electronic and protonic double excitations, and mixed electronic–protonic double excitations.
The NEO-CISD and NEO-CCSD methods, which include all single and double excitations, have been derived and implemented.30 Moreover, the NEO-BCCD50 and NEO-OOCCD51 methods, which use the coupled cluster with doubles ansatz in conjunction with optimized Brueckner orbitals or variationally optimized orbitals, respectively, have also been implemented. The NEO-CISD methods are not sufficiently accurate due to the absence of orbital relaxation.30 In contrast, the NEO-CCSD, NEO-BCCD, and NEO-OOCCD methods have been shown to provide accurate proton densities, proton affinities, and optimized geometries for the chemical systems studied.30,50,51 More recently, the NEO unitary coupled cluster (UCC) ansatz has been combined with the variational quantum eigensolver algorithm and formulated in the qubit basis to enable calculations of ground state energies and wave functions of multicomponent systems on quantum computers.52
An illustrative example is the application of the NEO-CCSD method to protonated water tetramers, treating all nine protons quantum mechanically.53 This application became computationally feasible by combining the NEO-CCSD method with density fitting53,54 for approximating the four-center two-particle integrals. The inclusion of zero-point energy contributions has been shown to change the relative stabilities of the four protonated water tetramer isomers at the conventional electronic CCSD with perturbative triples [CCSD(T)] level of theory, as shown in Fig. 2 (black and blue symbols).55 This effect is observed when the zero-point energies are computed via diagonalization of the mass-weighted Hessian matrix within the harmonic approximation and when anharmonic corrections are included with the vibrational second-order perturbative (VPT2) approach.53,55,56 NEO-CCSD single-point energy calculations on the optimized CCSD(T) geometries predict the correct ordering of the isomer energies,53 as shown in Fig. 2 (purple symbols). The NEO-CCSD method inherently includes the zero-point energy contributions from the quantized protons, and the contributions from the oxygen atoms to the relative zero-point energies of these isomers have been shown to be negligible.53 Thus, the NEO-CCSD method provides accurate results for this system without requiring the computationally expensive calculation of the Hessian matrix.
The NEO orbital optimized second-order Møller–Plesset perturbation theory (NEO-OOMP2) method51 offers a more computationally efficient alternative that approaches the accuracy of the NEO-CCSD method. Additional computational efficiency arises from the scaled-opposite-spin (SOS) strategy,57 where the opposite-spin and same-spin components of the second-order electron–electron perturbative term are scaled differently, thereby offering the opportunity to eliminate the computationally expensive same-spin component. The multicomponent extension of this strategy,51 denoted SOS′, also applies a scaling factor to the second-order electron–proton perturbative term. Although the NEO-OOMP2 method is stable and qualitatively reasonable for the systems studied, the NEO-SOS′-OOMP2 method is advantageous in that it produces proton densities and proton affinities that are nearly as accurate as the NEO-CCSD results.51 Moreover, the NEO-SOS′-OOMP2 method can be implemented in a manner that leads to N4 scaling rather than the N6 scaling of NEO-CCSD, where N is a measure of the system size.
The NEO-CCSD and NEO-SOS′-OOMP2 methods provide the foundation for additional developments of NEO ground state wave function methods. The implementation of analytical gradients and Hessians will enable geometry optimizations and the calculation of zero-point energy contributions, thermochemical properties, and MEPs at these correlated levels of theory. An advantage of the NEO-CCSD method is that it can be improved systematically by including higher-order excitations, starting with triple excitations. However, this approach will become computationally prohibitive for larger systems, even with density fitting methods. In these cases, the NEO-SOS′-OOMP2 method will be a more practical option, and reparameterization of the scaling factors for specific types of systems could provide additional flexibility. Embedding methods58 in which only a small portion of the system is treated at the NEO-CCSD level provide another viable option.
MULTICOMPONENT EQUATION-OF-MOTION METHODS FOR EXCITED STATES
Analogous to conventional electronic equation-of-motion CCSD (EOM-CCSD) methods,59,60 the NEO-EOM-CCSD method34 has been implemented for calculating excited electronic, protonic, and mixed electron–proton vibronic excitation energies in both the frequency and time domains.35 Figure 3 illustrates the application of this approach to the HCN molecule, where all electrons and the hydrogen nucleus are treated quantum mechanically. This approach produces qualitatively reasonable energies and intensities for the fundamental proton vibrational excitations (i.e., the hydrogen stretch and bend), as well as the combination band corresponding to simultaneous excitation of both bending modes and the overtone associated with the bending mode, as depicted in Fig. 3 (left side). In addition, this method can describe double excitations corresponding to the simultaneous excitation of an electron and a proton, such as an excited proton vibrational state within an excited electronic state, as shown in Fig. 3 (right side).35 However, the energies of these double excitations are much too high, presumably mainly due to the truncated coupled cluster ansatz.
An advantage of the NEO-EOM methods is that they can be improved systematically by including higher-order excitations, starting with triples, but the computational expense becomes prohibitive with current computational capabilities. Again, embedding methods may render the NEO-EOM methods computationally viable. Moreover, these methods play an important role for benchmarking other excited state NEO methods. Additionally, a NEO-EOM method has been formulated in the qubit basis to enable the quantum computation of excitation energies for multicomponent systems.52
In multicomponent DFT,22,25,61,62 more than one type of particle is treated quantum mechanically. Here, we discuss multicomponent DFT for electrons and protons. In the multicomponent Hohenberg–Kohn theorems,63 the total energy is a functional of the one-particle electronic and protonic densities. In the multicomponent Kohn–Sham formalism,22,25,61,62 the reference wave function is the product of electronic and protonic Slater determinants composed of electronic and protonic Kohn–Sham orbitals, respectively, as in Eq. (1). Optimization of the NEO-DFT energy with respect to the electronic and protonic orbitals leads to two sets of Kohn–Sham equations for the electrons and protons, respectively, as in Eq. (2). These equations are solved self-consistently, as described above for the NEO-HF method.
Within the NEO-DFT framework, the total energy functional is expressed as
where ρe and ρp are the electronic and protonic densities, respectively. Here, Eext[ρe, ρp] is the interaction of the electronic and protonic densities with the external potential created by the classical nuclei, Eref[ρe, ρp] includes the noninteracting kinetic energies of the electrons and quantum protons, as well as the classical electron–electron, proton–proton, and electron–proton Coulomb energies, and the last three terms correspond to different types of functionals. The electron–electron exchange–correlation functional, Eexc[ρe], can be chosen from the many available electronic functionals.64 The proton–proton exchange–correlation functional, Epxc[ρp], is typically chosen to be the Hartree–Fock exchange energy. For molecular systems, the proton densities are highly localized, and therefore, proton–proton exchange and correlation energies are negligible,17 although the diagonal exchange terms must be included to eliminate self-interaction terms. The electron–proton correlation functional, Eepc[ρe, ρp], has required the development of new types of functionals.
The epc family of functionals was developed by extending the Colle–Salvetti formalism65 to multicomponent systems. The multicomponent electron–proton wave function ansatz is the product of electronic and protonic wave functions that include all electron–electron and proton–proton exchange–correlation effects, respectively, multiplied by an explicit electron–proton correlation factor. This electron–proton correlation factor depends on an inverse correlation length parameter, which is chosen to be the geometric mean of the inverse Wigner–Seitz radii for an electron and a proton for the epc17 and epc19 functionals32,33,66 and the arithmetic mean of these radii for the epc18 functional.64 The epc17 functional has the following form:32
where a, b, and c are parameters that can be fit to properties of interest.
The epc17 and epc18 functionals32,33,64 rely on the local density approximation (LDA) and depend on only the electron and proton densities. They have been shown to be of similar accuracy, and the epc17 functional has been adopted for most applications because it is mathematically simpler. The epc19 functional66 relies on the generalized gradient approximation (GGA) and depends on the gradients of the electronic and protonic densities as well as the densities themselves. Thus, the epc19 functional is analogous to the well-known LYP functional67 for electron–electron correlation. The epc17 and epc19 functionals depend on three and five parameters, respectively, which were determined by fitting to proton densities and energies. In particular, the epc17-2 functional33 has been parameterized to produce an effective balance between accurate proton densities and energies. The NEO-DFT method in conjunction with these functionals has been shown to produce accurate proton densities, proton affinities, and optimized geometries for the chemical systems studied.32,33,64,66 These epc functionals have also been shown to be transferable across a wide range of electronic functionals.64
The implementation of NEO-DFT analytical gradients enables geometry optimizations. The geometries of protonated water tetramers were optimized using the NEO-DFT method with the B3LYP electronic functional67,68 and the epc17-2 electron–proton correlation functional, treating all nine protons quantum mechanically. The relative stabilities of the protonated water tetramer isomers computed at the NEO-DFT level are consistent with the NEO-CCSD results, as well as the conventional CCSD(T) results including zero-point energy contributions shown in Fig. 2. The lack of NEO-CCSD analytical gradients prevented geometry optimizations at that level of theory, and the NEO-DFT geometry optimizations allowed further analysis of geometric effects for this system. In terms of geometry optimizations, the constrained NEO-DFT (cNEO-DFT) method69 produces the same stationary points as those of the NEO-DFT method but also allows the calculation of an extended NEO potential energy surface that depends on the expectation values of the quantum nuclear positions. This type of extended NEO potential energy surface may provide useful conceptual insights.
Thus, the NEO-DFT method is computationally efficient and reasonably accurate for computing ground state molecular properties. The implementation of NEO-DFT analytical Hessians allows the calculation of transition states, zero-point energy contributions, thermochemical properties, and MEPs. In contrast to the wave function methods, however, the systematic improvement of the NEO-DFT method is not straightforward. The main source of uncertainty is the electron–proton correlation functional. Although the epc17 and epc19 functionals perform well for many applications, they still have limitations in terms of predicting properties that are sensitive to the subtle, quantitative aspects of the proton densities. The development of more accurate and robust electron–proton correlation functionals is an important direction to explore. Possible strategies include the multicomponent extension of the SCAN functional,70 which was designed to satisfy exact constraints, or machine learning approaches to functional development.71
Excited states may be studied with the NEO-TDDFT method,36,37 which is analogous to its conventional electronic counterpart.25 The linear response of the NEO Kohn–Sham system to perturbative external nuclear and electronic fields leads to a matrix equation for the linear-response NEO-TDDFT method.36 The solution of this matrix equation provides the excitation energies and corresponding transition amplitudes for electronic, protonic, and mixed electron–proton vibronic excitations in a single calculation. Because the derivation relies on the adiabatic approximation, where the kernel is assumed to be independent of frequency, this method captures only single excitations, which can be described as linear combinations of products of electronic and protonic determinants with only one determinant singly excited in each term. In addition to the excitation energies, the transition densities, transition dipole moments, oscillator strengths, and intensities can be calculated to characterize each excitation.37
The low-lying electronic and protonic excitations are separable for predominantly electronically adiabatic systems but become mixed for higher excitations and for nonadiabatic systems. The electronic excitation energies computed with NEO-TDDFT for the lower electronic states are similar to those computed with conventional electronic TDDFT,36,37 but higher electronic states exhibit differences due to vibronic mixing.38 The key advantage of the NEO-TDDFT method is that it also describes proton vibrational excitations and, in some cases, mixed electron–proton vibronic excitations. The proton vibrational excitation energies for the fundamental vibrational modes (i.e., the bend and stretch modes of a molecule such as HCN) have been shown to be accurate, often within 20 cm−1 of nearly numerically exact grid-based reference values computed for the ground state potential energy surface generated with the same electronic structure method.37 This accuracy carries over to systems with multiple protons, such as C2H2, where the six fundamental proton vibrational modes are linear combinations of single excitations associated with each proton.72 In contrast to the NEO-EOM-CCSD method,35 the combination bands and overtones are not described accurately, presumably due to the limitations of the electron–proton correlation functionals, and double excitations are absent due to the underlying adiabatic approximation.
The implementation of NEO-TDDFT analytical gradients allows the optimization of excited state geometries and the calculation of adiabatic excitation energies.73 Within the NEO framework, the adiabatic excitation energy is defined to be the energy difference between geometries optimized on the ground and excited state NEO electron–proton vibronic potential energy surfaces and thus implicitly includes the zero-point energies and associated anharmonicities of the quantum protons. The NEO-DFT analytical Hessian and the NEO-TDDFT semi-numerical Hessian can be used to compute the zero-point energy contributions from the classical nuclei within the harmonic approximation for the ground and excited state NEO vibronic surfaces. A combination of these zero-point energy contributions and the adiabatic excitation energies produces the 0–0 adiabatic excitation energies, which can be compared directly to experimental data.73 The NEO-TDDFT analytical gradients also enable the simulation of adiabatic dynamics on excited state vibronic surfaces and provide the foundation for nonadiabatic surface hopping simulations.
The NEO-DFT(V) method44 was developed to compute molecular vibrational frequencies that can be compared to experimental spectra. Because the NEO potential energy surface depends only on the classical nuclei, the NEO Hessian produces modes that depend only on the classical nuclei. The NEO-DFT(V) method couples the NEO-DFT Hessian for the classical nuclei with NEO-TDDFT vibrational excitation energies and transition dipole moments for the quantum nuclei to generate the molecular vibrational modes composed of both quantum and classical nuclear motions. This method entails the calculation of an extended Hessian that depends on the classical nuclear coordinates and the expectation values of the quantum proton coordinates, where the diagonal block matrix of the Hessian associated with the quantum protons is obtained from NEO-TDDFT calculations. Diagonalization of this extended Hessian provides the coupled vibrational modes depending on both classical and quantum nuclear coordinates. A schematic illustration of this approach for HCN is given in Fig. 4. The NEO-DFT(V) method has been shown to incorporate the significant anharmonic effects associated with the quantum protons into the molecular vibrational frequencies, leading to good agreement with experimental data as well as calculations that include anharmonic effects perturbatively.44,72
The NEO-TDDFT method has been shown to provide electronic, protonic, and vibronic excitations in a computationally efficient manner. As discussed above in the context of the NEO-DFT method, the development of more reliable electron–proton correlation functionals could improve the quantitative accuracy of the combination bands and overtones as well as the fundamental vibrational modes. For some systems, double excitations can be computed by combining the NEO-ΔSCF and NEO-TDDFT methods,17,74 although complications with spin contamination, triplet instabilities, and unstable solutions may arise. The exploration of alternative approaches for computing double excitations within the NEO-TDDFT framework is an important future direction. In terms of the NEO-DFT(V) method for calculating molecular vibrational frequencies, the anharmonic effects associated with the quantum nuclei are included through the NEO-TDDFT calculations, but the modes dominated by classical nuclei are still calculated within the harmonic approximation. As discussed below, this issue is addressed by the NEO-Ehrenfest dynamics approach for computing vibrational frequencies.39
NONEQUILIBRIUM DYNAMICS BEYOND THE BORN–OPPENHEIMER APPROXIMATION
An alternative to the linear-response NEO-TDDFT approach is the real-time NEO-TDDFT (RT-NEO-TDDFT) approach.38 In this case, substitution of the nuclear–electronic wave function of the form given in Eq. (1) into the time-dependent Schrödinger equation leads to the following time-dependent equations:
These electronic and protonic equations are strongly coupled because the electronic and protonic Fock matrices each depend on both the electronic and protonic orbital coefficients. After some rearrangements, these two sets of equations are propagated numerically to obtain the orbital coefficients as a function of time. The electronic and protonic excitation energies are computed from the Fourier transform of the time-dependent electronic and protonic dipole moments. The resulting excitation energies are in agreement with the linear-response NEO-TDDFT results.38
An advantage of the RT-NEO-TDDFT method is that it provides the nonequilibrium dynamics of the electronic–protonic system without invoking the Born–Oppenheimer separation. For example, an electric field pulse can be applied to a molecular system to simulate ultrafast vibrational excitation by a laser. The application of this approach to the HCN molecule illustrates that a resonant driving field leads to energy absorption that produces time-dependent oscillations of the dipole moment. These types of calculations could assist in the interpretation of time-resolved infrared spectroscopy.38
The RT-NEO-TDDFT method has also been applied to excited state intramolecular proton transfer (ESIPT) in o-hydroxybenzaldehyde (oHBA).38 In this system, photoexcitation induces intramolecular proton transfer between two oxygen atoms. Photoexcitation to the excited electronic state, corresponding to an S0 → S1 excitation, is modeled by an electronic transition from the highest occupied molecular orbital to the lowest unoccupied molecular orbital. In the initial studies, the transferring hydrogen nucleus was represented by three sets of electronic and protonic basis functions spanning the initial and final positions of the hydrogen nucleus. When these RT-NEO-TDDFT calculations were performed at the equilibrium ground state geometry, the proton did not transfer because the classical nuclei were fixed and unable to reorganize. When these calculations were performed at a geometry optimized in the excited state with the proton constrained to be bound to the donor oxygen atom, the proton transferred to the acceptor oxygen atom. These calculations highlight the strong coupling between the transferring proton and the other nuclei as well as the importance of excited state structural relaxation in these types of photoinduced reactions.
To simulate the nonadiabatic dynamics of the classical nuclei within the NEO framework, the RT-NEO-TDDFT method can be combined with the Ehrenfest dynamics approach.39,40 In this case, the quantum subsystem composed of the electrons and quantum protons is propagated according to Eq. (6), while the classical nuclei are propagated on the average surface determined by the time-dependent electronic–protonic Kohn–Sham wave function
where MI and are the masses and coordinates, respectively, of the classical nuclei. The use of fixed electronic and protonic basis function centers for the quantum protons requires a large number of basis function centers to describe a proton transfer reaction and prior knowledge of its trajectory. To address this issue, a semiclassical traveling proton basis function approach has been developed.39 In this approach, the basis function centers associated with the quantum protons follow the same equations of motion as the classical nuclei. Additional terms related to the time derivatives of the proton basis function centers must be added to the protonic equations of motion given in Eq. (6). The semiclassical traveling proton basis function method will converge to the exact result within the RT-NEO-TDDFT Ehrenfest treatment of the dynamics as the basis set approaches completeness, although it is approximate for incomplete basis sets. The molecular vibrational frequencies computed with the RT-NEO-TDDFT Ehrenfest method are consistent with those computed with the NEO-DFT(V) method for the chemical systems studied,39 although the NEO Ehrenfest approach includes anharmonic effects associated with all nuclei.
The application of the RT-NEO-TDDFT Ehrenfest method to ESIPT in oHBA enables the nonequilibrium dynamical simulation of the electrons and all nuclei beyond the Born–Oppenheimer approximation.40 The nonadiabatic effects between the electrons and quantum protons are included via the RT-NEO-TDDFT method, and the nonadiabatic effects between the classical nuclei and the electron–proton quantum subsystem are described with Ehrenfest dynamics. In this case, photoexcitation at the equilibrium ground state geometry leads to proton transfer because the classical nuclei are allowed to reorganize in the excited electronic state. A comparison of Ehrenfest dynamics with a quantum or classical treatment of the transferring hydrogen nucleus shows that the quantization of the hydrogen nucleus within the NEO framework leads to faster proton transfer, as shown in Fig. 5. The delocalization of the quantized proton allows the transfer to occur at a longer distance between the donor and acceptor oxygen atoms, thereby requiring less movement of the oxygen atoms prior to proton transfer. Moreover, the kinetic isotope effect for this proton transfer reaction is slightly larger with a quantum treatment of the transferring hydrogen or deuterium.
The simulation of nonadiabatic dynamics within the NEO framework opens up a wide range of opportunities. In addition to Ehrenfest dynamics, the NEO method can also be combined with other nonadiabatic dynamics methods, such as trajectory surface hopping.6 In this case, the fewest switches surface hopping algorithm could be applied to trajectories associated with the classical nuclei moving on the adiabatic vibronic surfaces computed with linear-response NEO-TDDFT. In principle, the NEO approach could also be combined with the multiple spawning method.7,8 Various technical challenges, such as the calculation of the nonadiabatic coupling vector between NEO-TDDFT vibronic states, would need to be addressed for these types of simulations. In addition, more robust strategies for moving the proton basis function centers could improve the accuracy and energy conservation of the nondiabatic dynamics, and spin-flip methods may be useful. The description of hydrogen tunneling within these nonadiabatic dynamics simulations also necessitates further developments, as discussed below.
The description of hydrogen tunneling requires the calculation of bilobal, delocalized hydrogen vibrational wave functions. The variationally optimized proton wave function or density computed at the NEO-HF or NEO-DFT level for a proton moving in a symmetric double-well potential is localized in one well even when the transferring hydrogen is represented by two basis function centers (i.e., a set of electronic and protonic basis functions centered at each minimum).75 Thus, a multireference approach is required to produce bilobal wave functions in this situation. Although the NEO complete active space SCF (NEO-CASSCF) approach was developed nearly two decades ago,16 extremely large active spaces are required to describe such bilobal wave functions because dynamical as well as non-dynamical correlation is important.75,76 Similarly, NEO multireference CI (MRCI) calculations require a very large number of determinants in the expansion.
As a computationally practical alternative, the NEO multistate DFT (NEO-MSDFT) method41 includes both dynamical and non-dynamical correlation and has been shown to be capable of describing hydrogen tunneling systems. In the NEO-MSDFT approach, two localized electronic–protonic wave functions are obtained with the NEO-DFT method using two basis function centers to represent the transferring hydrogen. A nonorthogonal configuration interaction approach is used to mix these two localized wave functions to produce bilobal, delocalized wave functions, as shown in Fig. 6. This method only requires the diagonalization of a 2 × 2 Hamiltonian matrix in the basis of the two localized NEO-DFT states, ΨI and ΨII,
where the diagonal matrix elements are the NEO-DFT energies of the localized states.
The off-diagonal element can be approximated by a physically motivated form inspired by the conventional electronic MSDFT method.77,78 In the NEO framework, this element is expressed as
Here, the first term is the off-diagonal matrix element computed at the NEO-HF level with the NEO Kohn–Sham determinants, SI,II is the overlap between the localized Kohn–Sham wave functions, and the terms in parentheses correspond to the NEO-DFT and NEO-HF energies of the localized states. A correction function is applied to SI,II in order to account for limitations of the electron–proton correlation functionals and the approximate form of the off-diagonal matrix element.41 This correction function requires two parameters that were determined for simple model systems and fixed thereafter. The NEO-MSDFT method has been shown to produce quantitatively accurate hydrogen tunneling splittings with errors less than ∼6 cm−1 for systems such as malonaldehyde, acetoacetaldehyde, and the benzyl-toluene radical complex.41
Although the NEO-MSDFT method is computationally efficient and quantitatively accurate in its current form, systematic improvement is not straightforward. More rigorous multireference wave function approaches, such as NEO-CASSCF or MRCI methods,16,76,79,80 are systematically improvable but are significantly more computationally expensive and may not be easily applied to large molecular systems. The importance of dynamical correlation most likely necessitates a method such as the analog of the conventional electronic second-order perturbative (CASPT2) method81 to obtain quantitatively accurate hydrogen tunneling splittings. In conjunction with analytical gradients, multireference methods such as NEO-MSDFT and NEO-CASSCF can also be used to generate MEPs that include hydrogen tunneling effects. In addition, such multireference methods could be combined with the NEO Ehrenfest or surface hopping algorithms for nonadiabatic dynamics simulations of hydrogen tunneling processes. As discussed in the context of other methods, embedding approaches may be helpful for applications of these multireference methods to larger systems.
This Perspective outlines the foundational NEO methods, including NEO-HF, NEO-CI, NEO-CCSD, NEO-OOMP2, NEO-EOM-CCSD, NEO-DFT, NEO-TDDFT, RT-NEO-TDDFT, NEO-Ehrenfest, NEO-MSDFT, and NEO-CASSCF, as depicted in Fig. 7. In each section, the limitations of the methods and the opportunities for future developments are discussed. Most of these methods are implemented in the Q-Chem 5.4 software package82 or a development branch. The real-time NEO methods are implemented in a development branch of the Chronus Quantum software package.83 Some of the NEO methods are also available in other quantum chemistry software packages.
These NEO methods will serve as the building blocks for constructing a framework of strategies aimed at diverse applications. For example, many of these NEO methods could be extended to periodic systems, allowing the study of nuclear quantum effects in solids and at interfaces. In principle, all nuclei could be treated quantum mechanically within the NEO framework,84 as long as the translations and rotations are treated appropriately. The advantages of including all nuclear quantum effects at the same level may warrant the added complexity. However, the treatment of only specified nuclei within the NEO framework is appealing for its simplicity, accessibility, and computational efficiency. Solvation effects can be included using either implicit or explicit solvation models. A particularly promising strategy relies on embedding85 or hybrid approaches that treat a relatively small portion of the system at a correlated NEO level and treat the environment at a lower level of quantum mechanics or with a molecular mechanical force field. Exploring these directions will open up possibilities for simulating a wide range of nonadiabatic processes in complex environments.
The NEO work discussed in this Perspective would not have been possible without key contributions from the following group members: Simon Webb, Chet Swalina, Michael Pak, Ari Chakraborty, Tanner Culpitt, Kurt Brorsen, Yang Yang, Fabijan Pavošević, Patrick Schneider, Coraline (Zhen) Tao, Ben Rousseau, and Qi Yu. The real-time NEO work was performed in collaboration with Xiaosong Li’s group, with key contributions from Luning Zhao and Andrew Wildman. I am also grateful for many useful discussions with John Tully. This work was supported mainly by the National Science Foundation under Grant No. CHE-1954348.
Data sharing is not applicable to this article as no new data were created or analyzed in this study.