13.5.2 Overview of Available NEO Models

The simplest method within the NEO framework is the Hartree-Fock (NEO-HF) method, where the total nuclear-electronic wavefunction is approximated as a product of electronic (${𝚽}_0^{\text{e}}({𝐱}_{\text{e}})$) and nuclear (${𝚽}_0^{\text{p}}({𝐱}_{\text{p}})$) Slater determinants composed of electronic and protonic spin orbitals, respectively:

Here, ${𝐱}_{\text{e}}$ and ${𝐱}_{\text{p}}$ are collective spatial and spin coordinates of the quantum electrons and protons. The NEO-HF energy for a restricted Hartree-Fock (RHF) treatment of the electrons and a high-spin open-shell treatment of the quantum protons is

The $i,j,\mathrm{\cdots }$, indices denote occupied spatial electronic orbitals, and the $I,J,\mathrm{\cdots }$, indices correspond to occupied spatial protonic orbitals. In Eq. (13.34), $h_{ij}^{\text{e}}$ and $(ij|kl)$ are conventional electronic core Hamiltonian and two-electron integrals, respectively, and the corresponding terms for quantum protons are defined analogously. The last term in Eq. (13.34) is the Coulomb interaction between the electrons and the quantum protons. The spatial electronic and protonic orbitals [${\psi }_i^{\text{e}}({𝐫}_{\text{e}})$ and ${\psi }_I^{\text{p}}({𝐫}_{\text{p}})$] are expanded as linear combinations of electronic or protonic Gaussian basis functions [${\varphi }_{\mu }^{\text{e}}({𝐫}_{\text{e}})$ or ${\varphi }_{{\mu }^{\prime }}^{\text{p}}({𝐫}_{\text{p}})$]:

The lower-case Greek letters without and with primes denote basis functions for electrons and protons, respectively, and $C_{\mu i}^{\text{e}}$ and $C_{{\mu }^{\prime }I}^{\text{p}}$ are electronic and protonic MO expansion coefficients, respectively.

Analogous to the conventional electronic Hartree-Fock method, the electronic and protonic coefficients are determined by variationally minimizing the energy in Eq. (13.34) via the self-consistent field (SCF) procedure. This procedure leads to a set of coupled electronic and protonic Roothaan equations:

where ${𝐒}^{\text{e}}$ and ${𝐒}^{\text{p}}$ are electronic and protonic overlap matrices, respectively. The electronic and protonic Fock elements in Eqs. (13.36a) and (13.36b) are given by

The electronic and protonic density matrix elements in Eqs. (13.37a) and (13.37b) are defined as

The generalization to the unrestricted Hartree-Fock (NEO-UHF) treatment of electrons is accomplished by introducing separate spatial orbitals for $\alpha$ and $\beta$ electron spins.

The analytical gradients of the NEO-HF energy ¹³³⁵ Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106. Link with respect to the classical nuclear coordinates (or coordinates of the centers of the quantum proton basis functions) are available. These gradients allow geometry optimizations within the NEO framework. The analytical Hessians of the NEO-HF energy with respect to the classical nuclear coordinates are also available. ¹¹³⁴ Schneider P. E. et al.
J. Chem. Phys.
(2021), 154, pp. 054108. Link The Hessians can identify whether the optimized geometries are minima or transition states on the ground state vibronic potential energy surface.

NEO density functional theory (NEO-DFT) ⁹⁶⁶ Pak M. V., Chakraborty A., Hammes-Schiffer S.
J. Phys. Chem. A
(2007), 111, pp. 4522. Link ^{, 213} Chakraborty A., Pak M. V., Hammes-Schiffer S.
Phys. Rev. Lett.
(2008), 101, pp. 153001. Link ^{, 214} Chakraborty A., Pak M. V., Hammes-Schiffer S.
J. Chem. Phys.
(2009), 131, pp. 124115. Link is an extension of DFT to multicomponent systems within the NEO framework. The Hohenberg-Kohn theorems have been extended to multicomponent systems, where the reference is expressed as the product of electronic and nuclear Slater determinants composed of Kohn-Sham orbitals. The NEO-DFT total energy is

In this equation, $E_{\text{ext}}[{\rho }^{\text{e}},{\rho }^{\text{p}}]$ is the interaction of the electronic and protonic densities with the external potential created by the classical nuclei, and $E_{\text{ref}}[{\rho }^{\text{e}},{\rho }^{\text{p}}]$ contains the electron–electron, proton–proton, and electron-proton classical Coulomb energies, as well as the noninteracting kinetic energies of the quantum particles. The terms $E_{\text{exc}}[{\rho }^{\text{e}}]$, $E_{\text{pxc}}[{\rho }^{\text{p}}]$, and $E_{\text{epc}}[{\rho }^{\text{e}},{\rho }^{\text{p}}]$ are the electron-electron exchange-correlation functional, the proton-proton exchange-correlation functional, and the electron-proton correlation functional, respectively. The quantities

