Local Enhancement of Dynamic Correlation in Excited States: Fresh Perspective on Ionicity and Development of Correlation Density Functional Approximation Based on the On-Top Pair Density.

We discuss the interplay between the nondynamic and dynamic electron correlation in excited states from the perspective of the suppression of dynamic correlation (SDC) and enhancement of dynamic correlation (EDC) effects. We reveal that there exists a connection between the ionic character of a wave function and EDC. Following this finding we introduce a quantitative measure of ionicity based solely on local functions without referring to valence bond models. The ability to recognize both the SDC and EDC regions underlies the presented method, named CASΠDFT, combining complete active space (CAS) wave function and density functional theory (DFT) via the on-top pair density (Π) function. We extend this approach to excited states by devising an improved representation of the EDC effect in the correlation functional. The generalized CASΠDFT uses different DFT functionals for ground and excited states. Numerical demonstration for singlet π → π* excitations shows that CASΠDFT offers satisfactory accuracy at a fraction of the cost of the ab initio approaches.

Excited states of a singlet multiplicity pose a particular challenge to both ab initio and density functional theory (DFT) methods because of their multireference character and the related diverse electronic structure.[1−3] In π-conjugated systems it is customary to use the classification into states of a covalent and ionic character which originates from the valence bond theory (VB) description.[4] The π → π* singlet excitations to ionic states are problematic even for methods which recover dynamic correlation effects, because their exact treatment requires capturing both the effect of the dynamic σ-polarization[5] and the resulting spatial contraction of the atomic p-orbitals.[6,7]

The complete active space (CAS) method[8] is well-suited to capture the basic electron effects of the lowest excitations in molecules. A relatively small active space is sufficient to account for configurations which are the most relevant for the studied excited state. However, the CAS excitation energy obtained as a difference between excited-state and ground-state self-consistent CAS (CASSCF) energies is impaired with the serious neglect of dynamic correlation. This neglect is more pronounced in the excited states than in ground states; thus, CASSCF consistently overestimates the excitation energies.

A viable way to efficiently account for the correlation energy is offered by DFT approximations. Over the years various approaches to merge wave function theory (WFT) with DFT have been developed.[9,10] Although several studies applied combined methods to excited states, they have exclusively employed ground-state functionals.[11−16] The extended Hohenberg–Kohn theorem implies that exchange–correlation energy functionals for excited states within a given symmetry are bound to be different from the ground-state one; however, this path seems hardly possible to follow in pure DFT.[17,18] A pragmatic answer may lie in the combined WFT-DFT methods which put less demand on their DFT part when applied to excited states. Thus, one can hope that in WFT-DFT methods dynamic correlation for various excited states can be represented with a common approximate functional.

In our previous work[19] we proposed the combined CASΠDFT method, a unique feature of which is that the modified correlation functional is sensitive to the nondynamic correlation included in the WF part. This is achieved through a correction function that involves the on-top pair density and locally adjusts the dynamic correlation energy to account for suppression of dynamic correlation (SDC) and enhancement of dynamic correlation (EDC). To date, the method has been applied to ground states for which a proper representation of the SDC effect is crucial.[20,21] A pilot study of prototypical singlet π → π* excitations[22] has indicated that a CASΠDFT description of excited states may be generally feasible.

In this Letter we recognize that excited-state wave functions feature spatial regions of enhanced dynamic correlation that are physically relevant (the electron density is high). We demonstrate a direct connection between the presence of the EDC regions and the occurrence of ionic structures in the wave function. On the basis of this observation, we introduce a quantitative measure of ionicity. Better understanding of the EDC effect also motivates an extension of the CASΠDFT method to excited states, which employs different correction functions for ground and excited states.

