Development of a Practical Multicomponent Density Functional for Electron-Proton Correlation to Produce Accurate Proton Densities

TL;DR

This paper introduces epc17, a new electron-proton correlation functional for multicomponent DFT, which accurately predicts proton densities and is efficiently implemented within the NEO framework for diverse applications.

Contribution

A novel electron-proton correlation functional, epc17, derived from Colle-Salvetti formalism, enhancing proton density accuracy in multicomponent DFT.

Findings

epc17 produces accurate proton densities

Efficient implementation within NEO-DFT framework

Promising for diverse quantum chemical applications

Abstract

Multicomponent density functional theory (DFT) enables the consistent quantum mechanical treatment of both electrons and protons. A major challenge has been the design of electron-proton correlation functionals that produce even qualitatively accurate proton densities. Herein an electron-proton correlation functional, epc17, is derived analogously to the Colle-Salvetti formalism for electron correlation and is implemented within the nuclear-electronic orbital (NEO) framework. The NEO-DFT/epc17 method produces accurate proton densities efficiently and is promising for diverse applications.

Equations

$$E[{\rho }_e,{\rho }_p]=$$

$$+(J_{ee}[{\rho }_e]+J_{pp}[{\rho }_p]+J_{ep}[{\rho }_e,{\rho }_p])$$

$$+(E_e^{xc}[{\rho }_e]+E_p^{xc}[{\rho }_p])+(E^{epc}[{\rho }_e,{\rho }_p]),$$

$$\displaystyle\bigg{(}-\frac{1}{2}\nabla^{2}_{e}+v_{e}(\mathbf{r}_{e})\bigg{)}\phi_{e}^{i}(\mathbf{r}_{e})=$$

$$\displaystyle\bigg{(}-\frac{1}{2m_{p}}\nabla^{2}_{p}+v_{p}(\mathbf{r}_{p})\bigg{)}\phi_{p}^{I}(\mathbf{r}_{p})=$$

$$v_e({\mathbf{r}}_e)=$$

$$+v_e^{xc}({\mathbf{r}}_e)+v_e^{epc}({\mathbf{r}}_e),$$

$$v_p({\mathbf{r}}_p)=$$

$$+v_p^{xc}({\mathbf{r}}_p)+v_p^{epc}({\mathbf{r}}_p).$$

$$\Psi ({\mathbf{x}}_{1,e},{\mathbf{x}}_{2,e},\dots {\mathbf{x}}_{N_e,e},{\mathbf{r}}_{1,p},{\mathbf{r}}_{2,p},\dots {\mathbf{r}}_{N_p,p})$$

$$=$$

$$\cdot \prod _{i\in N_e,I\in N_p}[1-\phi ({\mathbf{r}}_{i,e},{\mathbf{r}}_{I,p})],$$

$$\phi ({\mathbf{r}}_{i,e},{\mathbf{r}}_{I,p})=e^{-{\beta }^2(\mathbf{R})r^2}[1-\xi (\mathbf{R})(1-r)],$$

$$P_{2,ep}({\mathbf{x}}_{1,e}^{\prime },{\mathbf{r}}_{1,p}^{\prime };$$

$$\cdot [1-\phi ({\mathbf{r}}_{1,e}^{\prime },{\mathbf{r}}_{1,p}^{\prime })][1-\phi ({\mathbf{r}}_{1,e},{\mathbf{r}}_{1,p})],$$

$$E^{epc}=$$

$$=$$

$$-P_{2,e}^{FCI})]+\mathrm{Tr}[(-\frac{1}{2m_p}{\nabla }_p^2+v_p^{ext})(P_{1,p}-P_{1,p}^{FCI})]$$

$$+\mathrm{Tr}[V_{pp}(P_{2,p}-P_{2,p}^{FCI})]-\int \mathrm{d}{\mathbf{x}}_{1,e}\mathrm{d}{\mathbf{r}}_{1,p}\frac{1}{r}$$

$$\cdot [P_{2,ep}({\mathbf{x}}_{1,e},{\mathbf{r}}_{1,p})-P_{2,ep}^{FCI}({\mathbf{x}}_{1,e},{\mathbf{r}}_{1,p})].$$

$$E^{epc}=$$

$$\cdot [{\phi }^2({\mathbf{r}}_{1,e},{\mathbf{r}}_{1,p})-2\phi ({\mathbf{r}}_{1,e},{\mathbf{r}}_{1,p})].$$

$$\xi (\mathbf{R})\approx \frac{-9.2{\beta }^2(\mathbf{R})+14.1\beta (\mathbf{R})\sqrt{1.0{\beta }^2(\mathbf{R})-2.0\beta (\mathbf{R})+1.2}+12.0\beta (\mathbf{R})}{5.0{\beta }^2(\mathbf{R})-8.0\beta (\mathbf{R})+3.8}.$$

$$\beta (\mathbf{R})\propto \frac{1}{r_{s,ep}(\mathbf{R})}\propto {\rho }_e^{1/6}(\mathbf{R}){\rho }_p^{1/6}(\mathbf{R}).$$

$$\int \mathrm{d}\mathbf{r}$$

$$+\frac{2\pi }{3}({\nabla }_r^{2}F(\mathbf{r})){|}_{r=0}\int r^3{e}^{-{\beta }^2{r}^2}\mathrm{d}r+O(r^6),$$

$$E^{epc}\approx$$

$$+8\pi \int \mathrm{d}\mathbf{R}{P}_{2,ep}^{FCI}(\mathbf{R},0)\int \mathrm{d}r{e}^{-{\beta }^2(\mathbf{R})r^2}[1-\xi (\mathbf{R})(1-r)]r$$

$$-2\pi /3\int \mathrm{d}\mathbf{R}({\nabla }_r^2{P}_{2,ep}^{FCI}(\mathbf{R},\mathbf{r})){|}_{r=0}\int \mathrm{d}r{e}^{-2{\beta }^2(\mathbf{R})r^2}[1-\xi (\mathbf{R})(1-r){]}^2{r}^3$$

$$+4\pi /3\int \mathrm{d}\mathbf{R}({\nabla }_r^2{P}_{2,ep}^{FCI}(\mathbf{R},\mathbf{r})){|}_{r=0}\int \mathrm{d}r{e}^{-{\beta }^2(\mathbf{R})r^2}[1-\xi (\mathbf{R})(1-r)]r^3.$$

$$E^{epc}=-\int \mathrm{d}\mathbf{R}\frac{{\rho }_e(\mathbf{R}){\rho }_p(\mathbf{R})}{a-b{\rho }_e^{1/2}(\mathbf{R}){\rho }_p^{1/2}(\mathbf{R})+c{\rho }_e(\mathbf{R}){\rho }_p(\mathbf{R})},$$

Full text

I Introduction

Density functional theory (DFT) Hohenberg and Kohn (1964); Parr and Yang (1994) is a powerful computational tool in chemistry and physics. Typically DFT calculations are conducted in the context of the Born-Oppenheimer approximation, where the nuclei are treated as fixed point charges when solving the electronic Schrödinger equation. However, the quantum mechanical nature of the nuclei and the non-Born-Oppenheimer effects between electrons and nuclei are important for many challenging physical and chemical problems, such as proton-coupled electron transfer Hammes-Schiffer (2015). Multicomponent DFT, which enables the quantum mechanical treatment of more than one type of particle, such as electrons and protons, was developed to tackle these problems Capitani et al. (1982); Shigeta et al. (1998); Kreibich and Gross (2001); Chakraborty et al. (2008). Analogs to the Hohenberg-Kohn theorems Hohenberg and Kohn (1964) have been proven for multicomponent systems Capitani et al. (1982); Kreibich and Gross (2001); Kreibich et al. (2008), and the Kohn-Sham formalism Kohn and Sham (1965), in which the kinetic energy functional is approximated with electronic and nuclear orbitals, has been extended to multicomponent systems Shigeta et al. (1998); Udagawa and Tachikawa (2006); Pak et al. (2007); Kreibich and Gross (2001); Kreibich et al. (2008); Chakraborty et al. (2008, 2009). This paper centers on the nuclear-electronic orbital (NEO) Webb et al. (2002) implementation of multicomponent DFT, denoted NEO-DFT Pak et al. (2007); Chakraborty et al. (2008), applied to systems in which specific hydrogen nuclei are treated quantum mechanically, and at least two other nuclei are treated as classical point charges to avoid difficulties with translations and rotations.