are the electron and proton densities, respectively, and ${\psi }_i^{\text{e}}({𝐫}_1^{\text{e}})$ and ${\psi }_I^{\text{p}}({𝐫}_1^{\text{p}})$ are the electronic and protonic Kohn-Sham spatial orbitals, respectively. These orbitals are obtained by solving two sets of coupled Kohn-Sham equations for the electrons and quantum protons:

The effective potentials $v_{\text{eff}}$ and $v_{\text{eff}}$ are obtained by taking the derivative of the total energy expression in Eq. (13.39) with respect to electron density and proton density, respectively. Analogous to NEO-HF, these electronic and protonic Kohn-Sham orbitals are expanded as linear combinations of electronic or protonic Gaussian basis functions (${\varphi }_{\mu }^{\text{e}}({𝐫}_{\text{e}})$ and ${\varphi }_{{\mu }^{\prime }}^{\text{p}}({𝐫}_{\text{p}})$). The extension to open-shell electron systems is analogous to the NEO-UHF method.

The practical implementation of the NEO-DFT method requires an electron-electron exchange-correlation functional, a proton-proton exchange-correlation functional, and an electron-proton correlation functional. Any conventional electron-electron exchange-correlation functional can be used within the NEO-DFT framework. ¹⁵⁰ Brorsen K. R., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2018), 149, pp. 044110. Link Because the proton-proton exchange and correlation are negligible in molecular systems, only the exchange at the NEO-Hartree-Fock level is included to eliminate self-interaction error in the NEO-DFT method. A suitable electron-proton correlation functional is essential for obtaining accurate proton densities and energies, and the epc17-2 ¹⁴¹² Yang Y. et al.
J. Chem. Phys.
(2017), 147, pp. 114113. Link ^{, 151} Brorsen K. R., Yang Y., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2017), 8, pp. 3488. Link and epc19 ¹²⁵⁸ Tao Z., Yang Y., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 124102. Link functionals are designed to achieve this goal. These two functionals are based on the multicomponent extension of the Colle-Salvetti formalism. The epc17-2 functional is of the local density approximation (LDA) type with the functional form:

The epc19 functional is its multicomponent generalized gradient approximation (GGA) extension that depends on the electron and proton density gradients and is of the form:

In addition to the parameters $a$, $b$, and $c$ in the epc17-2 functional, ¹⁵¹ Brorsen K. R., Yang Y., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2017), 8, pp. 3488. Link the epc19 functional ¹²⁵⁸ Tao Z., Yang Y., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 124102. Link has the $d$ and $k$ parameters and also depends on the proton mass $m_{\text{p}}$. Analogous to the NEO-HF analytical energy gradients, the NEO-DFT analytical gradients ¹²⁵⁷ Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110. Link are available for these two functionals, allowing geometry optimizations on the ground state vibronic potential energy surface. The NEO-DFT analytical Hessians ¹²⁵⁷ Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110. Link are available for the epc17-2 functional or when no electron-proton correlation functional is used and allow characterization of the stationary points.

The NEO multistate DFT (NEO-MSDFT) method was developed to describe hydrogen transfer and hydrogen tunneling systems within the NEO framework. ¹⁴²⁴ Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106. Link ^{, 307} Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170. Link Analogous to the conventional electronic MSDFT method, ⁸⁹⁴ Mo Y., Bao P., Gao J.
Phys. Chem. Chem. Phys.
(2011), 13, pp. 6760. Link ^{, 396} Gao J. et al.
J. Phys. Chem. Lett.
(2016), 7, pp. 5143. Link ^{, 461} Grofe A. et al.
J. Chem. Theory Comput.
(2017), 13, pp. 1176. Link ^{, 815} Lu Y., Gao J.
J. Phys. Chem. Lett.
(2022), 13, pp. 7762. Link the NEO-MSDFT method linearly combines localized NEO-DFT states in a nonorthogonal configuration interaction (NOCI) scheme. ¹¹⁷⁷ Skone J. H., Pak M. V., Hammes-Schiffer S.
J. Chem. Phys.
(2005), 123, pp. 134108. Link This approach captures the delocalized, bilobal proton densities needed to describe hydrogen tunneling and avoids complications of local minima in orbital space.