Connection between Enhancement of Dynamic Correlation and Ionic Character of States. In a reliable combined WFT-DFT method, the correlation functional should be sensitive to the amount of electron correlation already accounted for by the assumed WFT. A local quantity which can convey information about electron correlation from WFT to DFT is the conditional density ρcond. It is defined as a ratio of the on-top pair density[23](where N stands for a number of electrons and σ indicates a spin coordinate) to the electron density, ρcond(r) = Π(r)/ρ(r), and yields electron density at r on condition that a reference electron occupies the same position. In the case of a single-determinantal closed shell wave function, the conditional density coincides everywhere with half of the density, i.e. ∀ρcond(r) = ρ(r)/2. If a wave function is multiconfigurational (MC), the corresponding conditional density can be less than, equal to, or greater than conventional density. In the first case, ρcond[ΨMC](r) < ρ(r)/2, electron correlation leads to instantaneous depletion of density; in the third case, ρcond[ΨMC](r) > ρ(r)/2, the conditional density is increased. Depletion of ρcond with respect to the conventional density ρ(r)/2 caused by nondynamic correlation leads to a local reduction of dynamic correlation. This effect has been called suppression of dynamic correlation (SDC).[19] In the spatial regions where the conditional density is increased, the dynamic correlation energy is also increased, and one refers to the enhancement of dynamic correlation effect (EDC). We arrive at a conclusion that the ratio of ρcond(r) to ρ(r)/2, given also by the on-top pair density (cf. eq )indicates either the SDC (X[ΨMC](r) < 1) or the EDC (X[ΨMC](r) > 1) effect. Despite the fact that the ratio X(r) has already been used in other methods combining WFT with DFT[9,11,23,24] and its behavior has been extensively studied,[25,26] no particular physical meaning has been assigned to spatial regions where X exceeds 1. In ref (23), the possibility that X (r) > 1 was considered an “interesting mathematical complication”. In ref (11), X(r) was truncated at 1 for practical reasons, i.e., to avoid complex numbers in the function responsible for effective spin polarization. Only recently has it been reported that spatial regions where X (r) > 1 may be related to ionic character of excited states.[19,22] Below we rigorously validate and generalize this observation, demonstrating that there exists a general relation between the ionic structure of a wave function and the presence of large-electron-density regions where X(r) exceeds the value of 1.

A model CAS wave function, the active part of which is ionic, may be represented as , where |...| denotes a Slater determinant. To reveal the ionic character, we switch to the basis of atomic-like orthogonal orbitals, and by employing the transformation and , we arrive atwhere the orbitals τ are localized on nuclei a and b, and α and β stand for spins. It is straightforward to show (cf. the Supporting Information for details) that for the CAS wave function with the active part given in eq , the ratio XCAS(r)is greater than 1 in the valence regions where τ(r)2/τ(r)2 ≪ 1, i.e. in the vicinity of a nucleus (analogously, the inequality holds in the vicinity of b, where τ(r)2/τ(r)2 ≪ 1). The function F is expressed by electron densities originating from the inactive and active parts of the wave function, F[ρ](r) = ρinact(r)2/2 + ρinact(r)ρact(r). Therefore, for ionic states X > 1 in regions of relatively high electron density. Because such a behavior of the ratio X indicates local ionicity resulting from “squeezing” of electrons (conditional density is greater than the conventional density) and enhancement of dynamic correlation, EDC, then the more general conclusion is that ionic structure of the wave function leads to EDC.

This finding is illustrated in Figure which presents the behavior of the X ratio for the 11B1 state of ethene, recognized as a paradigm of ionic states.[6] In accordance with eq , the ratio exceeds the value of 1, approaching the value of 2 in the vicinity of carbon atoms. Conversely, for the ground state we observe that X < 1 in out-of-core regions, which is a known general behavior of the ratio for ground (covalent) states.[25,26] Note that in both ground and excited states X = 1 in the core region.

Figure 1

X(r) ratio for the ground (11A) state and excited (11B1u) state of the ethene molecule. The molecule lies in the yz plane; carbon atoms are oriented along the z-axis and located at zC = ± 1.263 bohr. Left: curves plotted along the r = (x, 0, zC) line. Right: values of X(r) at r = (x, y, z), where z ∈ (zC – 1, zC + 1) points are depicted; “covalent” and “ionic” panels pertain to ground and excited states, respectively. Based on CAS(2,2)SCF calculations in the TZVP basis set.

In Figure we illustrate the behavior of the X(r) ratio obtained for the 11Bu states of butadiene and hexatriene molecules, described with CAS(4,4) and CAS(6,6) functions, respectively. Notice that in the case of butadiene all four carbon atoms are surrounded by the EDC regions, while for hexatriene it is true only for the two middle atoms. The model of ionicity assumed for a simple CAS(2,2) function extended to larger CAS functions explains these observations.

Figure 2

X(r) ratio along the r = (xC, yC, z) direction, where (xC, yC,0) is a position of the C carbon nucleus, plotted for the excited 11Bu states of E-butadiene (C4H6, left) and all-E-hexatriene (C6H8, right). The molecules lie in the xy plane. Based on CAS(4,4)SCF (left) and CAS(6,6)SCF (right) calculations in the TZVP basis set.

Consider first the CAS(4,4) 11Bu wave function of butadiene. The orbitals π and π* involved in the π → π* excitation are linear combinations of p carbon atoms orbitals. A closer inspection shows that π = p + p – p – p and π* = p – p – p + p (both coefficients multiplying orbitals and normalization factors are skipped for simplicity). The localized orbitals following from the π + π* and π – π* transformation, which would lead to eq if a two-configurational model were considered, are τ = p – p and τ = p – p. We conclude that X > 1 in the spatial regions surrounding all carbon nuclei, in agreement with Figure . For the excited hexatriene molecule, the CASSCF orbitals of interest π = p – p – p + p and π* = p – p + p – p, involved in the main single excitation, lead to τ = p – p and τ = p – p. According to eq , EDC should prevail either in the neighborhood of C1, C6 or C3, C4 nuclei. In contrast, in the vicinity of and nuclei the SDC effect (X < 1) dominates. This prediction is again confirmed in Figure . We conclude that ionic structures in the wave function are reflected in EDC, but for extended CAS functions describing ionic states, because of the presence of the covalent configurations next to the ionic ones in the wave function, the EDC effect competes with SDC. The overall effect of the latter may be locally greater than the former.

Having recognized that covalent and ionic states may be distinguished based on the role of the SDC and EDC regions and their importance in describing dynamic correlation energy, we introduce the ionicity index, IEDC, which measures the degree of ionicity of an ith state described with a wave function ΨThe index is then defined by the relative contribution of the region selected according to the X(r) > 1 criterion to the total correlation energy EcdDFTwhere X and the local quantities entering the correlation functional are computed for the wave function, Ψ. The definition is general, and any wave function model and correlation energy functional could, in principle, be employed in eqs and 6.

For small CAS(2,2) or CAS(4,4) functions, which capture the essential single excitation resulting in ionic character of the wave function, the introduced index is expected to be high. Including more active orbitals in CAS, as we have seen on the example of hexatriene, will result in obtaining regions in space where X < 1. This could hide away the EDC effect and result in a low value of the index, making it less useful in identifying ionic states. To avoid this problem, we propose to compute the X ratio defining the boundary of the integration in eq , as well as other local functions needed to compute IEDC after truncating the active space to only four most relevant natural orbitals, i.e. HOMO–1, HOMO, LUMO, and LUMO+1. Such effective projection on the CAS(4,4) wave function form, which in practice requires diagonalized one-electron density matrix with the natural orbitals ordered according to the occupation numbers, can be applied to any wave function (see the Supporting Information for details).

Extension of CASΠDFT for Excited States. The fact that ionic structure results in EDC, i.e. in appearance of spatial high-density regions where X(r) > 1, implies that special care must be taken of EDC in combined WFT-DFT methods. Recently, we have proposed a CASΠDFT method[19] which represents the total energy of the system as a sum of the CASSCF energy and the DFT-based dynamic correlation (cd) correctionThe DFT functional incorporates the Lee–Yang–Parr (LYP)[27] correlation functionaland is evaluated in a post-CAS fashion, i.e., proceeding a self-consistent CASSCF calculation. Therefore, the CASΠDFT method benefits from the fact that nondynamic correlation effects may be credibly represented already with small CAS in moderate basis sets. The key ingredient of CASΠDFT is the multiplicative P[X] function which adjusts the DFT correlation energy to account for SDC and EDC effects in spatial regions where its argument X(r), defined in eq , is smaller or greater than 1, respectively. In this work, we extend applicability of CASΠDFT to excited states by distinguishing between ground and excited states and introducing two dedicated correction functions Pgs[X] and Pexc[X], respectively. In the presented development both functions share the same formwhich assures independent representation of the SDC and EDC regions.

Expressions for both segments have been proposed in our earlier work.[19] For ground states, the b parameter in Pgs[X] has been fixed at bgs = ags – 1 to satisfy the PgsSDC[1] = 1 condition. The latter allows the DFT correlation to act at its full capacity when nondynamic correlation effects become negligible (X → 1). The single free parameter ags = 0.2 follows from a fit reproducing the H2 ground-state energy.[19] The form of PEDC[X] function was chosen to effectively account for the contribution from EDC regions in the 11Σ+ state of H2 upon bond dissociation.[19] Although the PEDC[X] segment of eq is more relevant for accurate description of the excited states rather than for the ground states, it should also be accounted for in the Pgs[X] correction. The parameter dgs is fixed at dgs = cgs – 1 to obey the continuity condition PgsEDC[1] = PgsSDC[1] and the parameters cgs = 2.6 and ggs = 1.5, successfully used in describing potential energy curves of molecules,[21,22] are retained.

In the Pexc[X] correction for excited states, both the asymptotic and continuity conditions at X = 1 are lifted. The parameters for the SDC and EDC segments were simultaneously fitted to reproduce the lowest π → π* excitations for four representative π conjugated molecules from the benchmark data set of ref (28): E-butadiene, cyclopentadiene, benzene, and all-E-hexatriene. The resulting parameters (aexc = 0.191856, bexc = 0.814889, cexc = 2.566480, dexc = 6.315147, and gexc = 1.501069) complete the definition of the correction function Pexc for excited states (note that they are close to their ground-state counterparts). The vertical excitation energy following from the proposed extended CASΠDFT method readswhere the indices 0 and i correspond to a ground state and an excited state of interest, respectively. Note that setting both P[X] correction functions to 1 reduces eq to evaluation of the excitation energy based on a simple addition of LYP correlation to the CASSCF energy, which we term the CAS+LYP method.

We validate the CASΠDFT performance for singlet vertical excitations predominantly of the π → π* type. We focus on a set of molecules comprising unsaturated aliphatic hydrocarbons as well as five- and six-membered aromatic hydrocarbons selected from the benchmark database of Schreiber et al.[28] (11 systems and 30 excitation energies in total). The CASΠDFT accuracy is compared against results from CASSCF, CAS+LYP, CASPT2,[29] and CC3[28] calculations. As a reference, we use the best theoretical estimates from the work of Schreiber et al.[28] All CASΠDFT calculations were performed with an in-house code. The 1- and 2-electron reduced density matrices, following from state-averaged CASSCF calculations, were acquired from the developer version of the Molpro[30] program. All systems were described with the TZVP basis set.[31] Geometries, active spaces, and the numbers of states in each symmetry used in the state-averaged (SA) CASSCF calculations (states of different spatial symmetry were state-averaged) are the same as in ref (28).

The analysis of individual excitation energies (Table ) and error statistics (Table ) reveals that the direct addition of LYP correlation in CAS+LYP brings virtually no improvement over regular CASSCF. Both methods exhibit mean absolute errors (MAE) which slightly exceed 1 eV and large standard deviations (SD) of 0.8 eV. In contrast, in CASΠDFT the LYP functional is modified with the introduced P[X] correction functions (eq ), which leads to a three-fold improvement with respect to CASSCF in terms of MAE (0.36 eV vs 1.05 eV, respectively) and reduction in standard deviation (from 0.8 eV in CASSCF to 0.5 eV in CASΠDFT). Compared to CASPT2 and CC3, CASΠDFT is slightly less reliable because of a larger spread of errors. Although the mean error (ME) of CASΠDFT amounts only to 0.08 eV, which matches the excellent CASPT2 result, the standard deviation in both CASPT2 and CC3 is close to 0.2 eV, whereas in CASΠDFT it reaches 0.5 eV. The accuracy of CASΠDFT calculations based on state-averaged CASSCF is especially good for the important lower excited states. For most of the two lowest excited states of 5-membered ring systems and linear polyenes, the errors stay below 0.4 eV. The 21A′ state of imidazole is an outlier reaching the error of 0.6 eV with respect to the reference (note that deviation for this state from the CC3 excitation energy is only 0.2 eV). For octatetraene, we observe larger deviations from the theoretical benchmark which amount to 0.6 and 0.8 eV for the 21A and 11B states, respectively. One should note, however, that in this case the disagreement between the best ab initio estimate[28] (4.45 eV for 21A and 4.7 eV for 11B) and experimental results (3.6 and 4.4 eV, respectively) is sizable.

