Multicomponent density functional theory: Including the density gradient in the electron-proton correlation functional for hydrogen and deuterium

Most quantum chemistry methods such as ab initio wave function theory and density functional theory (DFT) treat nuclei as fixed point charges and invoke the Born-Oppenheimer separation between the motions of the nuclei and the electrons. Such calculations cannot directly capture nuclear quantum effects, including zero-point energy, nuclear delocalization, and vibrational excitation energies, which are important for calculating properties such as vibrationally averaged geometries, geometric isotope effects, and electron-proton vibronic states in proton coupled electron-transfer systems.^1,2 Many approaches have been developed to describe these nuclear quantum effects to varying degrees. For example, zero-point energy and vibrational frequencies are commonly calculated by diagonalizing a Hessian matrix, which is an efficient approach but does not include anharmonic effects. A variety of perturbative and grid-based methods have been developed to include anharmonic effects in vibrational frequency calculations.^3–5 However, typically these methods do not include nuclear quantum effects during the optimization of geometries, the generation of reaction paths, or the propagation of reaction dynamics. Path integral molecular dynamics methods^6,7 incorporate nuclear quantum effects but require a significant amount of conformational sampling.

An alternative approach is the nuclear-electronic orbital (NEO) method,⁸ which treats electrons and select nuclei, typically protons, quantum mechanically on the same level. In the NEO method, a combined nuclear-electronic Schrödinger equation is solved with molecular orbital techniques. At least two classical nuclei are needed to eliminate difficulties with the translational and rotational degrees of freedom. The NEO method has the advantage of including the zero-point energy and nuclear delocalization of the quantum protons during the self-consistent field (SCF) procedure. Both NEO wave function methods^8,9 and NEO density functional theory (NEO-DFT) methods^10,11 have been developed. A practical advantage of NEO-DFT is that it has the formal scaling of conventional electronic DFT, which is O(N³) or O(N⁴), depending on the type of electronic functional, and density fitting methods¹² may improve the scaling. Other multicomponent DFT methods have also been explored,^13–18 but the main obstacle has been the lack of electron-proton correlation functionals that provide even qualitatively accurate proton densities, energies, and vibrational frequencies. Electron-proton correlation is highly significant because of the attractive interaction between an electron and a proton, and the neglect or inadequate description of this correlation leads to unphysical behavior due to severely overlocalized proton densities.¹⁹ In previous work, a number of important mathematical properties of the exact universal multicomponent density functional were derived.²⁰ 

Recently, a series of electron-proton correlation functionals were developed in the framework of NEO-DFT and have been shown to provide accurate proton densities,²¹ energies,²² geometries,²² and vibrational frequencies,²³ again demonstrating the importance of electron-proton correlation. These epc17 and epc18 functionals^21,24 were based on the Colle-Salvetti formalism,²⁵ which produced the popular Lee-Yang-Parr (LYP) electronic correlation functional,²⁶ with several essential differences required for describing electron-proton correlation. The epc17 and epc18 functionals could be parameterized to provide highly accurate proton densities or energies but not both simultaneously, although the epc17-2 and epc18-2 functionals achieved a reasonable compromise. In contrast to the LYP electronic correlation functional, the epc17 and epc18 functionals are based on the local density approximation (LDA) and thus depend only on the electron and proton densities. Herein, we derive an electron-proton correlation functional, denoted epc19, that also depends on the gradients of the electron and proton densities to further enhance the accuracy. Due to its dependence on the gradient terms, this functional belongs to the set of generalized gradient approximation (GGA) functionals and is analogous to the LYP electronic correlation functional.

This paper is organized as follows: Section II outlines the formalism underlying the epc19 functional and the parameterization procedure, with the full details provided in the supplementary material. Section III presents the performance of the epc19 functional for computing proton densities and energies, proton affinities, and geometry optimizations for both hydrogen and deuterium. A summary and discussion of future directions are provided in Sec. IV.

Multicomponent DFT, in which more than one type of particle is treated quantum mechanically, has been formulated as an extension of the Hohenberg-Kohn theorems.^{11,15,16,18,20,27} Within the NEO-DFT framework, the total energy can be expressed as a functional of the electron and proton densities, ρ^e and ρ ^p, respectively,

The quantities in each parenthesis are as follows: the noninteracting kinetic energy of the electrons and protons; the interaction of the external potential with the electrons and protons; the electron-electron, proton-proton, and electron-proton mean-field Coulomb interactions; the exchange-correlation energies of the electrons and protons; and the electron-proton correlation energy. The Kohn-Sham formalism can be used in conjunction with Eq. (1) to determine the electron and proton Kohn-Sham orbitals and densities, as well as the total energy, for a given classical nuclear configuration. Within this framework, existing electronic exchange-correlation functionals can be used.²⁴ Because the quantum nuclei are sufficiently localized in molecular systems, the proton exchange-correlation functional can be approximated by the diagonal Hartree-Fock exchange terms to eliminate self-interaction error.²⁸ 