Consider a system with $N$ transferring protons, where each proton moves in a double-well potential in the Born-Oppenheimer picture. In NEO-DFT calculations, the protonic density tends to localize in one well of a symmetric or nearly symmetric double-well potential instead of delocalizing over both wells. Quantizing each transferring proton with the NEO approach leads to $2^N$ diabatic NEO-DFT states. Each diabatic state has the protonic density of each transferring proton localized in one of the two wells of its corresponding double-well potential. In practice, higher-energy diabatic states can be excluded from the NOCI expansion. The set of all diabatic NEO-DFT states is $\{\text{|}{\tilde{\mathrm{\Psi }}}_0⟩,\text{|}{\tilde{\mathrm{\Psi }}}_1⟩,\mathrm{\dots }\text{|}{\tilde{\mathrm{\Psi }}}_n⟩\}$, where $n=2^N-1$. Each diabatic state is the product of a Kohn-Sham electronic and protonic determinant, as discussed in Section 13.5.2.2.

The adiabatic NEO-MSDFT states $\{\text{|}{\mathrm{\Psi }}_0⟩,\text{|}{\mathrm{\Psi }}_1⟩,\mathrm{\dots }\text{|}{\mathrm{\Psi }}_n⟩\}$ are linear combinations of all diabatic NEO-DFT states:

The coefficients in Eq. (13.44) are determined by solving the $2^N\times 2^N$ matrix eigenvalue problem

The overlap matrix $𝐒$ and effective Hamiltonian matrix $𝐇$ contain the overlap and couplings, respectively, between pairs of localized diabatic states. Note that the diagonal elements of $𝐒$ are unity, and the diagonal elements of $𝐇$ are the NEO-DFT energies of the diabatic states. For the sake of brevity, the analytical forms of the off-diagonal terms of these matrices are excluded here but are given in previous work. ³⁰⁷ Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170. Link

The limitations of the epc functionals and the resulting inaccuracies of the overlap between two localized NEO-DFT states, as well as the approximate form of the off-diagonal Hamiltonian matrix elements, can be accounted for by applying a simple correction function to the off-diagonal elements of the $𝐒$ matrix:

Users can control the values of $\alpha$ and $\beta$ used in NEO-MSDFT calculations, and their default values were determined though a fitting process discussed in our previous work. ¹⁴²⁴ Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106. Link Setting both these parameters to unity produces the original overlap matrix. However, in general the alpha and beta parameters should be kept at their default values unless another parameterization procedure is performed to determine alternative values.

A NEO-MSDFT calculation can be enabled by setting the $rem variable NEO_MSDFT = 1. Analytical gradients have also been implemented, ¹⁴²⁶ Yu Q., Schneider P. E., Hammes-Schiffer S.
J. Chem. Phys.
(2022), 156, pp. 114115. Link allowing users to optimize geometries on specified NEO-MSDFT surfaces by setting JOBTYPE = OPT. By default, analytical derivative couplings between NEO-MSDFT states are computed directly following an analytical gradient calculation, ¹⁴²⁵ Yu Q., Roy S., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2022), 18, pp. 7132. Link but this calculation can be disabled. Note that the only gradients and derivative couplings printed out are those pertaining to the ground and first-excited NEO-MSDFT states, as these are the only states of physical relevance (as explained elsewhere ³⁰⁷ Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170. Link ). Users can also request a semi-numerical Hessian calculation by setting IDERIV = 2. Further control of a NEO-MSDFT calculation is given by the $neo_msdft section, whose job-control options are given below:

As discussed in previous work, ¹⁴²⁴ Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106. Link ^{, 307} Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170. Link each quantum transferring proton in NEO-MSDFT is given two basis function centers: one center localized near the donor and another center localized near the acceptor. When requesting a NEO-MSDFT calculation, in the $molecule  section, each quantum transferring proton must have one center input as a standard H atom center and its other center input as a ghost H atom center. These centers must be listed consecutively in the $molecule section. Additionally, it is assumed that all ghost H atom centers in the $molecule section of a NEO-MSDFT calculation belong to a quantum transferring proton.

By default, all possible $2^N$ NEO-DFT diabatic states are included in the adiabatic state expansion of Eq. (13.44). If this is not desired, users can provide an input to the $neo_msdft_diabat_control section to specify which diabatic states to include. This input section takes in an array of Boolean values (input as ones or zeros), where each row is a unique string corresponding to a particular diabatic state. These Booleans refer to whether or not a quantum transferring proton is localized on each proton basis function center in the order in which they appear in the $molecule section. Thus, if the number of diabatic states the user wants to include is $m$ and the number of quantum transferring protons is $N$, the $neo_msdft_diabat_control section should contain an array with $m$ rows and $2N$ columns, where the $i$-th row specifies which proton basis function centers have a localized quantum proton in the $i$-th diabatic state. Note that other protons can be quantized with NEO that do not correspond to a transferring proton and therefore are represented by a single proton basis function center; these protons should not be included in the $neo_msdft_diabat_control section. Examples showing how to set up NEO-MSDFT calculations/input files with both one and multiple transferring protons can be found in Section 13.5.4.