Table 1

Vertical Excitation Energies (eV) and Index of Ionicity (in Percent)a

			CAS
molecule	state	IEDC	SCF	+LYP	ΠDFT	PT2	CC3	best est.
ethene	11B1u+	90	9.04	9.17	7.75	8.62	8.37	7.80
E-butadiene	11B1u+	80	7.88	7.90	6.51	6.42	6.58	6.18
	21Ag–	13	6.55	6.63	6.64	6.61	6.77	6.55
all-E-hexariene	11B1u+	54	7.30	7.31	5.08	5.35	5.58	5.10
	21Ag–	20	5.48	5.48	5.11	5.52	5.72	5.09
all-E-octatetraene	21Ag–	33	4.68	4.68	3.90	4.64	4.97	4.47
	11B1u+	55	6.69	6.69	3.84	4.70	4.94	4.66
cyclopropene	11B1	36	7.09	7.21	6.96	6.76	6.90	6.76
	11B2	64	8.11	8.17	6.73	7.06	7.10	7.06
cyclopentadiene	11B2+	75	7.03	7.07	5.41	5.52	5.73	5.55
	21A1–	28	6.54	6.60	6.70	6.48	6.61	6.31
	31A1+	84	10.46	10.51	9.14	8.39	8.69
norbornadiene	11A2	53	6.96	7.00	5.12	5.37	5.64	5.34
	11B2	46	8.60	8.64	7.21	6.12	6.49	6.11
	21B2	63	9.68	9.73	7.94	7.31	7.64
	21A2	36	9.79	9.84	8.10	7.42	7.71
benzene	11B2u−	22	4.83	4.83	4.48	5.04	5.07	5.08
	11B1u+	89	7.84	7.85	6.21	6.43	6.68	6.54
	11E1u+	65	9.23	9.24	7.36	7.09	7.45	7.13
	21E2g–	19	8.03	8.03	8.17	8.19	8.43	8.41
furan	11B2+	65	7.81	7.85	6.43	6.52	6.60	6.32
	21A1–	32	6.67	6.72	6.49	6.52	6.62	6.57
	31A1+	41	9.96	10.02	9.12	8.32	8.53	8.13
pyrrole	21A1–	31	6.53	6.57	6.26	6.30	6.40	6.37
	11B2+	65	7.65	7.69	6.78	6.33	6.71	6.57
	31A1+	44	9.41	9.45	8.84	8.06	8.17	7.91
imidazole	11A″	21	7.05	7.08	6.61	6.81	6.82	6.81
	21A′	32	6.76	6.79	6.81	6.58	6.58	6.19
	31A′	53	8.05	8.06	7.02	6.71	7.10	6.93
	21A″	26	8.49	8.53	8.16	7.90	7.93

All excitations are of the π → π* type except for excitations to the 11B1 state in cyclopropene (σ → π*) and to the 11A″ state in imizadole (n → π*). CASPT2 results are taken from ref (29); CC3 and “best est.” are from ref (28).

Table 2

Statistical Analysis of Singlet Excitation Energies for the Dataset of Table a

	CASSCF	CAS+LYP	CASΠDFT	CASPT2	CC3
