Nuclear–electronic orbital Ehrenfest dynamics

Across chemistry and physics, nuclear quantum effects play a vital role in many processes of great technological importance. Spectroscopic techniques probing quantized vibrational states have become essential tools to reveal the structure of molecules, especially when the molecules are undergoing rapid changes such as bond breaking and formation.^1,2 Hydrogen tunneling³ is another prominent nuclear quantum effect that is often important in proton transfer and proton-coupled electron transfer (PCET).^4–6 PCET is a key step in a wide range of chemical and biological reactions, including photosynthesis, respiration, DNA damage, and reactions in solar cells. In addition to hydrogen tunneling, carbon⁷ and oxygen⁸ tunneling have also been observed in experiments.

Due to the importance of nuclear quantum effects in spectroscopy and dynamics, developing accurate first principles theoretical approaches to model these effects is highly desirable. In order to accurately simulate processes such as PCET, which are often nonadiabatic, two effects need to be taken into account. First, the theoretical approach must incorporate the nuclear quantum effects, such as zero-point energy, tunneling, and quantized vibrational states, especially for the transferring protons. Second, the theoretical approach must capture the interplay between the nuclear wave function and the changes in the electronic configuration. In processes such as PCET, the widely used Born–Oppenheimer approximation breaks down due to the strong nonadiabatic effects between the electronic and nuclear degrees of freedom.⁶ 

Various approaches⁹ have been developed to simulate nonadiabatic dynamics, including the multiconfigurational time-dependent Hartree (MCTDH),^10,11 variational multiconfigurational Gaussian (vMCG),¹² full multiple spawning (FMS),^13,14 and exact factorization methods.^15,16 Given a sufficiently large number of grid points or basis functions along with accurate calculations of potential energy surfaces and non-adiabatic couplings, these approaches could produce highly accurate solutions to the time-dependent Schrödinger equation. However, most of these approaches become computationally expensive for relatively large chemical systems.

Another promising method that incorporates both nonadiabatic coupling and nuclear quantum effects in a computationally efficient manner is the nuclear–electronic orbital (NEO)^17–19 approach. The NEO approach is a direct generalization of standard electronic structure methods to incorporate nuclear quantum effects and nonadiabatic coupling. In this approach, electrons and certain key nuclei, typically protons, are treated quantum mechanically on equal footing, avoiding the Born–Oppenheimer separation between the specified nuclei and the electrons. Various electronic structure methods have been generalized to the NEO framework, including Hartree–Fock (HF),¹⁷ density functional theory (DFT),^20–23 second order perturbation theory (PT2),²⁴ coupled-cluster with singles and doubles (CCSD),^25–27 and multi-reference methods.^17,18 By properly treating electron–proton correlation, either by a correlated wave function approach or a density functional approach, accurate proton density distributions and proton affinities can be obtained. Excited state properties, such as vibrational excitation energies and oscillator strengths, can be modeled by linear response²⁸ and equation-of-motion (EOM)^25,29 formalisms. The NEO-TDDFT approach yields accurate predictions of vibrational frequencies for fundamental excitations.^28,30 A more recent study using the NEO-EOM-CCSD^26,29 approach also produces accurate excitation energies of overtones and combination bands.

Recently, we have extended NEO-HF and NEO-DFT to the real-time (RT)^31,32 domain by employing the time-dependent multicomponent Runge–Gross theory³³ within the adiabatic approximation.³⁴ The unique strength of the RT-NEO framework is that it enables simulations of nonequilibrium electronic–nuclear dynamics, as exemplified by excited state intramolecular proton transfer (ESIPT).³¹ The RT-NEO methods treat the quantum subsystem, including the electrons and selected protons, with the time-dependent variational principle, while the rest of the molecular system remains as frozen classical particles. To address the dynamical interplay between the quantum and classical subsystems, a coupled equation of motion that can take into account the mutual dynamical feedback between the quantum and classical degrees of freedom is needed. For example, in dynamical processes such as ESIPT, the frozen nuclei formulation ignores classical backbone displacements, which have been shown to be critical in modulating the characteristics and rates of proton transfer.³⁵ In addition, the vibronically nonadiabatic effects related to the dynamics of the classical subsystem and its coupling to the quantum electron–proton subsystem have been shown to be important in processes such as PCET.^6,36