The overlap correction of Eq. (13.46) is parameterized to produce accurate proton densities and tunneling splittings when the proton basis function centers are optimized for separate NEO-DFT diabatic states (see the procedures explained in our previous work). ¹⁴²⁴ Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 10106. Link ^{, 307} Dickinson J. A., Yu Q., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2023), 14, pp. 6170. Link However, users may wish to perform geometry optimizations on NEO-MSDFT surfaces instead. One can either optimize the positions of just the proton basis function centers for fixed geometries of classical nuclei or optimize the positions of both the classical nuclei and proton basis function centers simultaneously. Both cases require the user to set JOBTYPE = OPT, but in the former case, the $opt section must be used to freeze the classical nuclei during the optimization (see Section 9.4.3 for details). In either case, the input of the $molecule section serves as an initial guess for the optimization, so users must make an informed guess as to where to initially place each of the two proton basis function centers for each quantum transferring proton.

Bulk solvent effects can be directly incorporated into NEO calculations through the application of various implicit solvation models (Section 11.2) within the NEO framework. ¹³⁶⁷ Wildman A. et al.
J. Chem. Theory Comput.
(2022), 18, pp. 1340. Link The polarizable continuum model (PCM) constitutes one family of implicit solvation models and itself encompasses several different formulations: ⁷¹⁶ Lange A. W., Herbert J. M.
Chem. Phys. Lett.
(2011), 509, pp. 77. Link ^{, 529} Herbert J. M.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021), 11, pp. e1519. Link C-PCM, ¹²⁸² Truong T. N., Stefanovich E. V.
Chem. Phys. Lett.
(1995), 240, pp. 253. Link ^{, 72} Barone V., Cossi M.
J. Phys. Chem. A
(1998), 102, pp. 1995. Link IEF-PCM, ²²⁸ Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558. Link ^{, 171} Cancès E., Mennucci B.
J. Chem. Phys.
(2001), 114, pp. 4744. Link SS(V)PE, ²²⁸ Chipman D. M.
J. Chem. Phys.
(2000), 112, pp. 5558. Link etc. ⁵²⁹ Herbert J. M.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021), 11, pp. e1519. Link In the PCM approach, the solute molecule is placed in a cavity that is embedded in dielectric continuum solvent, and the cavity surface is discretized into $i$ tesserae grid points. The solvent response is represented by a partial charge $q_i$ centered at each tesserae grid point ${𝐬}_i$. ¹²⁷⁰ Tomasi J., Mennucci B., Cammi R.
Chem. Rev.
(2005), 106, pp. 2999. Link ^{, 529} Herbert J. M.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021), 11, pp. e1519. Link

For each SCF iteration, the current electronic and protonic densities, together with the fixed classical nuclei, define the solute’s charge distribution. This charge distribution gives rise to the solute’s electrostatic potential $V_i$ at each tesserae grid point:

The solute electrostatic potential is then used to compute $q_i$ using standard PCM methods. Once obtained, the set of tesserae charges is included as an additional one-electron (one-proton) contribution to the electronic (protonic) Fock or analogous Kohn-Sham matrix:

where $F_{\mu \nu }^{\text{e},\text{0}}$ and $F_{{\mu }^{\prime }{\nu }^{\prime }}^{\text{p},\text{0}}$ refer to the gas-phase, electronic [Eq. (13.37a)] and protonic [Eq. (13.37b)] Fock or analogous Kohn-Sham matrix elements, respectively.

NEO-PCM calculations involve iterative, self-consistent convergence of the nuclear-electronic wavefunction in the presence of the dielectric continuum solvent. Both NEO-HF and NEO-DFT PCM energies and analytic gradients are implemented. The calculation can be invoked by setting SOLVENT_METHOD = PCM in the $rem input section, alongside variables for NEO-SCF as described in Section 13.5.3.

In the simplest approach, the cavity surface is discretized into point charges. However, a more sophisticated approach utilizing Gaussian-smeared charges is also supported. ⁷¹³ Lange A. W., Herbert J. M.
J. Phys. Chem. Lett.
(2010), 1, pp. 556. Link ^{, 714} Lange A. W., Herbert J. M.
J. Chem. Phys.
(2010), 133, pp. 244111. Link ^{, 711} Lange A. W. et al.
Mol. Phys.
(2020), 118, pp. e1644384. Link ^{, 529} Herbert J. M.
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
(2021), 11, pp. e1519. Link Selection of these various PCM schemes and related variables can be set in the $pcm and/or $solvent input sections. An example on how to set up a solvated NEO calculation can be found in Section 13.5.4.