The development of the epc19 functional starts with the same analog of the Colle-Salvetti²⁵ wavefunction ansatz as used in the development of the epc17 and epc18 functionals,

where $xe$ and $xp$ represent the combined spatial and spin coordinates of all electrons and protons, respectively, and $ri,e$ and $rI,p$ represent the spatial coordinates of individual electrons and protons, respectively. The protons only adopt high spin configurations in separate orbitals to avoid unphysically large Coulomb repulsion, as is reasonable in molecular systems because each proton is spatially localized.²⁹ In Eq. (2), $ΨeFCI$ and $ΨpFCI$ are the electron and proton wavefunctions from the full configuration interaction (FCI) treatment of electron-electron exchange-correlation and proton-proton exchange-correlation, respectively, each of which includes only a mean-field Coulomb interaction between electrons and protons. Electron-proton correlation is included through the correlation factor φ(⁠$ri,e$⁠, $ri,p$⁠), which has the form

where $r=|r|$ = |$ri,e−rI,p$| and $R$ is the center of mass between the ith electron and the Ith proton. The electron-proton correlation energy is defined to be the difference between the exact total energy associated with Eq. (1) and the sum of the FCI energies corresponding to the electron and proton FCI wavefunctions with the double counting of the mean-field electron-proton Coulomb energy removed.

The form of $ξ(R)$ in Eq. (3) can be determined as a function of the inverse correlation length $β(R)$⁠, as shown previously.²¹ The inverse correlation length $β(R)$ is assumed to be inversely proportional to the electron-proton Wigner-Seitz radius $rs,ep$⁠, which is defined to be the geometric mean of the electron and proton Wigner-Seitz radii from a uniform electron or proton gas model: $rs,ep=(rs,ers,p)1/2$⁠. Using the expressions for the individual Wigner-Seitz radii, $rs,e(R)∝[ρe(R)]−1/3$ and $rs,p(R)∝[ρp(R)]−1/3$⁠, the inverse correlation length is expressed as $β(R)∝[ρe(R)]1/6[ρp(R)]1/6$⁠. This expression was used to develop the epc17 functionals.²¹ The alternative assumption that $β(R)$ is proportional to the arithmetic mean of the inverse electron and proton Wigner-Seitz radii was used to develop the epc18 functionals,²⁴ which perform similarly to the epc17 functionals. Note that the inverse of $β(R)$ can be interpreted as the length scale over which electron-proton correlation exerts a significant contribution. Thus, in contrast to prior efforts at developing electron-proton correlation functionals,^30,31 the inverse correlation length should depend on both the electron and proton densities, although the specific form is somewhat flexible.

Following additional approximations analogous to those invoked in the Colle-Salvetti formulation, the electron-proton correlation energy is expressed as

where $P2,epFCI(R,0)=ρeFCI(R)ρpFCI(R)≈ρe(R)ρp(R)$⁠. In the epc17 and epc18 functionals, the last two terms involving the second-order derivatives of the pair density were neglected, producing LDA functionals that depend only on the electron and proton densities. In the present work, we include all four of the terms in Eq. (4). After integration, we approximate the GGA functional with the following expression:

where m_p is the proton mass, and the second derivative terms can be converted to gradient terms through integration by parts. This functional has five parameters: a, b, c, d, and k. Note that setting $d=0$ recovers the functional form of epc17.

The GGA functional depends on the proton mass, in contrast to the LDA functionals developed previously. The explicit dependence of the epc19 functional on the proton mass suggests that a common set of parameters could potentially describe both hydrogen and deuterium. However, the validity of some of the approximations underlying these functionals may be different for hydrogen and deuterium and therefore could require the use of different parameters for describing hydrogen vs deuterium. Nevertheless, for simplicity, we use the same set of parameters to describe both isotopes with the epc19 functional. Although the energies and densities for deuterium are of similar accuracy as those for hydrogen, more quantitatively accurate results could potentially be obtained for deuterium with reparameterization of the epc19 functional.

Another potential complication with deuterium is that it is a boson rather than a fermion. For molecular systems, however, the nuclear orbitals are sufficiently localized, with only one deuterium per orbital in the lowest-energy state, that the energy expressions derived for hydrogen are also applicable to deuterium. Moreover, these bosonic effects impact the deuterium-deuterium exchange interactions but do not impact the validity of the epc19 functional. In this paper, these issues are not directly relevant because the molecules studied contain only a single deuterium.