In this work, we formulate the theoretical and mathematical underpinnings of the NEO-Ehrenfest approach in which the quantum subsystem modeled with the RT-NEO method is coupled to the classical subsystem via the Ehrenfest theorem. This approach takes into account the vibronic nonadiabaticity and the dynamical interplay between the quantum electron–proton subsystem and the classical subsystem within the NEO framework. In contrast to the conventional Ehrenfest approach,^37–39 where only electrons are considered in the quantum subsystem, the NEO-Ehrenfest method propagates the wave functions of electrons and quantum protons with the time-dependent Schrödinger equation, while the classical nuclei follow the Newtonian equations of motion with a mean-field potential, averaging over the ground and excited vibronic states. In this work, we also introduce a traveling proton basis approach with the aid of the time-dependent Schrödinger equation and semiclassical approximations.

This paper is structured as follows: We start with a brief overview of the Ehrenfest algorithm. Then, we present detailed expressions for the nuclear gradients of the energy for both NEO-TDHF and NEO-TDDFT. After presenting the gradient expressions, we summarize the algorithm for NEO-Ehrenfest dynamics. We then give an equation of motion for the proton basis function centers and discuss the corresponding modified Ehrenfest algorithm. The results are then presented for the fundamental vibrational frequencies of several small molecules containing one quantum proton. We conclude our discussions by pointing out the merits and shortcomings of the current method along with possible future directions.

The starting point of the derivation of Ehrenfest dynamics is a product ansatz⁴⁰ for the total wave function $Ψre,rp,R,t$⁠,

Here, Φ(r^e, r^p, t) is the wave function for the electrons with coordinates r^e and the quantum protons with coordinates r^p, and Θ(R, t) is the wave function for the other nuclei with coordinates R. E_q(t) is a phase factor defined by

where $Ĥq$ is the Hamiltonian for the electrons and quantum protons.

We now separate the equations of motion for the light (electrons and quantum protons) and heavy (other nuclei) particles and take the classical limit of the latter. We then arrive at the equations of motion for Ehrenfest dynamics,

where the dynamics of the classical nuclei follow the Newtonian equation of motion [Eq. (3)] with a mean-field potential generated by the electronic and protonic wave functions.

For NEO-TDHF and NEO-TDDFT, the quantum subsystem [Eq. (4)] can be formulated as two coupled equations of motion for the electronic and protonic orbitals,

where C^e(t) and C^p(t) are the time-dependent orbital coefficients for the electrons and protons, respectively, and F is the Fock/Kohn–Sham matrix in the orthonormal atomic orbital basis. Note that we will use primed notations (e.g., F′) for quantities in the nonorthogonal atomic orbital basis. In this work, the orthonormal basis is obtained with the Löwdin orthogonalization scheme.

The Ehrenfest dynamics procedure using RT-NEO-TDHF or RT-NEO-TDDFT requires computing the energy gradient with respect to the classical nuclear coordinates. In RT-NEO-TDHF and RT-NEO-TDDFT, the time-dependent wave function is in a superposition state, which leads to a gradient expression that differs from the expressions for equilibrium states, where the quantum degrees of freedom are in an eigenstate.^37–39 

We first write down the total energy expression for NEO-HF and NEO-DFT as

in which P^e′ and P^p′ are the total one-body density matrices for electrons and protons, respectively, evaluated in the non-orthogonal atomic orbital basis. h^e(p)′ is the one-body potential matrix for electrons (protons). J^ee(pp)′ and K^ee(pp)′ are the Coulomb and exchange matrices, respectively, for electrons (protons), evaluated using P^e′ (P^p′). J^ep′ accounts for the electron–proton Coulombic interaction. $Excee$ is the exchange-correlation energy for the electrons, and E_epc is the electron–proton correlation energy. V_NN is the repulsion energy between classical nuclei. ξ is a parameter that converts the energy between NEO-HF (ξ = 1) and pure NEO-DFT (ξ = 0). Note that this expression assumes a high-spin protonic configuration and only considers the proton–proton exchange interaction because proton–proton correlation has been shown to be negligible in molecular systems with localized proton orbitals.¹⁸ The high-spin protonic configuration results in a factor of 1/2 instead of 1/4 in the proton–proton exchange energy. In general, the proton–proton exchange term could be multiplied by the parameter ξ, and a proton–proton exchange-correlation energy could be included, analogous to the treatment of the electron–electron exchange-correlation energy.

The one-body potential matrices are

in which A is the index for classical nuclei and Z_A and R_A are their corresponding charges and positions. Here, m_p is the proton mass.

Assuming a closed-shell electronic configuration and a high-spin protonic configuration, the non-orthogonalized density matrices are computed as

in which C^e′ and C^p′ are the electronic and protonic orbitals in the non-orthogonalized basis. The orthogonalized density matrices are defined as

in which S^e(p) is the orbital overlap matrix for the electrons (protons).

The Coulomb and exchange matrices for electrons and protons are computed as

The electron–proton Coulomb matrix is computed as

