The recently developed real-time nuclear–electronic orbital (RT-NEO) approach provides an elegant framework for treating electrons and selected nuclei, typically protons, quantum mechanically in nonequilibrium dynamical processes. However, the RT-NEO approach neglects the motion of the other nuclei, preventing a complete description of the coupled nuclear–electronic dynamics and spectroscopy. In this work, the dynamical interactions between the other nuclei and the electron–proton subsystem are described with the mixed quantum–classical Ehrenfest dynamics method. The NEO-Ehrenfest approach propagates the electrons and quantum protons in a time-dependent variational framework, while the remaining nuclei move classically on the corresponding average electron–proton vibronic surface. This approach includes the non-Born–Oppenheimer effects between the electrons and the quantum protons with RT-NEO and between the classical nuclei and the electron–proton subsystem with Ehrenfest dynamics. Spectral features for vibrational modes involving both quantum and classical nuclei are resolved from the time-dependent dipole moments. This work shows that the NEO-Ehrenfest method is a powerful tool to study dynamical processes with coupled electronic and nuclear degrees of freedom.

## I. INTRODUCTION

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^{3} 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^{7} and oxygen^{8} 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.^{6}

Various approaches^{9} have been developed to simulate nonadiabatic dynamics, including the multiconfigurational time-dependent Hartree (MCTDH),^{10,11} variational multiconfigurational Gaussian (vMCG),^{12} 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),^{17} density functional theory (DFT),^{20–23} second order perturbation theory (PT2),^{24} 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^{28} 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^{33} within the adiabatic approximation.^{34} 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).^{31} 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.^{35} 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.

## II. THEORY

### A. Ehrenfest dynamics

The starting point of the derivation of Ehrenfest dynamics is a product ansatz^{40} for the total wave function $\Psi 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 $\u0124q$ 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.

### B. Ehrenfest energy gradient for nonequilibrium NEO-TDHF and NEO-TDDFT

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.^{18} 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.^{42} The gradient of **S**^{−1/2} is

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

where **s**_{i} and *σ*_{i} are the *i*th 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,

### C. NEO-Ehrenfest integration algorithm

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 *I*th 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.

### D. Traveling proton basis

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\u2032$ are the orbital coefficients and $RQ(t)$ is the time-dependent center for the *Q*th proton basis function.

Multiplying Eq. (26) to the left by $\varphi 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=\u2211iCQip*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 $R\u0307P\alpha $ term, which can be evaluated using various approaches.

If the basis function centers are frozen or time-independent, $R\u0307=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\u0303p$ 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.

Δt_{Nq} = Δt_{N}/n, Δt_{q} = Δt_{Nq}/m |

for t = 0, Δt_{N}, 2Δt_{N}, … do |

