Singling Out Dynamic and Nondynamic Correlation

TL;DR

This paper analyzes the correlation components of pair density to distinguish between dynamic and nondynamic correlation, revealing key factors for describing dispersion forces and guiding the development of advanced electronic structure methods.

Contribution

It introduces a method to classify molecular systems based on correlation components and identifies the key factors influencing dispersion forces and correlation types.

Findings

Correlation parts of pair density can be separated into short-range and long-range components.

Long-range asymptotics reveal the universal decay related to London dispersion forces.

Identification of correlation types aids in developing new electronic structure methods.

Abstract

The correlation part of the pair density is separated into two components, one of them being predominant at short electronic ranges and the other at long ranges. The analysis of the intracular part of these components permits to classify molecular systems according to the prevailing correlation: dynamic or nondynamic. The study of the long-range asymptotics reveals the key component of the pair density that is responsible for the description of London dispersion forces and a universal decay with the interelectronic distance. The natural range-separation, the identification of the dispersion forces, and the kind of predominant correlation type that arise from this analysis are expected to be important assets in the development of new electronic structure methods in wave function, density, and reduced density-matrix functional theories.

Equations

$${\rho }_2(1,2)=\frac{N(N-1)}{2}\int d_3\dots d_{\mathbf{n}}{\left|\Psi (1,\dots ,\mathbf{n})\right|}^2,$$

$$I({\rho }_2,s)=\int d_1{d}_2{\rho }_2(1,2)\delta (s-r_{12}),$$

$$V_{ee}\left[I\right]=\int ds\frac{I({\rho }_2,s)}{s}.$$

$$\Delta {\rho }_2^{\text{c}}(1,2)={\rho }_2(1,2)-{\rho }_2^{\text{HF}}(1,2).$$

$$h_{\text{c}}(s)=I(\Delta {\rho }_2^{\text{c}},s)=\int d_1{d}_2\Delta {\rho }_2^{\text{c}}(1,2)\delta (s-r_{12}).$$

$${\rho }_2^{\text{SD}}({\rho }_1,1,2)={\rho }_1(1){\rho }_1(2)-{\left|{\rho }_1(1;2)\right|}^2,$$

$${\rho }_2^{\text{HF}}(1,2)={\rho }_2^{\text{SD}}({\rho }_1^{\text{HF}},1,2),$$

$$\Delta {\rho }_2^{\text{c}}(1,2)$$

$${\rho }_2(1,2)={\rho }_2^{\text{HF}}(1,2)+\Delta {\rho }_2^{{\text{c}}_I}(1,2)+\Delta {\rho }_2^{{\text{c}}_{II}}(1,2),$$

$$h_{\text{c}}(s)=h_{{\text{c}}_I}(s)+h_{{\text{c}}_{II}}(s)=I(\Delta {\rho }_2^{{\text{c}}_I},s)+I(\Delta {\rho }_2^{{\text{c}}_{II}},s),$$

$$\hat{H}=-t\sum _{⟨\mu ,\nu ⟩,\sigma }\left({\hat{c}}_{\mu \sigma }^{†}{\hat{c}}_{\nu \sigma }+{\hat{c}}_{\nu \sigma }^{†}{\hat{c}}_{\mu \sigma }\right)+U\sum _{\mu }{\hat{\rho }}_{\mu \alpha }{\hat{\rho }}_{\mu \beta },$$

Full text

\alsoaffiliation

Kimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), Donostia, Euskadi, Spain \altaffiliation Equally contributed

\alsoaffiliationKimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), Donostia, Euskadi, Spain \alsoaffiliationInstitut de Química Computacional i Catàlisi (IQCC) and Departament de Química, University of Girona, C/ Maria Aurèlia Capmany, 69, 17003 Girona, Catalonia, Spain \altaffiliation Equally contributed

\alsoaffiliationKimika Fakultatea, Euskal Herriko Unibertsitatea (UPV/EHU), Donostia, Euskadi, Spain