The two-body integrals are defined in terms of the electronic and protonic basis functions as

The expression for a hybrid exchange-correlation functional within the generalized gradient approximation (GGA) for electrons is

where c is the fraction of the HF exchange.

The E_epc term can be treated with any available electron–proton correlation functional. In this work, we use the epc17-2 electron–proton correlation functional,^21,22

Here, the electron (proton) density is

The energy gradient with respect to nuclear coordinates can be computed as

The first term contains the gradients of the one- and two-body integrals and the nuclear repulsion energy, whereas the last two terms arise via the chain rule from the dependence of the atom-centered basis on the nuclear coordinates.

The first term depends on the nuclear coordinates explicitly,

where the gradients of J and K are

The gradients of $Excee$ and E_epc are

in which f and ϵ are the integrands of the electron–electron exchange-correlation functional and electron–proton correlation functional.

The last two terms in Eq. (16) can be evaluated as

where we have made the approximation that the density matrix in the orthonormal basis P does not depend on the nuclear coordinates. This approximation is common in conventional Ehfenfest algorithms, and for a finite atomic orbital basis, it has been shown to be valid in terms of energy conservation.^37–39,41 Note that due to the stationary condition of the energy on the density, coupled-perturbed equations do not need to be solved to obtain the first-order nuclear gradient.⁴² The gradient of S^−1/2 is

Substituting Eq. (21) into Eq. (20), we obtain

in which the gradient of S^1/2 can be computed as

where s_i and σ_i are the ith eigenvector and eigenvalue, respectively, of the overlap matrix S.

Combining Eqs. (17), (22), and (23), we arrive at the final working expression for the NEO-Ehrenfest energy gradients,

The NEO-Ehrenfest dynamics procedure requires the integration of the equations of motion for the classical nuclei and the RT-NEO-TDHF/TDDFT equations for the quantum subsystem. The velocity Verlet algorithm is used to propagate the classical degrees of freedom,

in which x_I, p_I, g_I, and M_I are the position, momentum, gradient, and mass of the Ith classical nucleus, respectively, and Δt_N is the integration time step.

Because electrons are much lighter than nuclei, a different smaller time step should be used to integrate the electronic degrees of freedom, and this time step is usually two or three orders of magnitude smaller than Δt_N. In order to account for the responses of the quantum protons and electrons to each other, we use the same time step to integrate the protonic orbitals. The time step used to integrate all quantum degrees of freedom is denoted as Δt_q.

In addition to Δt_N and Δt_q, we follow the triple-split approach developed by Li and Tully and introduce a third intermediate time step Δt_Nq.^37–39 The most computationally expensive step in the Ehrenfest dynamics is the evaluation of the energy gradient. The intermediate step Δt_Nq is introduced to account for the changes in classical nuclear coordinates in the Hamiltonian and the orthonormalization condition without recomputing the energy gradient. We use Δt_Nq with a value in between Δt_N and Δt_q as the time to re-evaluate the integrals and update the Fock matrix accordingly.

So far, we have assumed that the proton basis function centers are fixed parameters. However, since Gaussian basis functions are used to model proton wave functions and, typically, these basis functions are centered at a single position, they can only span a limited amount of coordinate space in the vicinity of the proton. Therefore, using fixed basis function centers is reasonable only when protons do not move too far from their initial positions, e.g., small oscillations around equilibrium driven by a weak external field. For extended dynamics in which the proton travels significantly through coordinate space, the use of fixed basis function centers will either fail to span the spatial extent of the dynamics or require more basis functions than that are computationally feasible. For example, in proton transfer processes, the protonic basis set must be able to span the space composed of all potential proton transfer pathways in order to correctly describe the protonic quantum dynamics. Using a limited stationary protonic basis set could result in qualitatively incorrect dynamics and even prevent proton transfer altogether. On the other hand, introducing a large number of protonic basis function centers to cover all potential pathways will make the computational approach prohibitively expensive. Thus, it is necessary to develop an algorithm that can dynamically adjust the proton basis function centers as the protons travel. In this section, we present an approach to achieve this goal based on the time-dependent Schrödinger equation and semiclassical approximations.

The equations of motion for the proton orbitals are governed by the time-dependent Schrödinger equation,

where $CQip′$ are the orbital coefficients and $RQ(t)$ is the time-dependent center for the Qth proton basis function.

Multiplying Eq. (26) to the left by $ϕQp$ leads to the equations of motion for the orbital coefficients,

where

and F^p′ is the protonic Fock matrix defined as

Equation (29) can be expressed in matrix form as

Then, Eq. (31) can be transformed to the orthonormal basis using

so that