We determined the epc19 parameters for hydrogen by minimizing the errors in the energies and the root-mean-square deviations (RMSDs) of the three-dimensional proton densities for HCN and FHF⁻. The reference proton vibrational ground state energies and proton densities were obtained from Fourier grid Hamiltonian (FGH)^32,33 calculations corresponding to the hydrogen moving on a three-dimensional potential energy surface. This surface was generated by moving the proton along a grid with 32 grid points per dimension and performing conventional DFT/B3LYP/def2-QZVP^26,34,35 single-point energy calculations at each grid point. The resulting proton vibrational ground state energies and proton densities are considered to be numerically accurate for these types of electronically adiabatic systems and serve as the benchmark. For the NEO-DFT calculations, the proton and all electrons were treated quantum mechanically. The geometries of HCNand FHF⁻ were optimized at the conventional DFT/B3LYP level with the def2-TZVPD³⁶ basis set, and the heavy atoms were fixed at these geometries for both the NEO calculations and the FGH calculations. All of the NEO calculations were performed with the B3LYP electronic functional and the def2-QZVP³⁵ electronic basis set, with an even-tempered 8s8p8d proton basis set.^21,37 During the parameterization, the basis functions associated with the quantum proton were centered at the equilibrium position of the hydrogen obtained through conventional DFT geometry optimization. Afterselecting the parameters for the epc19 functional, the calculations were performed again while variationally optimizing the basis function center positions to produce the energies and proton densities presented herein, with the exception of the proton affinities, which were computed with the basis function centers at the equilibrium hydrogen positions optimized with conventional DFT. All of these calculations were performed with an in-house developer version of GAMESS.³⁸ 

In addition to minimizing the errors in the energies and proton densities, we examined the shapes of the proton densities to avoid unphysical features and confirmed that the linear equilibrium geometries for HCN and FHF⁻ were retained. The epc19 functionalparameters a, b, c, d, and k determined from this procedure are 1.9, 1.3, 8.1, 1600, and 8.2, respectively. During the parameterization procedure, we found that smaller values of k tended to produce bent equilibrium geometries for HCN and/or FHF⁻. Given this relatively large value of k, which dampens the overall contribution from the gradient-dependent term in Eq. (5), a relatively large value of d was required to obtain accurate proton energies and densities. Although these parameters were determined by fitting only to data for hydrogen, this same set of parameters was used to describe deuterium as well. The parameter space is sufficiently extensive that multiple parameter sets could produce similar or even better results, but these parameters provide the desired level of accuracy for the present work.

The proton densities and energies computed with the various electron-proton correlation functionals for HCN and FHF⁻ are compared to each other and to the grid reference in Tables I and II, respectively. Table I provides the root-mean square deviations (RMSDs) for the NEO-DFT proton densities with respect to the grid-based reference densities. The epc19 functional produces proton densities reasonably comparable to those obtained with the epc17-1 functional, which was previously parameterized solely for accurate proton densities.²¹ The proton densities calculated with various electron-proton correlation functionals are depicted in Fig. 1, with one-dimensional slices corresponding to the proton density along the axis connecting the heavy atoms (on axis) or the proton density perpendicular to this axis (off axis). The neglect of electron-proton correlation, denoted NEO-DFT/no-epc, produces highly overlocalized proton densities, as also indicated by the large RMSDs for the NEO-DFT/no-epc proton densities given in Table I. The proton densities calculated with the epc17-1,²¹ epc17-2, and epc19 functionals are significantly more delocalized and agree much better with the grid-based reference, with slight variations in the shapes. Note that the proton densities computed with the epc17-1 functional are presented in Ref. 21 but are not shown in Fig. 1 for clarity. Table II provides the difference between the NEO-DFT energies and the grid-based reference energies. The epc19 functional produces energies for HCN and FHF⁻ comparable to those obtained with the epc17-2 functional, which has been widely used within the NEO-DFT framework because of its accurate energies. In contrast, the energies provided by the epc17-1 functional and without electron-proton correlation (i.e., no-epc) are physically unreasonable.

To further test the epc19 functional, we computed the proton affinities for a set of 23 molecules. The proton affinity is defined as the negative enthalpy corresponding to species A gaining a proton to produce species AH⁺. Within the NEO-DFT framework, the proton affinity is defined to be the difference between the conventional DFT energy of A and the NEO-DFT energy of AH⁺, with the addition of 5/2 RT to account for the conversion from energy to enthalpy and translational energy changes.²² NEO calculations inherently include the vibrational energy associated with the quantum protons, and the change in the vibrational energy associated with the classical nuclei is assumed to be negligible upon protonation.⁹ All of the proton affinity calculations were performed with the B3LYP electronic functional and the def2-QZVP³⁵ electronic basis set. The NEO calculations were performed with an even-tempered 10s10p10d proton basis set^22,37 centered at the equilibrium position of the proton obtained through a conventional DFT geometry optimization. Table III illustrates that the epc19 and epc17-2 functionals achieve the same level of accuracy for computing this set of proton affinities with a mean unsigned error (MUE) of 0.06 eV. Note that this MUE is within the experimental uncertainty of ∼0.09 eV, as reported by Hunter and Lias.³⁹ 