In NEO-SCF calculations, the nuclear solution corresponds to the variational nuclear wavefunction (density) for the specific method and basis sets. The nuclear orbitals in molecular systems are localized, and therefore the quantized nuclei can be treated as distinguishable when analyzing their position expectation values. These expectation values are generally not equal to the input coordinates of the NEO nuclear basis function centers. Control over the position expectation values of the nuclei can be achieved with the constrained (C)NEO approach, ¹⁴⁰² Xu X., Yang Y.
J. Chem. Phys.
(2020), 152, pp. 084107. Link based on the original NEO method ¹³³⁵ Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106. Link and in analogy to the conventional constrained DFT (Section 5.11). This is accomplished by imposing the following constraint:

where ${𝐑}_I$ is the chosen value for the position expectation value for nucleus $I$. For geometry optimizations and vibrational frequency calculations, typically ${𝐑}_I$ is the position of the NEO nuclear basis function center, as defined by the input coordinates. The constraint is imposed by using Lagrange multipliers. In this case, the protonic part of the coupled Kohn-Sham equations of NEO-DFT (Eq. (13.41b)) becomes:

where ${𝐟}_I$ is the Lagrange multiplier. In practice, the value of ${𝐟}_I$ is optimized numerically to satisfy the constraint in Eq. (13.49).

CNEO calculations can be invoked by setting CNEO = TRUE, and NEO = TRUE must also be set. Analytical CNEO gradients ¹⁴⁰³ Xu X., Yang Y.
J. Chem. Phys.
(2020), 153, pp. 074106. Link and Hessians ¹⁴⁰⁴ Xu X., Yang Y.
J. Chem. Phys.
(2021), 154, pp. 244110. Link are available, allowing geometry optimizations with JOBTYPE = OPT and vibrational frequency calculations with JOBTYPE = FREQ. CNEO can be used with Hartree-Fock or with DFT in conjunction with the epc17-2 or epc19 electron-proton correlation functional and the electron exchange-correlation functionals available in Q-Chem. CNEO is currently implemented with only one quantum nucleus.

The LR-NEO-TDDFT method ¹⁴¹³ Yang Y., Culpitt T., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2018), 9, pp. 1765. Link is a multicomponent extension of the linear response TDDFT method within the NEO framework. It allows the simultaneous calculation of the electronic and proton vibrational excitation energies. In the (LR)-NEO-TDDFT method, the linear response of the NEO Kohn-Sham system to perturbative external fields is computed. The NEO-TDDFT working equation is

where

Here, the occupied electronic orbitals are denoted with indices $i$ and $j$, whereas the unoccupied electronic orbitals are denoted with indices $a$ and $b$. The analogous upper case indices denote protonic orbitals. The solution of Eq. (13.51) provides the electronic and proton vibrational excitation energies $\omega$, as well as the transition excitation and de-excitation amplitudes, $𝐗$ and $𝐘$, respectively. Analogous to the TDDFT method, the Tamm-Dancoff approximation (TDA) can be imposed within the NEO framework, defining the NEO-TDDFT-TDA method that is represented by

The extension of the NEO-TDDFT and NEO-TDDFT-TDA approaches to open-shell electron systems is straightforward. ²⁷¹ Culpitt T. et al.
J. Chem. Phys.
(2019), 150, pp. 201101. Link NEO-TDHF and NEO-CIS have similar forms as NEO-TDDFT and NEO-TDA without electron-proton, electron-electron, or proton-proton correlation. The analytical gradients for NEO-CIS/NEO-TDA/NEO-TDHF/NEO-TDDFT are available, ¹²⁵⁷ Tao Z. et al.
J. Chem. Theory Comput.
(2021), 17, pp. 5110. Link enabling geometry optimizations on the excited state vibronic potential energy surfaces. For NEO-TDA and NEO-TDDFT, analytical gradients are available for the epc17-2 functional or when no electron-proton correlation functional is used. The transition densities can be analyzed to determine the percentages of electronic and protonic character for each vibronic excited state.

Real-time NEO time-dependent density functional theory (RT-NEO-TDDFT, or RT-NEO for brevity) is a multicomponent extension of conventional electronic RT-TDDFT within the NEO framework ¹⁴⁴⁰ Zhao L. et al.
J. Phys. Chem. Lett.
(2020), 11, pp. 4052. Link and is an approach that enables the modeling of nonequilibrium, non-Born-Oppenheimer nuclear-electronic quantum dynamics. This approach assumes the form of the single Slater determinant product ansatz of the NEO-DFT reference system, where ${𝐱}^{\text{e}}$ and ${𝐱}^{\text{p}}$ are collective spatial and spin coordinates of the quantum electrons and protons, ${𝐏}^{\text{e}}$ (or ${𝐏}^{\text{p}}$) is the single particle density matrix and ${𝐅}^{\text{e}}$ (or ${𝐅}^{\text{p}}$) is the Kohn-Sham matrix for the electrons (or protons), and $t$ is the time. The time-dependent Schrödinger equation