Compute g(t) at x(t) [Eq. (24) for each classical nucleus |

Compute $p(t+12\Delta tN)$ and x(t + Δt_{N}) [Eq. (25)] for each classical nucleus |

if traveling proton basis then |

Compute g(t) at $R(t)$ [Eq. (24) for proton basis center |

Compute $p(t+12\Delta tN)$ and $R(t+\Delta tN)$ [Eq. (25)] for proton basis center |

end if |

for j = 1, n do |

t′ = t + (j − 1)Δt_{Nq} |

$xt\u2032\u2009+\u200912\Delta tNq\u2009=\u2009x(t)\u2009+\u2009(j\u22121)\Delta tNq\u2009+\u200912\Delta tNqpt+12\Delta tN/M$ |

Compute h^{e}′, h^{p}′, (pq|rs), (pq|RS), S^{e}, and $[Se]1/2$ based on $x(t\u2032+12\Delta tNq)$ |

for i = 1, m do |

τ = t′ + (i − 1)Δt_{q} |

$Pe(p)\u2032(\tau )=[Se(p)]\u22121/2Pe(p)(\tau )[Se(p)]\u22121/2$ |

Form F^{e(p)}′(τ) using P^{e(p)}′(τ) and the integrals evaluated at $x(t\u2032+12\Delta tNq)$ |

Propagate P^{e}(τ) →P^{e}(τ + Δt_{q}) [Eq. (5)] |

if traveling proton basis then |

Propagate P^{p}(τ) →P^{p}(τ + Δt_{q}) based on Eq. (38) |

else |

Propagate P^{p}(τ) →P^{p}(τ + Δt_{q}) based on Eq. (5) |

end if |

end for |

end for |

end for |

Δt_{Nq} = Δt_{N}/n, Δt_{q} = Δt_{Nq}/m |

for t = 0, Δt_{N}, 2Δt_{N}, … do |

Compute g(t) at x(t) [Eq. (24) for each classical nucleus |

Compute $p(t+12\Delta tN)$ and x(t + Δt_{N}) [Eq. (25)] for each classical nucleus |

if traveling proton basis then |

Compute g(t) at $R(t)$ [Eq. (24) for proton basis center |

Compute $p(t+12\Delta tN)$ and $R(t+\Delta tN)$ [Eq. (25)] for proton basis center |

end if |

for j = 1, n do |

t′ = t + (j − 1)Δt_{Nq} |

$xt\u2032\u2009+\u200912\Delta tNq\u2009=\u2009x(t)\u2009+\u2009(j\u22121)\Delta tNq\u2009+\u200912\Delta tNqpt+12\Delta tN/M$ |

Compute h^{e}′, h^{p}′, (pq|rs), (pq|RS), S^{e}, and $[Se]1/2$ based on $x(t\u2032+12\Delta tNq)$ |

for i = 1, m do |

τ = t′ + (i − 1)Δt_{q} |

$Pe(p)\u2032(\tau )=[Se(p)]\u22121/2Pe(p)(\tau )[Se(p)]\u22121/2$ |

Form F^{e(p)}′(τ) using P^{e(p)}′(τ) and the integrals evaluated at $x(t\u2032+12\Delta tNq)$ |

Propagate P^{e}(τ) →P^{e}(τ + Δt_{q}) [Eq. (5)] |

if traveling proton basis then |

Propagate P^{p}(τ) →P^{p}(τ + Δt_{q}) based on Eq. (38) |

else |

Propagate P^{p}(τ) →P^{p}(τ + Δt_{q}) based on Eq. (5) |

end if |

end for |

end for |

end for |

## III. RESULTS AND DISCUSSION

The NEO-Ehrenfest algorithms are implemented in the Chronus Quantum^{43} open source package. NEO-HF/DFT(V)^{44} and second-order vibrational perturbation theory (VPT2)^{45} calculations were performed with a developer’s version of QChem 5.3.^{46} For all systems studied, the cc-pVDZ electronic basis set^{47} was used for non-hydrogen atoms, and the cc-pV5Z electronic basis set and PB4-F2 protonic basis set^{48} 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^{−9} 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,^{28} 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^{−5} Bohr in the *x*, *y*, or *z* direction. We then compute the time-dependent proton dipole moments as

where *α* ∈ {*x*, *y*, *z*} and $d\alpha p\u2032$ is the dipole integral matrix with elements $d\alpha ,PQp\u2032=\u2009\varphi Pp|r\alpha |\varphi Qp\u2009$. The proton vibrational frequencies are obtained by finding the peaks of the absorption cross section, computed as $\sigma (\omega )\u221d\omega \u2211\alpha =x,y,zIm[D\u0303\alpha (\omega )]$, where $D\u0303\alpha (\omega )$ 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.

In Tables I and II, we compare the predicted vibrational frequencies using NEO-TDHF/TDDFT-Ehrenfest, NEO-HF/DFT(V), and VPT2. The NEO-HF/DFT(V) approach obtains vibrational frequencies by diagonalizing an extended NEO Hessian matrix that includes second derivatives of the NEO energy with respect to the expectation values of the quantum nuclei as well as the classical nuclei.^{44,52} This approach incorporates the coupling between the classical and quantum nuclei and the anharmonic effects associated with the quantum nuclei, although anharmonicities associated with the classical nuclei are not taken into account. The VPT2 approach is expected to be a more accurate approach for predicting vibrational frequencies because it accounts for the anharmonic effects associated with all nuclei via perturbation theory, although it does not include nuclear wave function delocalization and zero-point energy in the geometry optimizations. In the NEO-Ehrenfest approach, the anharmonicity of all vibrations can be directly resolved through molecular dynamics, but its accuracy and effectiveness depend on the approximate functionals and basis sets.

. | vibrational mode . | FPB . | TPB . | NEO-HF(V)^{a}
. | NEO-HF(V)-aug^{b}
. | VPT2 . |
---|---|---|---|---|---|---|

HCN | CH stretch | 3620 | 3483 | 3618 | 3456 | 3514 |

CH bend | 1294 | 1131 | 1473 | 883 | 846 | |

CN stretch | 2446 | 2405 | 2417 | 2407 | 2397 | |

HNC | NH stretch | 4009 | 3867 | 3994 | 3841 | 3879 |

NH bend | 1214 | 1057 | 1257 | 594 | 485 | |

NC stretch | 2304 | 2296 | 2293 | 2290 | 2262 | |

FHF^{−} | FH stretch | 1741 | 1644 | 1728 | 1518 | 1636 |

FH bend | 1622 | 1482 | 1617 | 1417 | 1412 | |

FF stretch | 675 | 689 | 686 | 686 | 548 |

. | vibrational mode . | FPB . | TPB . | NEO-HF(V)^{a}
. | NEO-HF(V)-aug^{b}
. | VPT2 . |
---|---|---|---|---|---|---|

HCN | CH stretch | 3620 | 3483 | 3618 | 3456 | 3514 |

CH bend | 1294 | 1131 | 1473 | 883 | 846 | |

CN stretch | 2446 | 2405 | 2417 | 2407 | 2397 | |

HNC | NH stretch | 4009 | 3867 | 3994 | 3841 | 3879 |

NH bend | 1214 | 1057 | 1257 | 594 | 485 | |

NC stretch | 2304 | 2296 | 2293 | 2290 | 2262 | |

FHF^{−} | FH stretch | 1741 | 1644 | 1728 | 1518 | 1636 |

FH bend | 1622 | 1482 | 1617 | 1417 | 1412 | |

FF stretch | 675 | 689 | 686 | 686 | 548 |

^{a}

The NEO-HF(V) calculations were performed with the same basis sets as the NEO-Ehrenfest dynamics with the proton basis function centers placed at the variationally optimized position.

^{b}

For the NEO-HF(V)-aug calculations, the NEO-TDHF calculations used to obtain the hydrogen vibrational excitations in the extended Hessian were performed using the Cartesian cc-pVTZ electronic basis set for the heavy nuclei and the Cartesian cc-pV6Z^{*} (no h function) electronic basis set for the quantum proton, as well as the Cartesian 8s8p8d8f protonic basis set, and the proton basis function centers were placed at the expectation value determined by a full geometry optimization with these basis sets, as prescribed in previous studies.^{44,52}

. | vibrational mode . | FPB . | TPB . | NEO-DFT(V)^{a}
. | NEO-DFT(V)-aug^{b}
. | VPT2 . |
---|---|---|---|---|---|---|

HCN | CH stretch | 3844 | 3521 | 3910 | 3297 | 3319 |

CH bend | 1135 | 866 | 1297 | 797 | 744 | |

CN stretch | 2278 | 2222 | 2212 | 2191 | 2173 | |

HNC | NH stretch | 3970 | 3791 | 3991 | 3635 | 3604 |

NH bend | 1071 | 464 | 1214 | 583 | 458 | |

NC stretch | 2119 | 2097 | 2105 | 2101 | 2067 | |

FHF^{−} | FH stretch | 2260 | 2114 | 2271 | 1826 | 1645 |

FH bend | 1768 | 1451 | 1774 | 1286 | 1283 | |

FF stretch | 653 | 639 | 640 | 640 | 556 |

. | vibrational mode . | FPB . | TPB . | NEO-DFT(V)^{a}
. | NEO-DFT(V)-aug^{b}
. | VPT2 . |
---|---|---|---|---|---|---|

HCN | CH stretch | 3844 | 3521 | 3910 | 3297 | 3319 |

CH bend | 1135 | 866 | 1297 | 797 | 744 | |

CN stretch | 2278 | 2222 | 2212 | 2191 | 2173 | |

HNC | NH stretch | 3970 | 3791 | 3991 | 3635 | 3604 |

NH bend | 1071 | 464 | 1214 | 583 | 458 | |

NC stretch | 2119 | 2097 | 2105 | 2101 | 2067 | |

FHF^{−} | FH stretch | 2260 | 2114 | 2271 | 1826 | 1645 |

FH bend | 1768 | 1451 | 1774 | 1286 | 1283 | |

FF stretch | 653 | 639 | 640 | 640 | 556 |

^{a}

The NEO-DFT(V) calculations were performed with the same basis sets as the NEO-Ehrenfest dynamics with the proton basis function centers placed at the variationally optimized position.

^{b}

For the NEO-DFT(V)-aug calculations, the NEO-TDDFT calculations used to obtain the hydrogen vibrational excitations in the extended Hessian were performed using the Cartesian cc-pVTZ electronic basis set for the heavy nuclei and the Cartesian cc-pV6Z^{*} (no h function) electronic basis set for the quantum proton, as well as the Cartesian 8s8p8d8f protonic basis set, and the proton basis function centers were placed at the expectation value determined by a full geometry optimization with these basis sets, as prescribed in previous studies.^{44,52}

The vibrational frequencies in Tables I and II exhibit several general trends. First, for modes that are not directly related to the quantum proton, such as the C–N, N–C, and F–F stretch modes, NEO-Ehrenfest with or without a traveling proton basis yields similar results as NEO-HF/DFT(V) and VPT2. This behavior is expected because the locality of the quantum proton makes it unlikely to affect these modes significantly. Second, NEO-Ehrenfest with fixed proton basis functions produces similar frequencies for hydrogen stretch and bend modes as NEO-HF/DFT(V) with the same basis sets, and both significantly deviate from the VPT2 results. The use of larger basis sets with NEO-HF/DFT(V) leads to significantly more accurate hydrogen vibrational frequencies similar to those obtained with VPT2. NEO-Ehrenfest with a traveling proton basis decreases the predicted frequencies and brings them closer to the VPT2 results mainly because the traveling proton basis can more effectively explore regions outside the fixed basis function limit. The use of a larger traveling proton basis would lead to further improvements.

## IV. CONCLUSION

We have presented an Ehrenfest dynamics algorithm that utilizes real-time time-dependent nuclear–electronic orbitals to capture the interplay between nuclear and electronic dynamics. In the NEO-Ehrenfest approach, the quantum subsystem, including all electrons and selected protons, is simulated using the time-dependent Schrödinger equation, whereas the remaining molecular degrees of freedom move classically according to the mean NEO-TDHF/NEO-TDDFT potential. To account for a finite proton basis, a semiclassical approach is developed to allow the proton basis functions to travel during the dynamics. In benchmark studies of vibrational spectra for several small molecules, the NEO-Ehrenfest approach with traveling proton basis function centers yields reasonably accurate vibrational frequencies, and this level of accuracy will improve with larger electronic and nuclear basis sets.

NEO-Ehrenfest is a powerful method to resolve spectroscopic and dynamical properties of quantum mechanically coupled electrons and nuclei. Although the semiclassical treatment of the traveling proton basis has proven to be an effective method to account for the basis set incompleteness for computing vibrational spectra, it is still an approximate approach. It can be further improved with more rigorously derived approaches that are currently under investigation. Moreover, applications to proton transfer reactions will provide more stringent tests of the traveling proton basis function centers. This work forms the foundation for these future developments.

## SUPPLEMENTARY MATERIAL

See the supplementary material for the molecular geometries used in this work.

## ACKNOWLEDGMENTS

The development of quantum dynamics for studies of proton transfer was supported by IDREAM (Interfacial Dynamics in Radioactive Environments and Materials), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES). The development of the Chronus Quantum open source software package was supported by the National Science Foundation (Grant No. OAC-1663636). The development of the real-time electronic structure method was funded by the U.S. Department of Energy (Grant No. DE-SC0006863). The development of nuclear–electronic orbital time-dependent density functional theory methods was supported by the National Science Foundation (Grant No. CHE-1954348). Computations were facilitated through the use of advanced computational, storage, and networking infrastructure provided by the Hyak supercomputer system at the University of Washington, funded by the Student Technology Fee. L.Z. acknowledges the Dalton Postdoc Fellowship at the University of Washington for funding. We also thank John Tully for insightful discussions.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.