A major challenge in the Kohn-Sham implementation of multicomponent DFT is the design of accurate and practical electron-proton correlation functionals. When electron-proton correlation is neglected, the proton densities are highly over-localized, leading to unphysical results for most properties of interest. Efforts to approximate the electron-proton correlation functional by analyzing the connectivity between the electron-proton pair density and the total wave function or certain orbitals Kreibich and Gross (2001); Chakraborty et al. (2008); Sirjoosingh et al. (2011, 2012); Kreibich et al. (2008) have improved proton densities but are computationally expensive or not easily transferable to general systems. Recent attempts Imamura et al. (2008); Udagawa et al. (2014) to develop simple electron-proton correlation functionals related to the Colle-Salvetti formulation Colle and Salvetti (1975), which forms the basis of the Lee-Yang-Parr (LYP) electron-electron correlation functional Lee et al. (1988), relied on problematic approximations. For example, Ref. Imamura et al. (2008) has the incorrect sign of electron-proton correlation in the paper and incorrectly assumed a norm condition that leads to a singularity problem. Moreover, these previous attempts only applied those functionals as energetic corrections after the self-consistent-field (SCF) procedure. As a result, these approaches Imamura et al. (2008); Udagawa et al. (2014) do not influence the highly over-localized proton densities and thus are not useful for describing any properties that rely on the proton density. In this paper, we develop an electron-proton correlation functional based on the Colle-Salvetti formulation using well-justified approximations and apply this new functional within the SCF procedure, producing accurate proton densities at a low computational cost.

II Theory

Consider a multicomponent system with $N_{e}$ electrons and $N_{p}$ protons, where the subscripts $e$ and $p$ denote electrons and protons for all variables hereafter. Within the framework of multicomponent DFT, the total energy can be expressed as

[TABLE]

where the quantities in each parenthesis represent the non-interacting reference kinetic energy, external potential energy, mean-field Coulomb interactions, same-particle exchange-correlation (xc) energy, and electron-proton correlation (epc) energy, respectively. The electron and proton densities are defined in terms of the electron and proton orbitals, $\phi_{e}^{i}(\mathbf{r}_{e})$ and $\phi_{p}^{I}(\mathbf{r}_{p})$: $\rho_{e}(\mathbf{r}_{e})=\sum_{i}^{N_{e}}|\phi_{e}^{i}(\mathbf{r}_{e})|^{2}$, $\rho_{p}(\mathbf{r}_{p})=\sum_{I}^{N_{p}}|\phi_{p}^{I}(\mathbf{r}_{p})|^{2}$. The Kohn-Sham equations for electrons and protons obtained by minimizing the total energy are

[TABLE]

in atomic units with the effective potentials for the electrons and protons defined as

[TABLE]

Each potential term is obtained by taking the functional derivative of the corresponding energy term in Eq. (1) with respect to either the electron or the proton density. The electron-electron exchange-correlation functional is defined the same as in standard electronic DFT, and previously developed functionals can be used. Because quantum protons are spatially localized in molecular systems, the proton-proton exchange and correlation energies are negligible, and can be approximated with the diagonal Hartree-Fock exchange terms to eliminate self-interaction error. For systems with multiple protons, the proton orbitals can be assumed to be only singly occupied with high spin.

As discussed above, the epc functional $E^{epc}[\rho_{e},\rho_{p}]$ is the most critical term in the multicomponent Kohn-Sham formalism. The approach that completely neglects this term produces highly over-localized proton densities and therefore is not a useful tool for accurately describing the quantum behavior of protons. Herein we develop a multicomponent epc functional, denoted epc17, based on an analog of the Colle-Salvetti formulation but utilizing different approximations from previous attempts Imamura et al. (2008); Udagawa et al. (2014) in an effort to produce accurate proton densities. The epc17 functional is a multicomponent counterpart to the well-known LYP functional Lee et al. (1988) in conventional electronic DFT, excluding the density gradient terms in its current form, and it is implemented within the SCF procedure.

The derivation of our density functional begins with an analog of the Colle-Salvetti ansatz Colle and Salvetti (1975), which approximates the total wave function for a multicomponent electron-proton system as:

[TABLE]