can be propagated according to the set of multicomponent von Neumann equations:

Note that these equations are not propagated independently, but are actually coupled through Fock terms that depend on both particle densities. The electron and proton densities are evolved together in time. The resulting electronic and protonic dipole moments can be obtained at every time step and analyzed via post-processing to provide relevant spectral features. Fourier transformation of these signals from the corresponding time-domain to the frequency-domain can yield both electronic and vibrational spectra that have been shown to agree well with LR-NEO-TDDFT. In addition, this approach provides the nonequilbrium, real-time dynamics of the electronic and protonic densities.

Up to this point, the equations of motion presented above [Eq. (13.59a) and Eq. (13.59b) assume that the classical nuclei are held at a fixed geometry throughout the entire simulation. To combine the real-time dynamics of the electrons and protons together with the mean-field motion of the classical nuclei, the RT-NEO-TDDFT-Ehrenfest approach has been formulated. ¹⁴⁴¹ Zhao L. et al.
J. Chem. Phys.
(2020), 153, pp. 224111. Link In this approach, the classical nuclei are evolved according to Newtonian dynamics, while the nonequilibrium electrons and protons are still being propagated according to Eq. (13.59a) and Eq. (13.59b), respectively:

Here we denote $M_A$ and ${𝐑}_A$ to be the mass and the position, respectively, of the $A$th classical nucleus, and ${\nabla }_A\text{𝐸}$ is the abbreviated form of the total nuclear gradient expression. We support both of these approaches, i.e., RT-NEO with fixed classical nuclei and RT-NEO with moving classical nuclei, and we have designated those methods as METHOD = REALTIME and METHOD = EHRENFEST, respectively.

In addition to electrons, it is important that the proton dynamics is adequately described throughout the course of a given trajectory. This is especially important in cases where proton delocalization may start to occur, such as during a proton transfer reaction or when the geometry of the classical nuclei is no longer held fixed. To provide the flexibility for describing proton delocalization, one may use a larger protonic basis set and/or utilize multiple fixed proton basis (FPB) function centers placed at pre-designated positions. Alternatively, a traveling proton basis (TPB) approach has been formulated and will be made available in a future release. ¹⁴⁴¹ Zhao L. et al.
J. Chem. Phys.
(2020), 153, pp. 224111. Link Examples of using either one or multiple FPB function centers can be found in the examples section of this manual.

In certain applications, the requirement of using a small time step due to the fast electronic dynamics can prohibit long simulation times. For systems that are electronically adiabatic, the electronic Born-Oppenheimer (BO) approximation can be invoked by quenching the electrons to an instantaneous electronic ground state via the standard self-consistent-field (SCF) procedure instead of propagating Eq. (13.59a) at every time step:

This approach enables the use of a time step that is several orders of magnitude larger than that required for integrating Eq. (13.59a), therefore rendering longer simulation times more tractable. The electronic BO approximation can be invoked for RT-NEO with fixed classical nuclei or for RT-NEO with moving classical nuclei by toggling the methods METHOD = BO-REALTIME or METHOD = BO-EHRENFEST, respectively.

A RT-NEO calculation can be invoked by setting NEO_TDKS = TRUE in the $rem section, and additional job control is provided in its own input section. The input file for a RT-NEO propagation is illustrated in several examples in Section 13.5.4, however, the available controls are also discussed below. Note that the keywords and field types defined for $tdks (Section 7.4) are not the same parameters as those defined for the $neo_tdks, and vice versa. Only parameters/keywords documented below are supported. The parameters of the jobs are set using options specified in the $neo_tdks input section. The format of the $neo_tdks section is analogous to the $rem section:

$neo_tdks
  <Keyword>  <parameter/option>
$end

Note:  The following job control variables belong only in the $neo_tdks  section. Do not place them in the $rem section.

To perturb the system, a time-dependent electric field pulse $𝓔(t)$ can be applied via a dipole coupling term added to the electronic and/or protonic Fock matrices,

where ${𝐅}^{\text{e},\text{0}}$ (or ${𝐅}^{\text{p},\text{0}}$) is the unperturbed Kohn-Sham matrix and ${𝝁}^{\text{e}}$ (or ${𝝁}^{\text{p}}$) is the dipole moment matrix for the electrons (or protons). Note that the dipole moment (i.e., ${𝝁}^{\text{e}}$ and ${𝝁}^{\text{p}}$) have the opposite sign for electrons and protons. Tunable parameters for FIELD_TYPE are provided in the additional keywords found below. For the supported field types below, note that the electric field vector $𝓔=({ℰ}_x,{ℰ}_y,{ℰ}_z)$ is a quantity whose magnitude is controlled by FIELD_AMP and whose direction is controlled by FIELD_DIRECTION.

1. 1.

   Delta simulates a Dirac $\delta$-function kick with a field that is turned on only at time $t_0=0$:

   $$𝓔(t)={𝐀}_0\delta \left(t-t_0\right)$$

   (13.63)

   • •

     The amplitude of ${𝐀}_0$ and its direction are set using FIELD_AMP and FIELD_DIRECTION, respectively.

2. 2.

   Gaussian simulates the following impulse field that is turned on at all times $t\ge 0$:

   $$𝓔(t)={𝐀}_0\exp \left(\frac{{(t-t_{\text{peak}})}^2}{{\tau }^2}\right)\sin (\omega t)$$

   (13.64)

   • •

     The amplitude of ${𝐀}_0$ and its direction are set using FIELD_AMP and FIELD_DIRECTION, respectively.

   • •

     The center of the pulse $t_{\text{peak}}$ is set using the FIELD_PEAK keyword.

   • •

     The pulse half-width $\tau$ is set using the FIELD_TAU keyword.

   • •

     The pulse frequency $\omega$ is set using the FIELD_FREQUENCY keyword.

Within the NEO framework, select nuclei are treated quantum mechanically at the same level as the electrons. This removes the Born-Oppenheimer separation between the quantum nuclei and the electrons and naturally includes nonadiabatic effects between the quantum nuclei and the electrons. At the same time, quantizing the select nuclei gives rise to a potential energy surface with fewer nuclear degrees of freedom, which prevents a direct calculation of the vibrational frequencies of the entire molecule. Consequently, diagonalization of a coordinate Hessian in the NEO framework yields vibrational frequencies and accompanying normal modes of only the classical nuclei, with the quantum nuclei responding instantaneously to the motion of the classical nuclei. ¹¹³⁴ Schneider P. E. et al.
J. Chem. Phys.
(2021), 154, pp. 054108. Link Although the fundamental anharmonic vibrational frequencies of the quantum nuclei can be accurately obtained through LR-NEO-TDDFT, ²⁷¹ Culpitt T. et al.
J. Chem. Phys.
(2019), 150, pp. 201101. Link the couplings between the vibrations of the classical and quantum nuclei are missing. To obtain the fully coupled molecular vibrations, an effective strategy denoted NEO-DFT(V) was developed. ¹⁴¹⁴ Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167. Link ^{, 272} Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840. Link The NEO-DFT(V) method has been shown to incorporate key anharmonic effects in full molecular vibrational analyses and to produce accurate molecular vibrational frequencies compared to experiments. ¹⁴¹⁴ Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167. Link ^{, 272} Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840. Link

The NEO-DFT(V) method involves diagonalization of an extended NEO Hessian composed of partial second derivatives of the coordinates of the classical nuclei (${𝐫}_c$) and the expectation values of the quantum nuclei (${𝐫}_q$). This extended Hessian matrix is composed of three sub-matrices: ${\text{𝐇}}_0={({\partial }^{2}E/\partial {𝐫}_c^2)|}_{{𝐫}_c={𝐫}_q}$, ${\text{𝐇}}_1={\partial }^{2}E/\partial {𝐫}_q\partial {𝐫}_c$, and ${\text{𝐇}}_2={\partial }^{2}E/\partial {𝐫}_q^2$, where in each case, all other coordinates of the classical nuclei and expectation values of the quantum nuclei are fixed. The extended Hessian has the following structure:

where

The quantity $𝐑=d{𝐫}_q/d{𝐫}_c$ and the NEO Hessian matrix is ${𝐇}_{\text{NEO}}=d^{2}E/d{𝐫}_c^2$ (without the constraint that the expectation values of the quantum nuclei are fixed). In the expression for the ${𝐇}_{\text{2}}$ matrix, $𝐌$ is the diagonal mass matrix, and $𝛀$ is the diagonal matrix with elements corresponding to the squares of the LR-NEO-TDDFT fundamental vibrational frequencies. ¹⁴¹⁴ Yang Y. et al.
J. Phys. Chem. Lett.
(2019), 10, pp. 1167. Link $𝐔$ is a unitary matrix that transforms $𝛀$ to the target coordinate system and is approximated with the transition dipole moment vectors afforded by a LR-NEO-TDDFT calculation. ²⁷² Culpitt T. et al.
J. Chem. Theory Comput.
(2019), 15, pp. 6840. Link

Diagonalization of ${𝐇}_{\text{DFT(V)}}$ produces the fully coupled molecular vibrational frequencies including anharmonic effects associated with the quantum protons. The NEO-DFT(V) method is available for use with the epc17-2 functional or when no electron-proton correlation functional is used. The NEO-HF(V) method, which involves building the extended NEO-Hessian based on the NEO-HF Hessian and inputs from NEO-TDHF, is also available.

An alternative route for inclusion of correlation effects between quantum particles (i.e., electrons and protons) is with wave functions methods that are systematically improvable and parameter-free. ¹³³⁵ Webb S. P., Iordanov T., Hammes-Schiffer S.
J. Chem. Phys.
(2002), 117, pp. 4106. Link ^{, 978} Pavošević F., Culpitt T., Hammes-Schiffer S.
Chem. Rev.
(2020), 120, pp. 4222. Link Among the various developed multicomponent wave function methods, the NEO coupled cluster (NEO-CC) methods have been particularly successful. ⁹⁷⁷ Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338. Link ^{, 979} Pavošević F., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 074104. Link ^{, 980} Pavošević F., Rousseau B. J. G., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 1578. Link ^{, 981} Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631. Link The NEO-CC wave function is given by

where $\widehat{T}$ is the cluster operator that incorporates the correlation effects between quantum particles, and $|0^{\text{e}}{0}^{\text{p}}⟩$ is the NEO-HF reference wave function. In the NEO-CCSD method, ⁹⁷⁷ Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338. Link the cluster operator is given by

where $a^{\alpha }=a_{\alpha }^{†}=\{a_i^a,a_I^A,a_{ij}^{ab},a_{IJ}^{AB},a_{iI}^{aA}\}$ are the excitation operators expressed in terms of creation/annihilation ($a_p^{†}/a_p$) fermionic operators, and $\alpha$ is the excitation rank. Here, the $i,j,\mathrm{\dots }$ indices denote occupied electronic orbitals, the $a,b,\mathrm{\dots }$ indices denote unoccupied electronic orbitals, and the $p,q,\mathrm{\dots }$ indices denote general electronic orbitals. The protonic orbitals are denoted analogously using the capitalized indices. The unknown $t_{\alpha }$ wave function parameters (amplitudes) are determined by solving the set of nonlinear equations for each $\alpha$: ⁹⁷⁷ Pavošević F., Culpitt T., Hammes-Schiffer S.
J. Chem. Theory Comput.
(2018), 15, pp. 338. Link

In this equation, $H_{\text{NEO}}=h_q^p{a}_p^q+\frac{1}{2}{g}_{rs}^{pq}{a}_{pq}^{rs}+h_Q^P{a}_p^q+\frac{1}{2}{g}_{RS}^{PQ}{a}_{PQ}^{RS}-g_{qQ}^{pP}{a}_{pP}^{qQ}$ is the second-quantized NEO Hamiltonian, where $h_q^p=⟨q|{\widehat{h}}^{\text{e}}|p⟩$ and $g_{rs}^{pq}=⟨rs|pq⟩$ are conventional electronic core Hamiltonian and two-electron integrals, respectively. The remaining protonic ($h_Q^P$ and $g_{RS}^{PQ}$) and electron-proton ($g_{qQ}^{pP}$) integrals are defined analogously. Lastly, the NEO-CCSD energy is calculated from

To increase the computational efficiency and reduce the memory requirements for the NEO-CCSD method, the two-particle integrals can be approximated with the density fitting (DF) approximation, ⁹⁸¹ Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631. Link in which the two-particle four-center integrals are factorized into a sum of products of three-center and two-center two-particle integrals. In particular, the four-center two-electron integrals are approximated by

where $(\mu \nu |X)$ and $(X|Y)$ are three-center and two-center two-electron integrals, respectively. In this equation, $\mu ,\nu ,\mathrm{\dots }$ and $X,Y,\mathrm{\dots }$ indices denote electronic and auxiliary electronic basis functions, respectively. The four-center two-proton integrals are approximated analogously by

where primed indices denote protonic basis functions and $({\mu }^{\prime }{\nu }^{\prime }|X)$ and $(X^{\prime }|Y^{\prime })$ are three-center and two-center two-proton integrals, respectively. Finally, the four-center electron-proton integrals are approximated as

By employing the DF approximation, the memory requirements for storing four-center two-particle integrals are reduced from $N_{\text{bf}}^4$ to $N_{\text{bf}}^2\times N_{\text{aux}}$, where $N_{\text{bf}}$ and $N_{\text{aux}}$ are the number of electronic or protonic basis functions and auxiliary basis functions, respectively. ⁹⁸¹ Pavošević F., Tao Z., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2021), 12, pp. 1631. Link

Second-order Møller-Plesset perturbation theory provides a useful framework for wave function based correlation effects with lower cost than NEO-CC ¹²⁴⁴ Swalina C., Pak M. V., Hammes-Schiffer S.
Chem. Phys. Lett.
(2005), 404, pp. 394. Link . The NEO-MP2 correlation energy is composed of the electron-electron, proton-proton and electron-proton contributions: