Alternative forms and transferability of electron-proton correlation functionals in nuclear-electronic orbital density functional theory

An accurate inclusion of nuclear quantum effects in quantum chemistry calculations is essential for a meaningful description of a wide variety of chemical phenomena.^1,2 One promising direction for including nuclear quantum effects in quantum chemistry calculations is multicomponent self-consistent-field (SCF) methods.^3–5 In these multicomponent SCF methods, more than one type of particle is treated quantum mechanically, typically electrons and nuclei, using extensions of single-component methods developed for conventional electronic structure calculations. This study focuses on the nuclear-electronic orbital (NEO) method, a specific type of multicomponent SCF approach designed to study chemical systems.⁴ Although some multicomponent methods treat all nuclei quantum mechanically,^6–8 the NEO approach treats only specific protons quantum mechanically to avoid complications involving translations and rotations. The quantum protons are selected based on their importance in describing the chemical phenomena of interest.

Historically multicomponent SCF methods have provided a poor description of nuclear quantum effects, producing highly over localized proton densities and qualitatively inaccurate relative energies, thereby adversely affecting all other properties related to the nuclear quantum effects. A large contributing factor to this poor description is the opposite nature of the electron-proton and electron-electron Coulombic interactions. As electrons and protons attract each other, a correct description of electron-proton correlation is important for even a qualitatively correct description of the system. Owing to the importance of electron-proton correlation when both the electrons and protons are treated quantum mechanically, efforts have been made to include explicit electron-proton correlation directly in the wave function in NEO methods using the explicitly correlated Hartree-Fock^9–11 (XCHF) and reduced explicitly correlated Hartree-Fock^12–14 (RXCHF) approaches. While these XCHF and RXCHF approaches have proven to be more accurate for calculating nuclear densities than methods based on a mean field NEO Hartree-Fock wave function, they still do not provide sufficiently accurate three-dimensional proton densities and are extremely computationally demanding with a formal scaling of up to N¹² with respect to system size.

An especially promising direction within the NEO framework focuses on multicomponent density functional theory (DFT).^6,7,15 A multicomponent analog of the Hohenberg-Kohn theorem¹⁶ states that the total energy of a multicomponent system is a functional of the one-particle densities of each component,^6,17 leading to a multicomponent version of the Kohn-Sham equations. The only new term without a known form arising in the NEO-DFT approach is the electron-proton correlation functional. Similar to the electronic exchange-correlation functional in conventional electronic DFT,^18,19 the exact electron-proton correlation functional is unknown and therefore must be approximated.

Initial NEO-DFT electron-proton correlation functional development derived functionals starting from the explicitly correlated wave function ansatz underlying the NEO-XCHF method.^15,20–22 These functionals suffered from the same problems that plague the NEO-XCHF approaches in terms of an inaccurate description of the three-dimensional proton densities and a high computational scaling factor. More recently, an electron-proton correlation functional was derived from a multicomponent version of the Colle-Salvetti formulation of the correlation energy.^23–25 The Colle-Salvetti formulation of the electronic correlation energy²⁶ provides the starting point of the Lee-Yang-Parr (LYP) correlation functional in single-component electronic DFT,²⁷ and the new electron-proton correlation functional can be thought of as a multicomponent analog of the LYP functional. This electron-proton correlation functional, denoted epc17, has been shown to be able to accurately predict three-dimensional nuclear densities and proton affinities.^23,24 Additionally, the epc17 functional has the same formal scaling as single-component DFT and thus offers the potential for including the nuclear quantum effects of select nuclei for any system in which single-component DFT can be performed.

The derivation of the epc17 functional requires an approximation for a functional parameter β, which is interpreted as specifying the length scale over which electron-proton correlation effects play a significant role. An essential realization in the derivation of the epc17 functional was that the approximation for β should depend on both the electron and proton densities. This approximation for β differs from the LYP functional²⁷ and other previous attempts^28,29 at constructing an electron-proton correlation functional using the Colle-Salvetti formulation of the correlation energy, which all assumed that β depends on only the electron density. For the epc17 functional, β was approximated to be proportional to the geometric mean of the inverse Wigner-Seitz radii of the electron and proton. In this study, we explore alternative approximations for β and thereby develop the alternative epc18 functional, which provides equally accurate results as the epc17 functional. In addition, the previous NEO-DFT/epc17 calculations were performed in conjunction with only a single electronic exchange-correlation functional, namely, the B3LYP functional.^27,30,31 Herein we investigate the transferability of these types of NEO-DFT electron-proton correlation functionals to be used in conjunction with a wide range of other electronic exchange-correlation functionals. Our calculations confirm this transferability, which is critical for the general applicability of the NEO-DFT approach to a wide range of chemical systems.

In NEO-DFT, the total energy can be written in terms of the one-particle electron and proton densities, $ρe$ and $ρp$⁠, respectively, as

In this equation, E_ext is the interaction of the densities with the classical nuclei, E_ref contains both the non-interacting kinetic energies and the classical electrostatic interactions of the one-particle densities, E_exc is the electronic exchange-correlation energy, E_pxc is the proton exchange-correlation energy, and E_epc is the electron-proton correlation energy. The partitioning of the energy in this form has many advantages. The first two terms are similar to terms that arise in single-component electronic DFT and have known forms. For the third term, $Eexcρe$⁠, existing electronic exchange-correlation functionals can be used. For the fourth term, $Epxcρp$⁠, the quantum nuclei are sufficiently localized in molecular systems such that this term can be approximated by the diagonal exchange terms that cancel out the spurious self-interaction error in E_ref. Therefore, E_epc is the only new term without a known form arising in the NEO-DFT formalism. Similar to single-component electronic DFT, the proof of Eq. (1) does not provide the exact form of the E_epc functional.^6,17 As such, most recent NEO-DFT development has centered on new approximations to E_epc.

Applying the usual Kohn-Sham prescription to Eq. (1) leads to the following set of coupled equations in atomic units:

with

In these equations, unprimed (primed) indices represent electrons (quantum protons), the superscripts e (p) indicate electronic (proton) quantities, and m_p is the mass of the quantum proton. Moreover, N_e is the number of electrons, $r1e$ is the spatial coordinate of electron 1, and $ψie$ is the spatial orbital of electron i. The proton quantities are defined analogously. The total single particle densities are defined as

For the electronic effective potential $veffe$⁠, $vexte$ is the electronic external potential, $vJeee$ is the electron-electron Coulombic potential, $vxce$ is the electronic exchange-correlation potential, $vJpee$ is the electron-proton Coulombic potential, and $vepce$ is the electron-proton correlation potential for the electrons. The terms within the proton effective potential $veffp$ are defined analogously. The individual terms in Eq. (3) arise from the functional derivative of Eq. (1) with respect to either the electron or proton density. These equations are specific to a closed-shell electronic subsystem, where each electronic orbital is doubly occupied, and an open-shell proton subsystem, where each quantum proton occupies a different nuclear orbital.

As mentioned in the Introduction, recently a new type of electron-proton correlation functional, denoted epc17,²⁴ was derived and implemented. The epc17 functional is the first NEO-DFT functional capable of calculating accurate proton densities and proton affinities.^23,24 In the following, we briefly outline the derivation and the form of the epc17 functional. Additional details associated with the derivation of the epc17 functional are provided in Ref. 24.

The epc17 functional is derived using a multicomponent extension of the Colle-Salvetti²⁶ formulation for correlation energy, starting from a multicomponent wave function ansatz,

where $ΨFCIe$ is a full configuration interaction (CI) electronic wave function that contains all electron-electron correlation and $xe$ represents the spin coordinates of all electrons. Again, the corresponding proton quantities are defined analogously. Moreover, R is the center of mass for an electron-proton pair, which is effectively the proton coordinate, and r is the electron-proton distance for this pair. Here β and ξ are unspecified functions that are later assumed to have specific forms. The Colle-Salvetti prescription leads to an expression for the electron-proton correlation energy as

where $P2,epFCI$ is the second-order reduced electron-proton density matrix and $P2,epFCIR,0=ρe(R)ρp(R)$⁠.

The Colle-Salvetti formulation also enables the derivation of an approximate relationship between β and ξ.²⁴ This relationship, in combination with neglecting the third and fourth terms of Eq. (6), gives an approximate expression for the electron-proton correlation energy that depends only on β,

In order to define an electron-proton correlation functional, this expression for the electron-proton correlation energy must be defined in terms of the electron and proton densities. Similar to the derivation of the LYP electronic correlation functional,²⁷ which used a single-component formalism of the Colle-Salvetti electronic correlation energy as its starting point, we have interpreted β as the length scale over which correlation contributes significantly. For the LYP functional, β was approximated to be proportional to the inverse of the electronic Wigner-Seitz radius, r_s,e, corresponding to a uniform electron gas model,

The correlation length scale for electron-proton correlation should depend on both the electron and proton densities. Therefore, in the derivation of the epc17 functional, β was approximated to be proportional to the geometric mean of the inverse Wigner-Seitz radii corresponding to a uniform electron or proton gas model,

With this definition of β, the final form of the epc17 functional is

where a, b, and c are adjustable parameters. As the epc17 functional depends on only the electron and proton densities, in a technical sense, the epc17 functional may be viewed as a multicomponent local density approximation (LDA) functional. However, it should be noted that the epc17 functional is not derived from a multicomponent homogeneous gas reference system. Two sets of parameters for the epc17 functional were identified: the parameters used in epc17-1 are more accurate for calculating proton densities, while the parameters used in epc17-2 are more accurate for calculating energies. Note that the electronic correlation functionals based on the Colle-Salvetti formulation of the correlation energy include the terms in Eq. (6) that depend on the gradient of the density. Although Eq. (6) depends on the Laplacian of the density, the terms depending on the Laplacian of the density can be converted to terms depending on the gradient of the density through an integration-by-parts. It is expected that more accurate electron-proton correlation functionals would be obtained by including these gradient terms, possibly allowing for the accurate calculation of proton densities and energies with a single set of parameters. Research exploring this possibility is currently ongoing.

The dependence of β on both the electron and proton densities is a key difference between the epc17 functional and other electron-proton correlation functionals developed using an ansatz based on the Colle-Salvetti formulation,^28,29 which defined β to depend on only the electron density. While it appears that an electron-proton correlation functional must define β to depend on both the electron and proton densities, the definition of β in Eq. (9) is not unique. Two possible alternative definitions of β are that it is proportional to the arithmetic mean of the inverse Wigner-Seitz radii or proportional to the inverse of the arithmetic mean of the Wigner-Seitz radii of the electrons and protons,

Replacing the proportionalities with equalities, these two new definitions of β are substituted into Eq. (7), leading to two new LDA-type electron-proton correlation functionals with adjustable parameters a, b, and c. Note that when using the geometric mean, as in the epc17 functional, the inverse of the mean of the Wigner-Seitz radii and the mean of the inverse Wigner-Seitz radii of the electrons and protons are identical.

In partitioning the total energy of the system in Eq. (1), E_exc, E_pxc, and E_epc are defined to be independent.³² Previous work²¹ demonstrated that the electronic exchange-correlation and electron-proton correlation functionals are effectively independent for electron-proton correlation functionals derived from the NEO-XCHF wave function ansatz. In other words, these electron-proton correlation functionals are transferable for use in conjunction with different electronic exchange-correlation functionals. As discussed above, however, those previously developed electron-proton correlation functionals are not sufficiently accurate or computationally practical for studying molecular systems greater than diatomics with relatively small basis sets. In this study, we demonstrate the independence and transferability for the electron-proton correlation functionals derived from the multicomponent version of the Colle-Salvetti formulation, namely, the epc17 and epc18 functionals derived in Secs. II A and II B.

The accuracy of the electron-proton correlation functionals derived with the new forms for the correlation length β was tested for computing proton densities and proton affinities. The proton densities were analyzed by performing NEO-DFT calculations on HCN and FHF⁻ with the proton and all electrons treated quantum mechanically. The proton affinities were analyzed by performing NEO-DFT calculations on a test set of 23 molecules identical to that used in a previous study with the epc17-2 functional.²³ Both the proton density and the proton affinity calculations used the def2-QZVP electronic basis set³³ and the B3LYP electronic exchange-correlation functional,^27,30,31 as well as a series of other electronic exchange-correlation functionals for the proton affinity calculations. For consistency with previous work, the proton density calculations used an 8s8p8d even-tempered nuclear basis set with exponents ranging from $22$ to 32, and the proton affinity calculations used a 10s10p10d even-tempered nuclear basis set with exponents ranging from 2 to $322$⁠. For all of these calculations, the positions of the classical nuclei were obtained from a geometry optimization using conventional electronic DFT, and the centers of the nuclear and electronic basis functions associated with the quantum proton were placed at the position of the classical proton determined from this conventional electronic DFT geometry optimization.

The NEO-DFT/epc18 proton densities for HCN and FHF⁻ were compared to the densities previously computed with the NEO-DFT/epc17-1 method, the NEO-DFT/no-epc method, which does not include any electron-proton correlation,³⁴ and the Fourier Grid Hamiltonian (FGH) method.^35,36 In the FGH method, the conventional electronic DFT energy is computed for the proton positioned at each grid point of a three-dimensional grid, and subsequently the three-dimensional proton vibrational wave functions and energies are computed numerically. Although the NEO-DFT method is also expected to be applicable to nonadiabatic systems, this benchmarking is performed for electronically adiabatic systems to enable the use of the FGH method as a numerically exact reference. The parameters a, b, and c were fit using the new forms for the correlation length β in Eqs. (11) and (12) by a grid-based search in the parameter space aimed at minimizing the root-mean-square deviation (RMSD) between the calculated NEO-DFT proton density and the FGH proton density for the two molecules. In Sec. IV, we present results only for the arithmetic mean of the inverse Wigner-Seitz radii of the electron and proton, $βRarith-inv$⁠, as defined in Eq. (11), with the parameters a, b, and c determined to be 1.8, 0.1, and 0.03, respectively. This electron-proton correlation functional is denoted as epc18-1. The results for the inverse of the arithmetic mean of the Wigner-Seitz radii of the electron and proton, $βRinv-arith$⁠, as defined in Eq. (12), with the parameters a, b, and c determined to be 0.54, 0.0021, and 120.0, respectively, are presented in the supplementary material.

The NEO-DFT/epc18 proton affinities for a test set of 23 molecules²³ were compared to those computed with the NEO-DFT/epc17-2 method, the NEO-DFT/no-epc method, and experimental data. Under the ideal gas approximation and the assumption that the rotational energy of the molecule and the vibrational energy associated with the classical nuclei do not change upon protonation, the proton affinity of a chemical species A is defined as²³ 

where $EAH+$ is the energy of AH⁺ calculated with the NEO-DFT method and $EA$ is the energy of A calculated using conventional electronic DFT. The parameters a, b, and c were fit using the new forms for the correlation length β in Eqs. (11) and (12) by a grid-based search in the parameter space aimed at minimizing the RMSD between the NEO-DFT proton affinities calculated with the B3LYP electronic exchange-correlation functional for the test set of molecules and their corresponding experimental values. In Sec. IV, we present results for only the arithmetic mean of the inverse Wigner-Seitz radii of the electron and proton, $βRarith-inv$⁠, with the parameters a, b, and c determined to be 3.9, 0.5, and 0.06, respectively. This new functional is denoted as epc18-2. When β was approximated by the inverse of the arithmetic mean of the Wigner-Seitz radii of the electron and proton, $βRinv-arith$⁠, no suitable set of parameters was found that could predict even qualitatively accurate proton affinities. In addition, to demonstrate the independence of electron-proton correlation and electronic exchange-correlation energies, the proton affinities were calculated with eight different electronic exchange-correlation functionals at the NEO-DFT/no-epc, NEO-DFT/epc17-2, and NEO-DFT/epc18-2 levels.

Figures 1 and 2 show two one-dimensional slices of the proton density for FHF⁻ and HCN, respectively, calculated with the FGH, NEO-DFT/no-epc, NEO-DFT/ep17-1, and NEO-DFT/epc18-1 methods. For both FHF⁻ and HCN, the epc18-1 functional agrees well with the FGH method, which serves as the reference for these electronically adiabatic systems. Although not readily apparent from the one-dimensional slices of the proton density shown in Figs. 1 and 2, the epc18-1 functional provides a slightly more accurate proton density over three-dimensional space than the epc17-1 functional relative to the FGH method. The RMSDs in the proton density over three-dimensional space for the epc17-1 and epc18-1 functionals relative to the FGH method for the HCN molecule are 0.267 a.u. and 0.218 a.u., respectively, while for the FHF⁻ molecule these RMSDs are 0.120 a.u. and 0.072 a.u., respectively. Note that the previously developed epc17-1 functional was not optimized by quantitatively minimizing the RMSD. Nevertheless, the proton densities are of similar accuracy for the epc17-1 and epc18-1 functionals. As shown previously,²⁴ the NEO-DFT/no-epc method severely over localizes the proton density and does not provide even a qualitatively correct density.