where $\mathbf{x}_{i,e}$ denotes the combined electronic variable for the spatial coordinate $\mathbf{r}_{i,e}$ and spin $\sigma_{i,e}$, and only the spatial coordinate $\mathbf{r}_{I,p}$ is used for protons as the nuclei are assumed to be high spin. FCI denotes full configuration interaction for a particular type of particle, which takes into account all of the exchange and correlation effects between particles of this type but only includes a mean-field Coloumb potential from the other type of particle. Therefore, the product of the electron and proton FCI wave functions (i.e., Eq. (4) without the last correlation factor) is the mean-field FCI electron-proton wave function Cassam-Chenaï et al. (2017). The electron-proton correlation effects are included in the last Jastrow factor, with the function $\varphi(\mathbf{r}_{i,e},\mathbf{r}_{I,p})$ defined as

[TABLE]

where $r=|\mathbf{r}|=|\mathbf{r}_{i,e}-\mathbf{r}_{I,p}|$, and $\mathbf{R}$ is the center of mass for the two particles, which is approximately the same as $\mathbf{r}_{I,p}$ because of the significantly larger mass of the proton compared to the electron. In this equation, $\beta(\mathbf{R})$ and $\xi(\mathbf{R})$ are two functions that will be determined below. This ansatz is chosen because it fulfills the electron-proton cusp condition as $\mathbf{r}_{i,e}\rightarrow\mathbf{r}_{I,p}$ Bingel (1967) and reduces to the exact mean-field FCI electron-proton wave function when the electrons and protons are far apart.

Analogous to the original Colle-Salvetti formulation for electron-electron correlation, the electron-proton pair density matrix, which was the second-order reduced density matrix in the electronic case, can be approximated by the product of the mean-field pair density matrix and a correlation factor,

[TABLE]

where the mean-field electron-proton pair density matrix is the product of the associated electron and proton first-order reduced density matrices: $P_{2,ep}^{FCI}(\mathbf{x}_{1,e}^{\prime},\mathbf{r}_{1,p}^{\prime};\mathbf{x}_{1,e},\mathbf{r}_{1,p})=P_{1,e}^{FCI}(\mathbf{x}_{1,e}^{\prime};\mathbf{x}_{1,e})P_{1,p}^{FCI}(\mathbf{r}_{1,p}^{\prime};\mathbf{r}_{1,p})$. The electron-proton correlation energy is formally defined as the difference between the exact total energy associated with the wave function in Eq. (4) and the energies associated with the electron and proton mean-field FCI wave functions after removing the double counting of the mean-field electron-proton Coulomb interaction energy:

[TABLE]

Here $V_{ee}$ and $V_{pp}$ are two-particle Coulomb repulsion operators, $P_{1,e}$ and $P_{2,e}$ are the electron first- and second-order reduced density matrices associated with the total wave function in Eq. (4), and $P_{1,p}$ and $P_{2,p}$ are the proton first- and second-order reduced density matrices defined analogously. Moreover, the pair density $P_{2,ep}(\mathbf{x}_{1,e},\mathbf{r}_{1,p})$ denotes the diagonal elements of the electron-proton pair density matrix. Assuming that the electron and proton first- and second-order reduced density matrices are reasonably well approximated by those from mean-field FCI, the electron-proton correlation energy reduces to the last term in Eq. (7):

[TABLE]

Equating the electron and proton densities associated with the full wave function in Eq. (4) to those associated with the mean-field FCI electron-proton wave function, which is a direct consequence of the above approximation for the second-order reduced density matrices, leads to an approximate relationship between the undetermined functions $\xi(\mathbf{R})$ and $\beta(\mathbf{R})$:

[TABLE]

In the original electronic Colle-Salvetti formulation, $\beta(\mathbf{R})$ is assumed to be inversely proportional to the Wigner-Seitz radius $r_{s,e}\propto\rho_{e}^{-1/3}(\mathbf{R})$ from a uniform electron gas model. In the multicomponent electron-proton case, however, $\beta(\mathbf{R})$ must depend on both the electron and proton densities. Thus, we define a geometric mean Wigner-Seitz radius $r_{s,ep}=(r_{s,e}r_{s,p})^{1/2}$, where $r_{s,p}$ is defined analogously to $r_{s,e}$, and assume $\beta(\mathbf{R})$ to be inversely proportional to the local mean $r_{s,ep}(\mathbf{R})$

