Enhancing the applicability of multicomponent time-dependent density functional theory

Time-dependent density functional theory (TDDFT) is of interest in computational chemistry due to its ability to calculate excited state properties.^1–6 The conventional electronic formulation of linear response TDDFT in terms of the density matrix response to a small perturbative potential can be extended to multicomponent systems, where more than one type of particle is treated quantum mechanically.^7,8 Recently, such an extension was developed within the nuclear–electronic orbital (NEO)^9–13 framework, denoted NEO-TDDFT.⁸ In contrast to conventional electronic TDDFT, NEO-TDDFT allows the calculation of both electronic and vibrational excitations by treating all electrons and select nuclei, typically protons, quantum mechanically. This ability to compute electronic and proton vibrational excitations in a single calculation has implications for a wide range of applications, such as photoinduced proton-coupled electron transfer (PCET) reactions.¹⁴ NEO-TDDFT has also been shown to be an important component in a recently developed strategy for the calculation of molecular vibrational frequencies that incorporate significant anharmonic effects.¹⁵ 

This communication builds upon our previous research on NEO-TDDFT in several key aspects. In the initial NEO-TDDFT paper,⁸ the errors in the proton vibrational frequencies were ∼100–250 cm⁻¹, which is not satisfactory for vibrational analysis. In the present work, the electronic and nuclear basis sets are increased to achieve significant improvements in the quantitative accuracy. Moreover, the previous work addressed only systems with closed-shell electronic subsystems, whereas the present work extends the approach to open-shell electronic subsystems, which are directly relevant to PCET reactions. In addition, the present work presents not only excitation energies but also nuclear transition densities, which are important for characterizing the proton vibrational excitations used in the recently developed molecular vibrational analysis scheme.¹⁵ All of these capabilities will be essential for applications to a wide range of chemical and biological systems.

In the NEO framework, both the electronic and protonic density matrices respond to small perturbations in their respective potentials, and the description of this linear-response process leads to NEO-TDDFT. The detailed derivation of NEO-TDDFT is presented in Ref. 8. Here, we simply present the working equation that is solved for the electronic and protonic excitation energies

where the matrix elements in Eq. (1) are

Here, δP denotes the change in the density matrix and ε denotes an orbital energy. The superscripts e and p denote electrons and protons, respectively. The lower case indices i and j denote occupied electronic orbitals, while the indices a and b denote virtual electronic orbitals. The upper case indices are defined analogously for protonic orbitals. Furthermore, the electronic exchange-correlation energy functional, the protonic exchange-correlation energy functional, and the electron-proton correlation energy functional are given by E_exc, E_pxc, and E_epc, respectively. The solution of Eq. (1) produces the excitation energies ω and is valid for both closed-shell and open-shell systems, although the technical implementation of the open-shell equations is slightly more complex.

In the current formulation of linear response NEO-TDDFT, which invokes the adiabatic approximation, only single excitations can be captured, and, in principle, these excitations could be of electronic, protonic, or mixed electron-proton vibronic character.⁸ However, for electronically adiabiatic systems, typically the excitations are either electronically or protonically dominated and thus can be described as pure electronic or vibrational excitations to a reasonable approximation. The character of the excitation can be evaluated by examination of the corresponding eigenvector: electronic excitations are dominated by X^e, and protonic excitations are dominated by X^p. The eigenvectors are subject to the orthonormalization condition

Equations (2)–(6) contain functional derivative terms stemming from the electron-proton correlation functional $Eepcρe,ρp$ that arises due to the multicomponent nature of the system. To date, several electron-proton correlation functionals have been developed within the NEO framework.^10–13 The epc17-2 functional,^11,12 which has been shown to provide accurate proton affinities and reasonable proton densities, is used in this communication, and the results are compared to those obtained without any electron-proton correlation functional, denoted no-epc.

Previously, the efficacy of NEO-TDDFT for calculating accurate proton vibrational excitation energies was tested on the molecular systems FHF⁻ and HCN.⁸ Here, we present additional calculations on these systems utilizing larger electronic and protonic basis sets, thereby significantly improving the quantitative accuracy of the NEO-TDDFT method. Additionally, two electronically open-shell systems, HOOH⁺ and HNF₂⁺, are studied to illustrate the applicability of this method to open-shell systems. This work focuses on the proton vibrational excitation energies because, as demonstrated previously, the electronic excitation energies calculated with NEO-TDDFT are nearly identical to those calculated with conventional electronic TDDFT.⁸ For benchmarking purposes, the proton vibrational excitation energies are compared to those obtained using the Fourier grid Hamiltonian (FGH) method^16,17 at the DFT/B3LYP/def2-QZVP^18–20 level with a three-dimensional grid spanning the relevant region for the density of the quantum proton. These grid-based calculations are numerically exact for these types of electronically adiabatic systems.

The geometries of all molecules were optimized using conventional electronic structure theory with the B3LYP electronic exchange-correlation functional^18,19 and the def2-QZVP basis set.²⁰ All electrons and one proton were treated quantum mechanically on the same level in the NEO-TDDFT calculations, and all nuclei except the quantum proton were fixed to their optimized values for the NEO-TDDFT and FGH calculations. For the NEO-TDDFT calculations, the cc-pVDZ or cc-pVTZ electronic basis set²¹ was used for heavy atoms, and the cc-pV5Z or cc-pV6Z electronic basis set^21,22 was used for the quantum hydrogen. The previously used⁸ even-tempered 8s8p8d proton basis set with exponents spanning the range from 2$2$ to 32 was expanded to include f type functions using the same range of exponents, with an even tempered 8s8p8d8f basis set yielding converged results for the excitations of interest. All calculations presented herein were performed using a developer version of the GAMESS program.²³ 

The results of the NEO-TDDFT calculations compared to the grid-based reference for the four molecules of interest are given in Table I. The results obtained with NEO-TDDFT/epc17-2 using the 8s8p8d8f/cc-pV6Z basis sets are within 30 cm⁻¹ of the grid-based reference values for most excitation energies. The exceptions are deviations of ∼100–150 cm⁻¹ observed for the FHF⁻ stretch and the HNF₂⁺ asymmetric bend. Note that these errors are significantly lower than those obtained from the initial presentation of the NEO-TDDFT method.⁸ The agreement of the HOOH⁺ and HNF₂⁺ results with the grid reference also demonstrates the capability of NEO-TDDFT to capture proton vibrational excitations in open-shell systems, which has not been investigated previously.

In addition to Table I, a more detailed analysis of the dependence of the NEO-TDDFT proton excitation energies on the electronic and nuclear basis sets for the quantum proton is provided in the supplementary material. A comparison of the results obtained with the 8s8p8d8f/cc-pV6Z basis sets to the results obtained with the 8s8p8d/cc-pV5Z basis sets (Table I and Table S1) demonstrates an improvement in excitation energies by ∼100–300 cm⁻¹ for most modes when f proton basis functions and more electronic basis functions on the quantum proton are included. Greater flexibility in the quantum proton basis sets may be required to describe excited proton vibrational states because they contain nodes and are more delocalized than the ground state. The results obtained with the 6s6p6d6f/cc-pV6Z nuclear/electronic basis sets, for which the two most diffuse nuclear basis functions were eliminated, were found to be virtually identical to the results obtained with the 8s8p8d8f/cc-pV6Z nuclear/electronic basis sets (Tables S1 and S2). A comparison of the results obtained with the 6s6p6d6f/cc-pV6Z and 6s6p6d6f6g/cc-pV6Z nuclear/electronic basis sets (Table S2) indicates that the addition of 6g shells has negligible impact on the proton vibrational excitation energies, with differences of only ∼5–15 cm⁻¹ observed. Note that a larger electronic basis set on hydrogen is required for NEO calculations compared to grid-based calculations because a single basis function center for hydrogen is used in NEO calculations, whereas the electronic basis functions are located at each grid point for grid-based calculations.