\alsoaffiliationIKERBASQUE, Basque Foundation for Science, 48013 Bilbao, Euskadi, Spain.

1 Introduction

Electron correlation being the holy grail of electronic structure methods, it has been the subject of extended analysis. 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14 The solution of quantum many-body problems hinges on the type of correlation present in the system, and one of the most practical classifications consists in the separation between dynamic- and nondynamic-correlation-including methods. Indeed, there are accurate methods to study systems with one predominant correlation type, but systems presenting both correlation types pose one of the greatest current challenges in electronic structure theory.15, 16, 17

The attempt at taking the best of both worlds has led to a resurgence of interest in hybrid schemes,18 merging methods that recover different correlation types.19, 20, 21 Among hybrid implementations, the most successful one is based on the range separation of electron correlation,18, 22, 23 using a mixing function to combine approximations that account for short-range dynamic correlation —such as density functional approximations— with approaches providing correct long-range asymptotics. The performance of these methods pivots on the choice of the function combining the two approaches, which provides a natural splitting of the Coulomb interaction and thus the pair density.24 In range-separation approximations, the typical choice is the error function that, in turn, depends on an attenuating parameter, which is both system- and property-dependent.25, 26 Even though the methods are chosen according to their ability of recovering dynamic and nondynamic correlation, the range-separation of the pair density has not been motivated by the correlation type present in the system, risking double counting of electron correlation.

Thus far, there has been very few attempts to separate dynamic and nondynamic correlation,2, 6, 7, 8, 9, 27, 4, 5, 13 most of them based on energy calculations. The lack of a physically sound separation of dynamic and nondynamic correlation precludes individual treatment of these effects. We analyze the decomposition of the pair density into three components: the uncorrelated reference and two correlation terms. The latter two behave differently with respect to large changes of the first-order reduced density matrix (1-RDM), permitting the identification of systems with prevalent dynamic or nondynamic correlation.2, 6, 7, 11, 12 Some of us have recently used a similar strategy to obtain scalar 11 and local 12 measures of dynamic and nondynamic electron correlation from a two-electron model. The intracule of the correlation components of the pair density yields a two-fold separation of the Coulomb hole in terms of correlation type and interelectronic range. These components of the pair density display a simpler mathematical form than the total pair density, one of them being dominant at short ranges and one with prevailing long-range contributions. This feature is particularly convenient for the design of energy functionals in wavefunction, density and density matrix functional theories. As a result of this separation, we will clearly identify the part of the pair density that is responsible for the correct description of van der Waals interactions and unveil a universal condition it should satisfy.28 To our knowledge, the latter is the only known condition of the pair density that can be employed to design methods including van der Waals interactions.

2 Theoretical background

Let us consider the pair density of a $N$-electron system described by the $\Psi(\mathbf{1},\ldots,\mathbf{n})$ wavefunction,

[TABLE]

where numerical variables ($\mathbf{1},\mathbf{2},\ldots$) refer to space and spin coordinates. Upon integration over its coordinates, the pair density can be reduced to the intracule density, which only depends on the interelectronic range separation, $s$,

[TABLE]

where $r_{12}$ is the Euclidean distance between the electrons at $\mathbf{1}$ and $\mathbf{2}$. The intracule density is the simplest function in terms of which we can express the Coulomb interaction energy,

[TABLE]

The electron correlation contents of the pair density can be determined by the difference between the actual pair density and an uncorrelated reference, which here we choose to be the Hartree-Fock (HF) one,

[TABLE]

The intracule of this function is Coulson’s Coulomb hole,29

[TABLE]

In order to split the correlation part of the pair density (Eq. 4) we employ an approximate pair density, the single-determinant (SD) ansatz of the pair density,1

[TABLE]

where $\rho_{1}(\mathbf{1};\mathbf{2})$ is the 1-RDM and $\rho_{1}(\mathbf{1})\equiv\rho_{1}(\mathbf{1};\mathbf{1})$ is the electron density. Substituting $\rho_{1}$ by the HF 1-RDM in Eq. 6, yields the HF pair density, i.e.,