A significant advantage of the NEO approach over conventional electronic structure methods is that nuclear quantum effects such as proton delocalization and zero-point energy can be included in geometry optimizations. To investigate the accuracy of the epc19 functional for geometry optimizations, we calculated the equilibrium F–F bond length for FHF⁻ with the NEO-DFT method using the epc19 functional and compared it to the equilibrium bond length obtained with conventional electronic DFT, the grid-based reference, and NEO-DFT/epc17-2. The grid-based reference equilibrium F–F distance was determined using the FGH method at the B3LYP/def2-QZVP level of theory to compute the energy for a range of F–F distances with a step size of 0.0005 Å in the region near the minimum. The NEO-DFT equilibrium F–F distances were obtained with a step size of 0.0001 Å using the same electronic functional and basis set with the 10s10p10d proton basis set. As shown previously,²² conventional DFT geometry optimization underestimates the equilibrium F–F distance. Table IV illustrates that both NEO-DFT/epc17-2 and NEO-DFT/epc19 correctly describe the increase in the equilibrium F–F distance by ∼0.02–0.03 Å when the nuclear quantum effects of the proton are included.

Tables I and II also provide the errors in the proton densities and energies for DCN and FDF⁻, where hydrogen has been replaced by deuterium. Without parameterization for deuterium, the NEO-DFT/epc19 method produces similar errors in the energies and densities for deuterium as for hydrogen. Analogously, even though the epc17-2 functional was parameterized for hydrogen, the NEO-DFT/epc17-2 method produces similar errors in the energies and densities for both hydrogen and deuterium. However, the epc19 functional continues to provide significantly more accurate nuclear densities than does the epc17-2 functional, as illustrated by Fig. 2. The qualitative increase in the equilibrium F–F distance for FDF⁻ compared to conventional DFT is reproduced by both the epc17-2 and epc19 functionals, although they overestimate the equilibrium distance by 0.0055 Å and 0.0164 Å, respectively. Moreover, both functionals qualitatively reproduce the geometric isotope effect, where the equilibrium F–F distance is smaller for FDF⁻ than for FHF⁻, although the quantitative differences are underestimated. These results suggest that reparameterizing for deuterium may be important for capturing subtle quantitative geometric isotope effects.

In this paper, we derived and parameterized the GGA version of the epc17 functional within the NEO-DFT framework by including the gradient terms with respect to both the electron and proton densities. The resulting epc19 functional provides comparable energies to the epc17-2 functional and significantly more accurate proton densities, achieving the objective of producing accurate energies and densities simultaneously. In addition, the epc19 functional achieves a similar level of accuracy as the epc17-2 functional for the proton affinities associated with a set of 23 molecules and reproduces the increased equilibrium F–F distance for FHF⁻ upon inclusion of nuclear quantum effects. These combined results demonstrate an overall improvement over the previous LDA-type electron-proton correlation functionals when the gradient terms are included. However, we emphasize that the parameterization depends on the weighting of these various properties, and a different weighting could improve any one of these properties, possibly at the expense of another property. The objective of this paper is to provide the form of a GGA functional that may be parameterized to optimize the properties of interest.

These functionals were also applied to deuterated systems to enable the investigation of hydrogen/deuterium isotope effects. Without further parameterization, the epc17-2 and epc19 functionals each produced energies and densities for deuterium that are of similar accuracy as those computed for hydrogen with the same functional. As observed for hydrogen, the nuclear densities for deuterium are significantly more accurate when computed with the epc19 functional than with the epc17-2 functional. Moreover, both functionals qualitatively reproduced the geometric isotope effect for FHF⁻ vs FDF⁻, where the equilibrium F–F distance is greater for hydrogen than for deuterium, although quantitative accuracy may require separate parameters for deuterium.

This work illustrates that the form of the epc19 functional is able to capture the essential aspects of electron-proton correlation. However, more systematic parameterization procedures with larger datasets will be necessary to develop more quantitatively accurate functionals for both hydrogen and deuterium. Future studies will focus on applying the NEO-DFT method to systems for which nuclear quantum effects such as zero-point energy and nuclear delocalization are important. The use of this method to investigate geometric isotope effects, as well as equilibrium and kinetic isotope effects, is another promising direction.

See the supplementary material for detailed formalism underlying the epc19 functional.

The authors thank Tanner Culpitt, Patrick Schneider, Dr. Fabijan Pavošević, and Benjamin Rousseau for helpful discussions. This work was supported by the National Science Foundation (Grant No. CHE-1762018).