The NEO-TDDFT results obtained without an electron-proton correlation functional (no-epc) using the 8s8p8d8f/cc-pV6Z basis sets are similar to those obtained with the epc17-2 functional, where the largest differences of ∼60–70 cm⁻¹ are observed for the stretch modes of FHF⁻, HOOH⁺, and HNF₂⁺. This similarity indicates that although electron-proton correlation is important for ground state properties such as accurate nuclear densities¹¹ and proton affinities,¹² it is not as essential for the calculation of accurate proton vibrational excitation energies for the systems studied. However, the electron-proton correlation functional may be important for systems with strong nonadiabatic couplings between electronic and/or vibrational states.

For a finite basis set, the values of the NEO-TDDFT excitation energies also depend on the positions of the nuclear/electronic basis function centers for the quantum protons. For ground state NEO-DFT calculations, the basis function center positions associated with the quantum protons can be optimized variationally, but for excited state NEO-TDDFT calculations, these positions are not as well-defined, although the dependence of the excitation energies on these positions is expected to be negligible for a sufficiently large basis set. In principle, the nuclear/electronic basis function centers could be varied so that the excitation energy is minimized for a particular excited state. However, such a state-dependent procedure is expensive and not theoretically justified because the variational theorem does not apply to excitation energies computed with TDDFT. Consequently, a well-defined protocol is required for the positioning of the nuclear/electronic basis function centers in NEO-TDDFT calculations. Two physically motivated, well-defined protocols are as follows: (1) position the basis function centers at the XH bond distance obtained from a conventional electronic DFT geometry optimization or (2) position the basis function centers at the expectation value of the nuclear density obtained from a NEO-DFT calculation. In Table I, all calculations were performed with the nuclear/electronic basis functions centered at the XH bond distance optimized with conventional electronic DFT for each molecule. The two protocols for positioning the quantum proton basis function centers lead to a difference of ∼10–20 cm⁻¹ for most proton vibrational modes (Table S3). This difference is small, and when comparing across molecules and even modes within a given molecule, there is no clear distinction as to which protocol is superior. Moreover, the dependence of the NEO-TDDFT excitation energies on the basis function center positions is only observed for a finite basis set and would vanish for a complete basis set.

Calculations of the proton vibrational excitation energies were also performed using the Tamm-Dancoff approximation within the NEO framework, denoted NEO-TDDFT-TDA, given by

where A^e, A^p, and C are defined in Eqs. (2), (4), and (6), respectively. Our previous work⁸ illustrated that the proton vibrational excitation energies computed with NEO-TDDFT-TDA are much less accurate than those computed with NEO-TDDFT. As shown in Table S4, the NEO-TDDFT-TDA results are still much less accurate, with errors of ∼2000 cm⁻¹, even with the larger nuclear and electronic basis sets used here. The magnitude of the differences in the proton vibrational excitation energies computed with NEO-TDDFT-TDA compared to those computed with NEO-TDDFT is significantly larger than the typical differences exhibited between electronic excitation energies computed with conventional electronic TDDFT-TDA compared to those computed with TDDFT.²⁴ 

Analysis of the various contributions to the working equations elucidates the main reason that NEO-TDDFT-TDA is not as reasonable an approximation for proton vibrational excitation energies as conventional electronic TDDFT-TDA is for electronic excitation energies. In conventional electronic TDDFT, the electronic orbital energy differences appearing in the diagonal part of the A matrix dominate the entire equation. These orbital energy differences are reasonable zeroth-order approximations to the excitation energies, and the kernel parts in both A and B simply act as perturbative corrections. In NEO-TDDFT, however, the proton orbital energy differences are far from being reasonable approximations to the proton vibrational excitation energies. As a result, the kernel parts in both A and B cannot be treated as small perturbative corrections in NEO-TDDFT, and they need to be accurately included to obtain even qualitatively accurate proton vibrational excitation energies. Specifically, the kernel parts of A^p have similar magnitudes to those of the proton orbital energy differences for the systems studied. A more detailed analysis of the different contributions is complicated because the proton vibrational excitations are not dominated by a single transition but rather are linear combinations of many transitions. Furthermore, we have evidence that both the B^e and C blocks are essential for obtaining accurate proton vibrational excitation energies. This evidence is based on NEO-TDDFT-TDA calculations with no electron-proton correlation functional, where the B^p block is identically zero and the results are comparable to those obtained with the epc17-2 functional. The qualitative inaccuracy of the proton vibrational excitation energies obtained using NEO-TDDFT-TDA with and without the epc17-2 functional implies that the B^e and the C blocks are important for attaining the accuracy of NEO-TDDFT.