all systems/states (24)
ME	1.06	1.10	0.08	0.09	0.27
MAE	1.12	1.15	0.37	0.17	0.27
SD	0.84	0.84	0.49	0.24	0.18
MAX	2.48	2.53	1.10	0.82	0.63
covalent states (9)
ME	0.11	0.14	–0.05	0.09	0.24
MAE	0.26	0.28	0.30	0.18	0.24
SD	0.30	0.30	0.40	0.22	0.23
MAX	0.57	0.60	0.62	0.43	0.63
ionic states (15)
ME	1.61	1.65	0.14	0.09	0.29
MAE	1.61	1.65	0.39	0.17	0.29
SD	0.44	0.43	0.53	0.25	0.14
MAX	2.48	2.53	1.10	0.82	0.57

Mean errors (ME), mean absolute errors (MAE), standard deviations (SD), and maximum signed errors (MAX), in eV, with respect to the “best est.” results of ref (28). States for which IEDC > 40% were qualified as “ionic”. Numbers in parentheses denote the cardinality of a given set.

Ionic States and the Index of Ionicity. Below we examine the novel index of ionicity, cf. eq , based on the EDC phenomenon and use it to verify the accuracy of the CASΠDFT method applied to states of ionic character. The values of IEDC presented in Table have been obtained with the ΠDFT model, ϵcDFT(r) = Pexc[X](r)ϵcLYP[ρ](r), the natural choice in this work. In the Supporting Information, we present the results obtained with LYP[27] and PBE[32] correlation functionals, ϵcDFT(r) = ϵcLYP[ρ](r) and ϵcDFT(r) = ϵcPBE[ρ](r), respectively, leading practically to the same ionicity index values as those obtained with ΠDFT, which manifests the universality of the index.

Excited states of alternant hydrocarbons may be grouped according to the minus and plus alternacy symmetry of the pertinent Hückel or Pariser–Parr–Pople (PPP) Hamiltonians.[33] In the VB picture, the singlet minus state has a covalent character while the plus state is ionic.[34−36] Gauging the ionicity at the IEDC value of 40% (in Table entries exceeding this value are marked in bold), one observes a striking agreement between the plus ionic states and the high value of the index. Another observation is that also for the remaining states classified as ionic according to the IEDC > 40% criterion, the CASSCF excitation energy error is greater than 1 eV (CASSCF excitation energies deviating by more than 1 eV from the reference are in bold in Table ). The poor performance of CASSCF for ionic excited states is also visible in Table ; the corresponding MAE amounts to 1.6 eV (for covalent excitations it is 0.3 eV). While adding the LYP correlation to CASSCF brings no improvement either for ionic or covalent excited states, the CASΠDFT method reduces the error of CASSCF for ionic excitations by a factor of 4 to an acceptable value of 0.4 eV. For covalent excitations, CASΠDFT retains the small MAE of CASSCF.

In summary, by analyzing wave functions of several representative excited states we have rigorously shown a direct relation between their ionic structure and regions of enhancement of dynamic correlation. Such regions are not meaningless peculiarities but dominate in the case of excited states of ionic character. This development has led us to extend the CASΠDFT method to excited states by developing a correction function designed for singlet excited states. The proposed extension has been applied to the π → π* excitations of organic molecules and proven reliable for the lowest excitations. The average absolute error is only 0.1 eV greater than that of the CC3 method, although the spread of errors from the latter method is 0.3 eV smaller than from CASΠDFT. It must be stressed that the computational effort to compute the ΠDFT correlation energy is negligible in comparison to the cost of the best-performing ab initio methods. This is because it involves solely local functions and, unlike in other WFT-DFT approaches, two-electron functions, for example two-electron integrals, are not involved. In practice, the cost of computing ΠDFT energy is as low as that of LYP energy (notice that ΠDFT requires construction of the on-top pair density in the active space, see the Supporting Information, which adds marginally to the overall cost).

Another achievement of this work is the introduction of the quantitative index of the wave function ionicity. While both the EDC-derived measure of ionicity and its counterpart based on the VB theory agree, the former identifies states as “ionic” only from the relative importance of the fundamental effects of suppressed and enhanced dynamic correlation. It is therefore straightforward to apply the ionicity index to any wave function, without the need to refer to the VB picture.

One can expect that the novel physical insight into the nature of ionic states, unraveled by understanding of the unknown features of the X ratio and their consequences for dynamic correlation, together with promising results obtained with the CASΠDFT method may spark further development of combined WFT-DFT methods dedicated to excited states.