[TABLE]

which does not account for electron correlation. However, $\rho_{2}^{\text{SD}}(\rho_{1},\mathbf{1},\mathbf{2})$ can be regarded as an approximation to the actual pair density; an approximation which does not account for dynamic correlation either at short 30 or at long range.28 Figure 1 depicts the two paths of arriving at the exact $\rho_{2}(\rho_{1},\mathbf{1},\mathbf{2})$ from $\rho_{2}^{\text{SD}}(\rho_{1}^{\text{HF}},\mathbf{1},\mathbf{2})$, either straightforwardly or through the intermediate SD approximation. The latter path defines the decomposition of the correlation part of the pair density,

[TABLE]

$\Delta\rho_{2}^{\text{c}_{I}}(\mathbf{1},\mathbf{2})$ will be large only if the HF 1-RDM and the actual 1-RDM are significantly different and, in such case, the system will be affected by nondynamic correlation. Indeed, the wavefunction of systems dominated by dynamic correlation can be described by a large expansion of Slater determinants with one of them (the HF one) having an expansion coefficient very close to one.31 Therefore, these systems are characterized by a 1-RDM that retains the shape of the HF 1-RDM. Conversely, the wavefunction of nondynamic-correlated systems can be written as a shorter expansion of Slater determinants, but in this case the HF determinant has an expansion coefficient that is qualitatively smaller than one.31 Since the 1-RDM is determined by the square of the expansion coefficients, we expect systems affected by nondynamic correlation to display large $\Delta\rho_{2}^{\text{c}_{I}}(\mathbf{1},\mathbf{2})$. Some authors have used similar arguments to use the electron density (the diagonal part of the 1-RDM) as a means to define dynamic and nondynamic correlation energy.2, 6 In this work, we prefer to employ the 1-RDM because the cases of spin entaglement would not be regarded as nondynamic correlation if only density differences were considered. Indeed, in the stretched H2 molecule, the HF electron density is qualitatively similar to the exact one, whereas there are large and notorious differences between the exact and the HF 1-RDMs.

The magnitude of $\Delta\rho_{2}^{\text{c}_{I}}(\mathbf{1},\mathbf{2})$ can be thus regarded as a measure of nondynamic correlation but it can also be interpreted as the correlation retrieved by using the actual 1-RDM rather than the HF one to construct the pair density. Conversely, $\Delta\rho_{2}^{\text{c}_{II}}(\mathbf{1},\mathbf{2})$ does not depend on the differences between $\rho_{1}^{\text{HF}}$ and $\rho_{1}$, but on the validity of the SD approximation. Note that $\Delta\rho_{2}^{\text{c}_{II}}$ coincides with the cumulant of the pair density.32, 14 The intracule functions of $\rho_{2}^{\text{SD}}(\rho_{1},\mathbf{1},\mathbf{2})$ and the exact pair density, $\rho_{2}(\rho_{1},\mathbf{1},\mathbf{2})$, display the same asymptotic behavior 33 and, therefore, $\Delta\rho_{2}^{\text{c}_{II}}$ is dominated by the short-range component. Interestingly, $\Delta\rho_{2}^{\text{c}_{I}}(\mathbf{1},\mathbf{2})$ is the long-range-dominant component of the correlated part of the pair density (Eq. 4) because the HF and the exact 1-RDM can differ substantially at large separations, for instance, in the presence of entanglement. On the contrary, $\Delta\rho_{2}^{\text{c}_{I}}(\mathbf{1},\mathbf{2})$ displays very small values at small interelectronic distances mostly due to the opposite-spin part of this term.

The current partition,

[TABLE]

provides a natural range separation of the pair density that can be employed to split the Coulomb hole into two correlation components,

[TABLE]

naturally yielding a separation of electron correlation by range. We will show that the decay of $I(\Delta\rho_{2}^{\text{c}_{II}},R)$ is universal and it corresponds to a characteristic signature of London dispersion forces ($R$ being the distance between two atoms in the molecule).