Multiplying Eq. (35) by C^p† and using $PPQp=∑iCQip*CPip$⁠, we obtain

and its complex conjugate,

Subtracting Eq. (37) from Eq. (36), we arrive at the equation of motion for the time-dependent proton density matrix,

The dependence of the time evolution of the proton density matrix on the time-varying basis functions is embedded in the τ term, which is evaluated as

From Eqs. (38) and (39), we can see that the atomic basis function centers enter the equation of motion for the proton density with an associated momentum via the $ṘPα$ term, which can be evaluated using various approaches.

If the basis function centers are frozen or time-independent, $Ṙ=0$ and τ = 0. As a result, Eq. (38) becomes Eq. (5), the equation of motion for proton basis functions with fixed centers. In this work, we introduce a semiclassical approximation as the equation of motion for the basis function centers,

in which $m̃p$ is the fictitious proton mass and the energy gradient is taken with respect to the proton basis function center. Note that as the proton basis function center moves, the associated electronic basis functions also travel along with the protonic basis functions. The semiclassical approximation introduces an arbitrary parameter, as it moves the proton basis function center with a fictitious mass. In this work, the proton mass is used. The choice of the fictitious mass is based on computations of vibrational frequencies of selected small molecules using NEO-Ehrenfest, showing that using the proton mass in Eq. (40) produces overall the best results.

The Algorithm presents the computational strategy for NEO-Ehrenfest dynamics with and without a traveling proton basis. For dynamics with a traveling proton basis, the quantum degrees of freedom are propagated using Eq. (38) instead of Eq. (5). In addition, one needs to compute the gradient and update the position and momentum not only for the classical nuclei but also for the proton basis function centers. Note that Eq. (38) is only applied to the protonic degrees of freedom to account for the incompleteness of the protonic basis. In this work, we do not use the traveling basis equation of motion for the electronic density because there is a sufficient number of electronic basis functions to describe electronic fluctuations.

The NEO-Ehrenfest algorithms are implemented in the Chronus Quantum⁴³ open source package. NEO-HF/DFT(V)⁴⁴ and second-order vibrational perturbation theory (VPT2)⁴⁵ calculations were performed with a developer’s version of QChem 5.3.⁴⁶ For all systems studied, the cc-pVDZ electronic basis set⁴⁷ was used for non-hydrogen atoms, and the cc-pV5Z electronic basis set and PB4-F2 protonic basis set⁴⁸ were used for hydrogen. The cc-pVDZ electronic basis set was used for all atoms in the VPT2 calculations. All NEO-Ehrenfest dynamics simulations were propagated for 100 fs, with Δt_N = 0.1 fs, Δt_Nq = 0.01 fs, and Δt_q = 0.001 fs. The energy conservation error for computing these vibrational spectra is on the order of 10⁻⁹ Hartree. The molecular geometries for HCN, HNC, and FHF⁻ were obtained from the corresponding NEO-HF and NEO-DFT ground state geometry optimizations (see the supplementary material for Cartesian coordinates). All the NEO-DFT calculations used the B3LYP^49–51 electronic functional and the epc17-2^21,22 electron–proton correlation functional. The initial nuclear and electronic basis function centers for the quantum proton were positioned at the optimized hydrogen position. Linear dependencies are removed from the basis set during the orthogonalization process to avoid singularity in the overlap matrix.

Although vibrational spectra can be easily obtained from response theories,²⁸ reproducing these spectra serves as a useful benchmark study to validate the NEO-Ehrenfest dynamics introduced herein. We apply the NEO-Ehrenfest algorithm to compute molecular vibrational frequencies of HCN, HNC, and FHF⁻. The vibrational frequencies are obtained by first converging NEO-HF/DFT at the optimized geometry and then performing NEO-Ehrenfest simulations by initially perturbing the positions of N, C, and F by 10⁻⁵ Bohr in the x, y, or z direction. We then compute the time-dependent proton dipole moments as

where α ∈ {x, y, z} and $dαp′$ is the dipole integral matrix with elements $dα,PQp′=ϕPp|rα|ϕQp$⁠. The proton vibrational frequencies are obtained by finding the peaks of the absorption cross section, computed as $σ(ω)∝ω∑α=x,y,zIm[D̃α(ω)]$⁠, where $D̃α(ω)$ is the Fourier transform of the time-dependent proton dipole moment. The vibrationalfrequencies of modes that depend only on heavy atoms, including the C–N, N–C, and F–F stretch modes, are computed by Fourier transform of the time-dependent variation of the corresponding bond length. Note that the frequencies corresponding to predominantly heavy atom and mixed modes can also be obtained from the Fourier transform of the time-dependent proton dipole moment for modes with non-zero oscillator strength.