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.

The essential $rem variables for toggling a NEO-HF calculation are NEO = TRUE and METHOD = HF. Supported JOBTYPE’s include SP, FORCE, OPT, TS, and FREQ. Consult Section 13.5.4 for the full list of variables as well as see the examples in Section 13.5.6.

NEO density functional theory (NEO-DFT) ⁹⁹³ Pak M. V., Chakraborty A., Hammes-Schiffer S.
J. Phys. Chem. A
(2007), 111, pp. 4522. Link ^{, 223} Chakraborty A., Pak M. V., Hammes-Schiffer S.
Phys. Rev. Lett.
(2008), 101, pp. 153001. Link ^{, 224} 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 ^{, 161} 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 essential $rem variables for toggling a NEO-DFT calculation are NEO = TRUE and one of the two options: (1) specify METHOD or (2) specify EXCHANGE and CORRELATION. We inherit Q-Chem’s conventional, electronic DFT library, and therefore METHOD, EXCHANGE, and CORRELATION should be set to a functional listed in Section 5.3. Although it has not been tested across all possible combinations and permutations, NEO-DFT should be able to support the vast majority of the electronic XC functionals listed. Supported JOBTYPE’s include SP, FORCE, OPT, TS, and FREQ. Consult Section 13.5.4 for the full list of variables as well as see the examples in Section 13.5.6. Lastly, the electron-proton correlation functional should be set to one of the supported types:

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 ^{, 317} 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 ^{, 410} Gao J. et al.
J. Phys. Chem. Lett.
(2016), 7, pp. 5143. Link ^{, 477} Grofe A. et al.
J. Chem. Theory Comput.
(2017), 13, pp. 1176. Link ^{, 833} 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 in addition to the variables in the NEO-DFT subsection (see  13.5.2.2):

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 ^{, 317} 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.6.

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 ^{, 317} 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.

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

CNEO calculations can be invoked by setting CNEO = TRUE in addition to the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2):

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. Note that using CNEO with multiple quantum protons invokes the nuclear Hartree product approximation (Section 13.5.2.5).

In NEO-HF, the total wavefunction is a product of electronic and nuclear Slater determinants (Section 13.5.2.1); relatedly, in NEO-DFT, the reference is a product of Slater determinants (Section 13.5.2.2). Because the nuclear orbitals in molecular systems are highly localized and the overlap between them is effectively zero, the nuclear exchange interactions are numerically insignificant, ¹⁰⁰⁶ Pavošević F., Culpitt T., Hammes-Schiffer S.
Chem. Rev.
(2020), 120, pp. 4222. Link and we can treat the quantum nuclei as distinguishable. This means that we can approximate the nuclear wavefunction as a Hartree product instead of a Slater determinant: ⁵⁶ Auer B., Hammes-Schiffer S.
J. Chem. Phys.
(2010), 132, pp. 084110. Link

Here, the product is over $N_{\text{p}}$ quantum protons, ${𝐱}_{\text{p}}$ are the collective spatial and spin coordinates for the quantum protons, and ${\chi }_I^{\text{p}}({𝐱}_I^{\text{p}})$ is the spin orbital for the $I$th quantum proton. The total wavefunction is then (see also Eq. (13.33)):

Note that a high-spin configuration is assumed for the protonic wavefunction, and therefore the protonic spin orbitals can be replaced with spatial orbitals multiplied by the same spin function. For NEO-HF, Eq. (13.49) leads to $N_{\text{p}}$ coupled one-nucleus Roothaan equations instead of one $N_{\text{p}}$-nucleus Roothaan equation (Eq. (13.36b)). For NEO-DFT, it leads to $N_{\text{p}}$ coupled one-nucleus Kohn-Sham equations (Eq. (13.41b)). Making use of the localized nature of the quantum nuclei, the spatial nuclear orbitals for each nucleus (Eq. (13.35)) are expanded in the local basis set for that nucleus (i.e., different nuclei do not share basis functions).

The nuclear Hartree product approximation can be activated by setting NEO_N_HARTREE_PROD = TRUE in addition to the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2).

It is recommended to use the simultaneous NEO-SCF algorithms for faster convergence (Section 13.5.4.3). By default, the Coulomb and exchange self-interaction terms for each nucleus are computed. Similarly to NEO-SCF calculations with one quantum proton, the evaluation of these terms can be turned off with NEO_VPP = FALSE to further accelerate the calculations. Note that this does not change the NEO-SCF energy and the occupied orbitals, but it does affect the nuclear virtual orbitals. ¹²⁷⁸ Swalina C., Pak M. V., Hammes-Schiffer S.
Chem. Phys. Lett.
(2005), 404, pp. 394. Link Supported values of JOBTYPE are the same as for CNEO (Section 13.5.2.4) and include SP, FORCE, OPT, TS, and FREQ.

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.

The LR-NEO-TDHF (or TDDFT) method is activated by the following keyword:

Furthermore, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included. Supported JOBTYPE’s include SP and OPT. Additional customization can be controlled by the following $rem variables:

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:

The quantities $M_A$ and ${𝐑}_A$ are the mass and the position of classical nucleus $A$. 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 Section 13.5.6.

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: ⁷⁸⁸ Li T. E., Hammes-Schiffer S.
J. Chem. Phys.
(2023), 158, pp. 114118. Link

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-HF or RT-NEO-TDDFT calculation can be invoked by setting NEO_TDKS = TRUE in the $rem section. In addition, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included, and note that the JOBTYPE $rem variable is not used as it’s not applicable:

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.6, 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.

By quantizing select nuclei, NEO 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 ^{, 282} 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 ^{, 282} 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. NEO-SCF(V) calculations can be invoked by setting NEO_SCFV = TRUE:

In addition, the variables for NEO-HF 13.5.2.1 (or NEO-DFT 13.5.2.2) should be included, and JOBTYPE should be set to FREQ. Furthermore, the NEO-HF(V) method involves building the extended NEO-Hessian based on the NEO-HF Hessian and inputs from NEO-TDHF, and therefore the variables specific to LR-NEO-TDHF (TDDFT) ( 13.5.2.6) can be also be used to tune the calculation.

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 ^{, 1006} 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 ^{, 1007} Pavošević F., Hammes-Schiffer S.
J. Chem. Phys.
(2019), 151, pp. 074104. Link ^{, 1008} Pavošević F., Rousseau B. J. G., Hammes-Schiffer S.
J. Phys. Chem. Lett.
(2020), 11, pp. 1578. Link ^{, 1009} 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