In addition to calculating the proton vibrational excitation energies, we implemented the methodology to visualize the nuclear transition densities associated with these excitations. If an excitation is dominated by protonic transitions, then the protonic transition density can be constructed according to

where the indices I and A denote protonic occupied and virtual orbitals, respectively.⁸ The terms $XIAp$ and $YIAp$ are elements of the NEO-TDDFT eigenvector associated with a given excitation. Protonic transition densities can be used to characterize the NEO-TDDFT excitations and link them descriptively to their normal mode counterparts. As an example, protonic transition densities for HNF₂⁺ and FHF⁻ are plotted in Figs. 1 and 2, respectively. Note the p-type shape of the transition density for each excitation and its spatial orientation, which corresponds to the analogous mode for a given vibrational excitation with all heavy nuclei fixed.

The proton vibrational excitation energies calculated with NEO-TDDFT, as well as the grid-based references, are determined for fixed heavy nuclei. As a result, they do not correspond to physically meaningful molecular vibrations, which are expected to include combined motions of the quantum proton(s) and the heavy nuclei. Recently, our group has developed an approach that enables the calculation of the vibrational frequencies of an entire molecule within the NEO framework.¹⁵ This strategy entails constructing and diagonalizing an extended NEO Hessian matrix that depends on the expectation values of the quantum nuclei as well as the coordinates of the classical nuclei. The elements of this extended Hessian matrix are computed with substantial input from NEO-TDDFT, and the accuracy of the resulting molecular vibrational frequencies depends strongly on the accuracy of the underlying NEO-TDDFT calculations. Thus, the benchmarking and extensions of NEO-TDDFT presented herein are critical for the calculation of accurate molecular vibrational frequencies within the NEO framework. Note that this formulation of NEO-TDDFT is not expected to describe double excitations, and investigation of multicomponent methods for describing double excitations is an important direction for the future.

In this communication, we have investigated the electronic and nuclear basis set convergence for proton vibrational excitation energies computed with NEO-TDDFT and have extended this method to allow applications to open-shell systems as well as the calculation of nuclear transition densities. The expansion of the nuclear basis to include f protonic basis functions as well as a larger electronic basis set associated with the quantum hydrogen produces proton vibrational excitation energies that are mostly within ∼30 cm⁻¹ of the grid reference for the molecules studied. Additionally, NEO-TDDFT was shown to be effective in calculating proton vibrational excitation energies for open-shell systems. Finally, an approach for computing the nuclear transition densities associated with the proton vibrational excitations was implemented and used to visualize these excitations. These nuclear transition densities are useful for characterizing the proton vibrational excitations and determining the spatial orientation of the corresponding vibrational mode. These capabilities will be essential for the calculation of molecular vibrational frequencies for complex molecules including multiple quantum protons within the NEO framework. An advantage of this general approach is that nuclear quantum effects such as proton delocalization, anharmonicity, and zero point energy are included during the geometry optimizations, and typically the most significant anharmonic effects are incorporated into the molecular vibrational frequency calculations.

See supplementary material for more data on the impact of electronic and nuclear basis sets on proton vibrational excitation energies.

This work was supported primarily by the National Science Foundation Grant No. CHE-1762018. Work by Y.Y. was partially supported as part of the Center for Light Energy Activated Redox Processes (LEAP), an Energy Frontier Research Center funded by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences (BES), under Award No. DE-SC0001059. Work by F.P. was supported by the U.S. Department of Energy, Office of Science, Offices of Basic Energy Sciences and Advanced Scientific Computing Research, Scientific Discovery through the Advanced Computing (SciDAC) program.