[TABLE]

Note that Refs. Imamura et al. (2008); Udagawa et al. (2014) assumed that $\beta(\mathbf{R})\propto\rho_{e}^{1/3}(\mathbf{R})$ in deriving their expressions for electron-proton correlation energy; however, such neglect of the dependence on the proton density does not seem to be physically reasonable.

Although $\xi(\mathbf{R})$ and $\beta(\mathbf{R})$ have been determined, the epc functional in Eq. (8) is still not practical because it involves an expensive six-dimensional integration over both electron and proton coordinates. We continue following the strategy of the Colle-Salvetti formalism and transform the coordinates {$\mathbf{x}_{1,e},\mathbf{r}_{1,p}$} to {$\mathbf{R}$,$\mathbf{r}$}. After making use of the mathematical relation

[TABLE]

and Eq. (5), the electron-proton correlation energy becomes

[TABLE]

In this paper, we neglect the last two terms involving the second-order derivatives of the pair density and leave the inclusion of these terms for future development. Therefore, we are able to obtain a simple form for the epc functional. Substituting Eq. (9) into the first two terms of Eq. (12) and integrating leads to the final form of our epc17 functional, which we approximate as:

[TABLE]

with three parameters to be determined below. Although the derivation of the epc17 functional is based on the wave function in Eq. (4) and its associated density matrices, these are never explicitly constructed. This simple type of functional has only two fundamental variables, namely electron density $\rho_{e}$ and proton density $\rho_{p}$. Therefore, it is a local density approximation (LDA) form for both electrons and protons. Note that $\rho_{e}$ is the total electron density for both closed-shell and open-shell cases, and $\rho_{p}$ is the total proton density, where the protons are assumed to be high spin.

III results

We have implemented the epc17 functional in GAMESS Schmidt et al. (1993) and applied this approach to the molecules \ceFHF- and \ceHCN, which have been extensively studied for testing purposes Sirjoosingh et al. (2015); Brorsen et al. (2015); Culpitt et al. (2016). The heavy nuclei are treated as point charges fixed to the positions obtained by optimizing the geometry at the conventional electronic DFT level, and the hydrogen nucleus, as well as all electrons, are treated quantum mechanically. The def2-QZVP Weigend and Ahlrichs (2005) electronic basis set and an even-tempered 8s8p8d proton basis set centered at the optimized hydrogen position with exponents running from $2\sqrt{2}$ to $32$ are used for all calculations. The B3LYP functional Becke (1988); Lee et al. (1988) is used for electron-electron exchange-correlation. Our epc17 functional with parameters $a$, $b$, and $c$ optimized to be 2.35, 2.4, and 3.2, is used for electron-proton correlation. Note that the parameters are fit to accurate proton densities for \ceFHF-, analogous to the Colle-Salvetti electron correlation formula, where the parameters were fit to the properties of the He atom. As for most parameterized electronic exchange-correlation functionals, the values of the parameters could significantly influence the final results. This parameter dependence will be investigated more thoroughly in future work.

This NEO-DFT/epc17 method is compared to the NEO-DFT/no-epc method, which includes no electron-proton correlation, as well as a grid-based method Marston and Balint-Kurti (1989) that is considered to be a reference for these electronically adiabatic systems. In the grid-based method, the total electronic energy is calculated using conventional electronic DFT/B3LYP for the hydrogen nucleus positioned at each grid point on a three-dimensional grid spanning the relevant region for the proton density, and the three-dimensional Schrödinger equation is solved numerically for the proton using the Fourier grid Hamiltonian method Marston and Balint-Kurti (1989). The proton density is defined to be the square of the proton wave function calculated in this manner. The grid method is numerically exact when the Born-Oppenheimer approximation is valid, as for these model systems.