3 Results and Discussion

In the following we introduce five selected examples that illustrate the effectiveness of the current scheme to separate the correlation part of the Coulomb hole at different ranges and how the long-range of $\Delta\rho_{2}^{\text{c}_{II}}$ can be used to identify and characterize van der Waals interactions.

The Hydrogen Molecule34.— At the equilibrium geometry, $h_{\text{c}_{II}}(s)$ dominates over $h_{\text{c}_{I}}(s)$ at all interelectronic distances $s$, as shown in the left panel of Fig. 2, whereas $h_{\text{c}_{I}}(s)$ increases importantly as the bond is stretched, in line with the expected increase of nondynamic correlation. The most likely distribution of the electron pair at large bond lengths corresponds to one electron sitting at each atom and, accordingly, the intracule density peaks around the bond-length distance. At the dissociation limit, the long-range part of the Coulomb hole is completely determined by $h_{\text{c}_{I}}(s)$ because one isolated electron cannot give rise to dynamic correlation. Hence, the unrestricted HF calculation of H2 produces Coulomb hole components that are not distinguishable from FCI.35 A simple interpretation is also obtained from valence bond theory: at large separations, the exact pair density is entirely described by covalent components, whereas the HF pair density contains equally contributing ionic and covalent terms. $h_{\text{c}_{II}}(s)$ removes the ionic contribution (i.e., removes contributions keeping the electrons at short distances), whereas $h_{\text{c}_{I}}(s)$ adds the missing covalent contribution (i.e., adds contributions placing one electron in each atom); in accord with the results plotted in the r.h.s. of Fig. 2 (see also Supp. Material).

The Hubbard Dimer.— The Hubbard dimer is the simplest model of interacting particles in a lattice and conceivably the most studied model for testing methods at different regimes.36, 37 We employ the one-dimension Hamiltonian of the Hubbard model,

[TABLE]

where $\mu$ and $\nu$ denote the sites, $\sigma$ the spin polarization ($\alpha$ or $\beta$), $\hat{c}^{\dagger}_{\mu\sigma}$ and $\hat{c}_{\mu\sigma}$ are creation and annihilation operators of one electron with spin $\sigma$ in site $\mu$, and $\hat{\rho}_{\mu\sigma}$ stands for a one-particle number operator with spin $\sigma$ acting on site $\mu$. $t$ is the hopping parameter and $U$ is the on-site interaction parameter. These parameters control the electron correlation within the Hubbard model, small (large) $U/t$ inducing dynamic (nondynamic) correlation. Hence large $U/t$ values prompt the electrons to distribute among the sites to minimize the electron repulsion. Fig. 3 presents plots of the Coulomb hole at various values of $U/t$ for the two-electron two-site Hubbard model in real space.36 At low $U/t$ values, the system is barely affected by correlation, thus dynamic correlation dominates (small $h_{\text{c}_{I}}(s)$ and large $h_{\text{c}_{II}}(s)$) and the electron pairs distribute equally between on-site and intersite components. As $U/t$ grows, nondynamic correlation dominates and $h_{\text{c}_{I}}(s)$ becomes more important, being the prevailing contribution between sites.

The He series.34— The He isoelectronic series is perhaps the simplest series of systems dominated by dynamic correction.38 As the atomic number $Z$ increases, the electron correlation of He($Z$) tends to a constant and the exact electron density barely distinguishes from the HF one. In Fig. 4 we observe that $h_{\text{c}_{I}}(s)$ decreases with the atomic number $Z$ and, hence, $h_{\text{c}_{II}}$ completely takes over.