The mean unsigned errors (MUEs) for the proton affinities of the test set of molecules are presented in Fig. 3 for eight different electronic exchange-correlation functionals in conjunction with the NEO-DFT/no-epc, NEO-DFT/epc17-2, and NEO-DFT/epc18-2 electron-proton correlation functionals. The ranges of proton affinity unsigned relative errors for the NEO-DFT/no-epc, NEO-DFT/epc17-2, and NEO-DFT/epc18-2 electron-proton correlation functionals in combination with the B3LYP electronic exchange-correlation functional are 0.69–0.91 eV, 0.09–0.22 eV, and 0.06–0.24 eV, respectively. The ranges of relative errors with the B3LYP functional are representative of the other electronic exchange-correlation functionals tested. Additionally, the errors in proton affinity for each molecule with each combination of electronic exchange-correlation and electron-proton correlation functional are provided in the supplementary material. Overall, the accuracy of the calculated proton affinities of the epc18-2 functional is similar to that of the epc17-2 functional. Note that the previously developed epc17-2 functional was determined by fitting to the experimental proton affinity of HCN rather than minimizing the MUE for the entire data set. Nevertheless, the proton affinities are of similar accuracy for the epc17-2 and epc18-2 functionals. As shown previously for the B3LYP functional,²³ a correct description of electron-proton correlation plays a critical role in obtaining accurate proton affinities using multicomponent SCF methods. Figure 3 shows that the neglect of electron-proton correlation results in large errors relative to experiment for all electronic exchange-correlation functionals.

For the epc17-2 and epc18-2 functionals, excluding the SVWN electronic exchange-correlation functional,^37,38 the error in proton affinity is of similar magnitude for all electronic exchange-correlation functionals, ranging from 0.06 eV to 0.14 eV for the epc17-2 electron-proton correlation functional and from 0.06 eV to 0.11 eV for the epc18-2 electron-proton correlation functional. These calculations demonstrate the independence of E_exc and E_epc in Eq. (1). For the SVWN functional, the MUEs for the proton affinities are significantly higher at 0.30 eV and 0.26 eV for the epc17-2 and epc18-2 electron-proton correlation functionals, respectively. Given that the MUE for the SVWN electronic exchange-correlation functional is also larger without inclusion of any electron-proton correlation (no-epc), these results indicate that the SVWN functional is a less accurate electronic exchange-correlation functional for these systems. Overall, the results in Fig. 3 show that any modern electronic exchange-correlation functional may be used in combination with the epc17 and epc18 classes of electron-proton correlation functionals.

The design of alternative electron-proton correlation functionals and the analysis of the transferability of such functionals within the NEO-DFT framework were explored. A new form of the parameter associated with the correlation length was utilized in the derivation of an electron-proton correlation functional within the Colle-Salvetti formulation. The correlation length is represented as the arithmetic mean of the inverse Wigner-Seitz radii of the electron and proton, in contrast to a previous approximation of the correlation length as the geometric mean of these inverse Wigner-Seitz radii. The resulting functional, denoted as epc18, performs similarly to the epc17 functional for the three-dimensional proton densities and the prediction of proton affinities relative to experiment. Furthermore, the transferability of the epc17 and epc18 functionals for use with a series of diverse electronic exchange-correlation functionals was demonstrated in the context of proton affinity calculations. This transferability is important for the general applicability of the NEO-DFT approach.

The development of accurate electron-proton correlation functionals within the framework of multicomponent DFT is still in its infancy. The fundamental characteristics of electron-proton correlation functionals and the requirements for the accurate prediction of various physical properties are not yet well understood. Thus, the exploration of different avenues for electron-proton correlation functional development is critical. Although the epc18 functionals developed herein perform similarly to the previously developed epc17 functionals, further studies could identify advantages of one functional over the other in certain cases. Another promising direction is the development of electron-proton correlation functionals beyond the local density approximation, where terms depending on the gradients of the densities are included. Analogous to electronic exchange-correlation functionals,⁴⁵ these functionals may enable a single parameterization to provide more accurate energies and densities simultaneously. Understanding the impact of different functional forms of the correlation length within the Colle-Salvetti formulation will be important for the design of these gradient-type functionals as well. Additionally, understanding the interplay between electron-proton and electron-electron correlation has significant implications for the general applicability of multicomponent DFT.

See supplementary material for quantum proton density of FHF⁻ and HCN calculated with a β factor as defined in Eq. (12). Proton affinity data for all combinations of electron-proton correlation functional and electronic exchange-correlation functional studied.

This material is based upon work supported by the National Science Foundation Grant No. CHE-1361293 and by the Center for Chemical Innovation of the National Science Foundation (Solar Fuels, Grant No. CHE-1305124). We also thank Dr. Yang Yang for helpful discussions.