Figure 1 shows the on-axis and off-axis proton densities for \ceFHF-, where the on-axis proton density is the slice along the axis connecting the two heavy nuclei, and the off-axis proton density is a slice perpendicular to this axis passing through the midpoint of the F—F bond. The NEO-DFT/no-epc method predicts a highly over-localized proton density, whereas the NEO-DFT/epc17 method is in excellent agreement with the grid-based reference method. Despite visual appearances, all of these proton densities are normalized in three-dimensional space, where multiplication by the volume element more heavily weights the regions further from the origin. In contrast, the results from previous papers were rescaled to be normalized in one dimension for the on-axis density Sirjoosingh et al. (2015); Culpitt et al. (2016); Brorsen et al. (2015) but contained significant errors because of a poor description of the off-axis proton density. The current epc17 functional corrects this problem and exhibits excellent agreement for both on-axis and off-axis proton densities with proper three-dimensional normalization. To our knowledge, no other electron-proton correlation functional can describe the proton densities with proper normalization in three dimensions even qualitatively accurately.

Without further parameterization, we applied the NEO-DFT/epc17 method to \ceHCN, an asymmetric case where the hydrogen atom is located at one end of the molecule. The results for this molecule are shown in Figure 2. The NEO-DFT/epc17 method also provides greatly improved results compared to the NEO-DFT/no-epc method for both on-axis and off-axis proton densities, although minor differences can be observed. Note that altering the $a$, $b$, and $c$ parameters within the epc17 functional did not lead to significant improvement, implying transferability across different molecules for this functional, although additional molecules will be studied for further testing. Inclusion of the last two terms in Eq. (12) may further improve this type of asymmetric case.

We investigated the reason for the more delocalized proton densities obtained with the epc17 functional by plotting the epc proton potential $v^{epc}_{p}$ as a function of the electron and proton densities in Figure 3. The ranges of the variables $\rho_{e}$ and $\rho_{p}$ are chosen to be in the physical regime for hydrogen atoms in a typical molecular environment. This plot shows that when the proton density increases, the epc potential for the proton also increases. Driven by this extra potential, the proton tends to become more delocalized to minimize the area with high proton density and hence decrease the potential energy. This property explains why such a simple LDA functional can effectively delocalize the proton density.

Thus, this LDA type of electron-proton correlation functional should be viewed as a first important step, with the understanding that gradient corrections will be critical for future functional development. Moreover, because the high and low density limits of electron-proton correlation are unknown, the parameterized form developed herein focuses on the physically meaningful regions for molecular systems, and the asymptotic high and low density limits are most likely not correct. Further analysis of these limits and development of functional forms that satisfy such constraints are other directions for future development.

In addition, similar to the LYP functional, the epc17 functional does not explicitly include the kinetic correlation effects. These effects could be included with an adiabatic connection formula, Langreth and Perdew (1975); Imamura et al. (2002) as in a previous electron-proton correlation functional from our group. Sirjoosingh et al. (2011, 2012) However, because the epc17 functional is parameterized to fit numerically exact densities or energies, these effects are implicitly included in a semiempirical manner. Another challenge in the development of electronic functionals has been the accurate inclusion of exchange terms for electrons. In contrast to electronic functionals, however, electron-proton functionals do not need to consider exchange between an electron and a proton because they are not identical particles. Given the extensive literature aimed at improving electronic functionals, future efforts will explore alternative models and treatments for electron-proton correlation functionals.

IV Conclusions

In this paper, we developed an electron-proton correlation functional analogous to the LYP electron correlation functional based on the Colle-Salvetti formulation. The implementation of this epc17 functional within the SCF procedure for multicomponent DFT produces accurate proton densities that are in good agreement with reference data for two representative molecules. To our knowledge, such accurate proton densities have not been achieved with any previous multicomponent DFT formulations. Moreover, the NEO-DFT/epc17 method is computationally inexpensive with the same formal scaling as conventional electronic DFT and is therefore promising for a wide range of future applications. For example, recent work based on the epc17 functional has been used to calculate proton affinities and to investigate the impact of proton delocalization on optimized geometries. Brorsen et al. (2017); Foo In addition, the NEO-DFT/epc17 method will enable the inclusion of nuclear quantum effects in calculations of pKa’s, reaction paths, and reaction dynamics, as well as tunneling splittings and vibronic couplings. Therefore, this development lays the foundation for a wide range of potential new directions.

Acknowledgements.

This work was supported by the National Science Foundation Grant No. CHE-1361293.

Supplementary material

See supplementary material for the complete derivation of the epc17 functional.