$N$* hydrogen atoms*.34— The size consistency of our approach and its ability to measure spin entanglement is examined in Fig. 5. We have plotted the Coulomb hole of the $N$-vertex polyhedron resulting from $N$ hydrogen atoms separated by 10$\AA\,$ from the center of the polyhedron. At these large separations, the hydrogen atoms only interact to each other through entanglement and this is the only term that remains in the cumulant,4 (i.e. in $\Delta\rho_{2}^{\text{c}_{II}}$), which shows a linear behavior with $N$ (see Fig. 5). As in previous systems, $h_{\text{c}_{II}}$ is short ranged and its contribution to the energy grows linearly with $N$. These systems can be classified as nondynamic correlated because $h_{\text{c}_{I}}$ is mostly long ranged and peaks at the same positions of the intracule density maxima. The planar $D_{4h}/D_{2h}$ potential energy surface of H4 has also been used for discriminating between dynamic and nondynamic correlation 39 and is given in the Supp. Material.

van der Waals (vdW) Interactions.34— Fig. 6 includes plots of the Coulomb hole of the helium dimer. $h_{c}$ compares satisfactorily to earlier calculations.40 The dynamic long-range interaction between the two noble-gas atoms is reflected by the second peak of the intracule density, whereas the interaction of the electron pair within each helium shows in the first peak. Regardless the bond length, $h_{\text{c}_{II}}$ dominates, indicating that the correlation is dynamic and mainly affects the electron pair within each He. Unlike H2, there is very little long-range nondynamic correlation in this system; however, at all distances, the long-range part of $h_{\text{c}_{II}}$ peaks around the bond-length distance (see the inset plots of Fig. 6). The plot in Fig. 7 presents $h_{\text{c}_{II}}(R)$ against the bond length, $R$, revealing a $R^{-3}$ decay. It is a textbook fact that the pairwise vdW energy decays like $R^{-6}$.41 Using perturbation theory, we have recently proved that the vdW contribution to $h_{\text{c}_{II}}$ should actually decay like $R^{-3}$, the integration of $h_{\text{c}_{II}}(s)/s$ over $s$ yielding a fraction of the Coulombic interaction (Eq. 3) due to London dispersion forces and, therefore, decaying as $R^{-6}$.28 Fig. 7 includes plots for other noble-gas dimers, which also satisfy this property. Most density functional theory (DFT) practitioners add ad hoc empirical corrections to the energy for vdW interactions and, therefore, they only shift the relative energies of different conformers, yet the electronic structure of the system is not completely considered.42 The present separation into correlation regimes unveils the target part of the pair density and the Coulomb hole, i.e., the long-range component of $h_{\text{c}_{II}}(s)$, which should be improved in order to incorporate the description of London dispersion forces and avoid the latter problem, thus opening a door to the accurate account of these forces within DFT and reduced density matrix functional theory (RDMFT).

In conclusion, Eqs. 8 and 9 represent a separation of the pair density and the Coulomb hole into components dominated by short- and long-range interactions. This result is expected to be important in the development of new hybrid electronic structure methods that can be employed in RDMFT 43, 21 and other computational approaches. For instance, the HF reference in Fig. 1 can be replaced by the Kohn-Sham system to adapt the present idea to DFT. It can be shown that the exchange-correlation functional can be entirely written in terms of the Kohn-Sham orbitals, $\Delta\rho_{2}^{\text{c}_{I}}$, and $\Delta\rho_{2}^{\text{c}_{II}}$. Hence, a template to construct density functional approximations, where the correlation components are treated separately, arises. Such possibility is already being explored in our laboratory.

Acknowledgements

We thank D. Casanova, X. Lopez, P. Salvador, and specially, P.M.W. Gill and J.M. Ugalde for helpful insights. This research has been funded by Spanish MINECO/FEDER Projects CTQ2014-52525-P, PGC2018-098212-B-C21, and EUIN2017-88605. We acknowledge doctoral grants BES-2015-072734 and FPU-2013/00176, and the funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement (No. 660943).

Supporting Information Available: Analysis of the $D_{2h}/D_{4h}$ H4 planar molecule and the hydrogen molecule from ionic and covalent contributions. Full plots of H2, the two-site Hubbard dimer, the He isoelectronic series, and the $N$ hydrogen atoms model.