Excited-State Polarizabilities: A Combined Density Functional Theory and Information-Theoretic Approach Study

1. Introduction

Molecular polarizability, especially the static electric dipole polarizability (α), is a fundamental physicochemical property. It reflects the change of a molecule’s dipole moment in a linear-response manner, as resulted from an external electric field perturbation. [1] The experimental determination of excited-state electro-static properties are mainly the Stark spectroscopy or electronic absorption/emission [2,3] method and the flash photolysis time-resolved microwave-conductivity (FP-TMRC) [4,5] technique.

In classical physics, the polarizability can be approximately obtained in terms of the volume of a system. [6,7] For example, many strong correlations have been observed for both atoms and molecules. [8,9,10,11,12,13,14,15,16] It is worthwhile to mention that Tkatchenko and Scheffler (TS) [17] proposed to use atomic volumes and atomic polarizabilities to predict the ground-state polarizabilities for small molecules. Recent progress can be found in ref 18. However, its performance for excited-state systems has not been reported.

In quantum mechanics, the polarizability can be obtained by iteratively solving the coupled-perturbed Hartree–Fock (CPHF) equation [19,20] or its Kohn–Sham DFT (density functional theory) counterpart. [21] Of note, this protocol requires a sufficiently large basis set with polarization and diffuse functions and huge computational costs. Note that the computational barriers can be partly overcome by using some linear-scaling methods. [22,23,24] In addition, machine learning (ML)-based [25,26,27] methods and a regression-based [28] model have been applied to predict the S₀ polarizabilities. It is worthy to note that the polarizability can be related to the band gap of HOMO (highest occupied molecular orbital) and LUMO (lowest unoccupied molecular orbital) in an inverse manner. [29,30,31]

In the literature, only a few studies [32–38] have been reported for the excited-state polarizabilities. This is likely because that accurate predictions of excited-state geometries and molecular properties of large molecules is a tough nut to crack, especially when there are perturbations such as external fields.

Following our previous work where the information-theoretic approach (ITA) quantities are employed to predict the S₀ polarizabilities of both small and large molecules, [39,40] here we aim to predict the S₁ polarizabilities of 2-(2′-hydroxyphenyl)benzimidazole (HBI, 1) and its derivatives as shown in Figure 1. For 1, it is well-documented [41] that the S₀ (Figure 1a) intramolecular proton transfer (IPT) reaction is difficult to take place and the S₁ (Figure 1b) or T₁ (triplet, not shown) intramolecular proton transfer (ESIPT) process can easily happens. Thus, in this work, only S₀ and S₁ are considered to reduce the computational cost without compromising the results much. We have found that with the S₀ or S₁ electron densities as input, ITA quantities can be in good correlations with the excited-state polarizabilities. When the transition densities are considered, both S₀ and S₁ polarizabilities can be in good relationships with ITA quantities. Furthermore, excitation and emission energies can be predicted based on multiple linear regression equations of ITA quantities. For the first time, we have applied the ITA quantities to predict the excited-state polarizabilities. It is anticipated that this protocol can be readily applied to condensed-phase systems.

2. Results

Shown in Table 1 are the correlation coefficients (R²) between the S₀ polarizabilities (α_iso) and ITA quantities, molecular volumes, and quadrupole moments, which are obtained at the CAM-B3LYP/6-311+G(d) level. It is clear from Table 1 that Ghosh−Berkowitz−Parr (GBP) entropy (S_GBP), 2nd and 3rd relative Rényi entropy (^rR₂ and ^rR₃), information gain (I_G), G₁, G₂, and G₃, and quadrupole moments (Θ_iso) are in strong linear relationships with α_iso, with R² > 0.8. However, molecular volumes are in only reasonaly good coorelation with α_iso, with R² = 0.618, indicating that it is not a good descriptor of α_iso. Of note, the G₃ data have been shown to be strongly correlated with α_iso for various systems, among which are 30 planar or quasi-planar bases, [39] 20/40/8000 amino acids/dipeptides/tripeptides, [39] and so on. [40] Also, the Θ_iso values can be in good relationship with α_iso and its theoretical rational can be found in ref 40. However, a solid and sound theoretical verification between G₃ and α_iso is staill lacking. Overall, the strong correlations can serve as an argument that our computational results are convincing.

Collected in Table 2 are the S₁ polarizabilities, the S₁ ITA quantities including Shannon entropy (S_S), Fisher information (I_F), 2nd and 3rd relative Rényi entropy (^rR₂ and ^rR₃), G₂ and G₃, molecular volumes (Vol), and quadrupole moments (Θ_iso), which are obtained at the TD-CAM-B3LYP/6-311+G(d) level. Also given in Table 2 are the correlation coefficients (R²) between the S₁ polarizabilities and other quantities at S₁. Note that some ITA quantities, well-defined at S₀, are numerically ill-behaved at S₁ and thus missing. One can see that ^rR₂, ^rR₃, G₂, G₃, Vol, and Θ_iso are in good correlations with α_iso at S₁, with R² > 0.8. It is intriguing to note that at S₁, G₃ is still in good correlation with α_iso. This is the first time to observe such a phenomenon. However, admittedly, the theoretical foundation lags behind the numerical evidence introduced in this work. Moreover, we have found that for Vol, the correlation coefficient is much stronger at S₁ (0.908) than that at S₀ (0.618). One possible reason is that the excited-state relaxation expands the volume and polarizability space. Finally, one can discover in Column 4 of I_F, 2, 8, and 13 seem to have abnormal values compared with the others. They are all Br-containing, indicating that there may be some regions for heavy atoms where density gradients are numerically ill-behaved at S₁. Similar results can also be observed for I_F at S₀ (see Table 1). Overall, we have unraveled that excited-state densities and molecular properties are mutually entangled.

Now we have shown that ITA quantities can be correlated with α_iso either at S₀ or S₁. It is natural to ask if one can use ITA quantities at S₀ to predict α_iso at S₁. The answer is definitely yes! In Table 3, We have tabulated the correlation coefficients (R²) between the α_iso values at S₁ and ITA quantities, molecular volumes, and quadrupole moments at S₀ as introduced in Table 1. More details can be found in Table S1. Except Shannon entropy (S_S), Fisher information (I_F), and molecular volumes (Vol), the other quantities at S₀ are in good relationships with α_iso at S₁, with R² > 0.8. Moving forward, we ask if one can use the transition density (matrix) as input for ITA quantities to correlate with α_iso either at S₀ or S₁. The answer is again yes! Shown in Table 4 are the strong correlations between α_iso either at S₀ or S₁ and ITA quantities with the transition density matrix. Except Shannon entropy (S_S) and Fisher information (I_F), the other ITA quantities are in strong correlations with α_iso either at S₀ or S₁, with R² > 0.8. Implication of this part is straightforward that electron-density-based quantities can be used to predict the excited-state properties, such as molecular polarizabilities.

Next, we will compare the α_iso data (either at S₀ or S₁) predicted by the TS formulas as with conventional results as reference. Employing the original Tkatchenko–Scheffler (TS) formula [17] on top of Becke [42] or Hirshfeld [43] partitions, the α_iso data (either at S₀/S₁) are either strongly underestimated or overestimated, with MUE(%) up to –24.90/–21.21 and 6.62/10.82, respectively, as shown in Table 5. It is found that a mean value can reduce the MSE(%) to 16.00/15.23. Moreover, with the new TS formula, [18] the results are not improved but worsened as shown in Table 6. Taken together, we have found that the TS formulas have large room to improve in predicting the S₁ polarizabilities.

Finally, we have found that both excitation and emission energies can be predicted on top of multiple linear regression equations of ITA quantities. For example, one can use the transition density matrix as input for ITA quantities to correlate with the excitation energies. Similarly, if the S₁ densities are used for ITA quantities, the emission energies can be predicted. Based on the two regression equations,

λ_fit = 0.32*S_S + 0.0027*I_F ‒ 0.31*S_GBP ‒ 0.85*^rR₂ + 2.87*^rR₃ ‒ 27.99*I_G + 0.33*G₁ ‒ 0.050*G₂ + 0.0099*G₃,

and

λ_fit = 0.085*S_S + 0.00088*I_F + 1.36*^rR₂ ‒ 1.42*^rR₃ ‒ 0.080*G₂ + 0.038*G₃,

we have obtained that the MUEs (mean unsigned error) and the MSEs (mean signed error) are –0.04/0.20 eV and 0.00/0.22 eV for excitation and emission energies, respectively. This indicates that the inaccuracy of this protocol is comparable to that of underlying approximations of DFT. [44,45].

3. Discussion

To accurately and efficiently predict the excited-state polarizabilities is an ongoing issue. Solving standard CPHF/CPKS equations are computationally intensive and the computational costs can be intractable for macromolecular systems. Other algorithms and models available in the literature are normally concerned with the ground-state polarizabilities. Within this context; we proposed to apply some density-based ITA quantities to correlate with α_iso at S₀/S₁. This is inspired by our previous work on predicting the ground-state polarizabilities for small and macromolecular systems. Our tentative results have shown that the protocol should be a promising theoretical tool. More systems along this line need to be considered to make this protocol more robust and applicable. We have to point out when the system under study becomes larger and larger; the molecular wavefunctions (thus electron density) are a tough nut to crack; sometimes computationally intractable. Under these circumstances; we have to resort to linear-scaling electronic structure methods; such as GEBF (generalized energy-based fragmentation method), [46–49] where only small subsystems of a few atoms or groups are treated.

Next, we will look into the TS method, as mentioned previously. We have already found that based upon the Hirshfeld or Becke partition scheme, the original TS formula has an unsatisfactory performance either by overestimating or underestimating the S₁ polarizabilities. Apparent reduction of the deviations can be obtained by averaging the two sets of results. The reason behind is unclear at the moment. From the original formula,

$${\alpha }_{\mathrm{m}\mathrm{o}\mathrm{l}}^{\mathrm{T}\mathrm{S}-\mathrm{o}\mathrm{l}\mathrm{d}}=\sum _{\mathrm{A}}{\alpha }_{\mathrm{A}}^{\mathrm{e}\mathrm{f}\mathrm{f}}=\sum _{\mathrm{A}}{\alpha }_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}\left(\frac{V_{\mathrm{A}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}{V_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}\right)$$

one can easily argue that the weights $\left(\frac{V_{\mathrm{A}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}{V_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}\right)$ may be the root cause of its poor performance, mainly because the atomic polarizabilities ${\alpha }_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}$are experimentally determined and computationally verified, as summarized in ref 50. In the same spirit, a revised TS formula,

$${\alpha }_{\mathrm{m}\mathrm{o}\mathrm{l}}^{\mathrm{T}\mathrm{S}-\mathrm{n}\mathrm{e}\mathrm{w}}=\sum _{\mathrm{A}}{\alpha }_{\mathrm{A}}^{\mathrm{e}\mathrm{f}\mathrm{f}}=\sum _{\mathrm{A}}{\alpha }_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}{\left(\frac{V_{\mathrm{A}}^{\mathrm{e}\mathrm{f}\mathrm{f}}}{V_{\mathrm{A}}^{\mathrm{f}\mathrm{r}\mathrm{e}\mathrm{e}}}\right)}^{4/3}$$

has witnessed improved performance of predicting the S₀ polarizabilities. However, its predicting power has been shown to be far from satisfactory for macromolecules. In this work, we further corroborate that both the original and new TS formulas fail to give a satisfactory description of the S₁ polarizability. This indicates that the two formulas are oversimplified and may be system-dependent. Overall, the volume-based are inferior to the density-based ITA quantities.

4. Materials and Methods

4.1. Information-Theoretic Approach Quantities

Shannon entropy S_S [51] and Fisher information I_F [52] are two cornerstone quantities in information theory. They are defined as Equations (1) and (2), respectively.

$$S_{\mathrm{S}}=-\int \rho \left(\mathbf{r}\right)\mathrm{l}\mathrm{n}\rho \left(\mathbf{r}\right)d\mathbf{r}$$

(1)

$$I_{\mathrm{F}}=\int \frac{{|\nabla \rho (\mathbf{r}\left)\right|}^2}{\rho \left(\mathbf{r}\right)}d\mathbf{r}$$

(2)

where $\rho \left(\mathbf{r}\right)$ is the electron density and $\nabla \rho \left(\mathbf{r}\right)$is the density gradient. The physical picture of $S_{\mathrm{S}}$and $I_{\mathrm{F}}$ is clear; the former measures the spatial delocalization of the electron density and the latter gauges the sharpness or localization of the same.

Except the total density, more ingredients, such as kinetic-energy density, can be used to defined an ITA quantity. With both the electron density and the kinetic energy density, Ghosh, Berkowitz, and Parr developed a formula for entropy (S_GBP), [53]

$$S_{\mathrm{G}\mathrm{B}\mathrm{P}}=-\int \frac{3}{2}k\rho \left(\mathbf{r}\right)\left[c+\mathrm{l}\mathrm{n}\frac{t(\mathbf{r};\rho )}{t_{\mathrm{T}\mathrm{F}}(\mathbf{r};\rho )}\right]d\mathbf{r}$$

(3)

where t(r; ρ) and t_TF(r; ρ) represent the non-interacting and Thomas–Fermi (TF) kinetic energy density, respectively. Here k, c, and c_K are three constants [k, the Boltzmann constant, c = (5/3) +ln(4πc_K/3), and c_K = (3/10)(3π²)^2/3]. Full integration of the kinetic energy density t(r; ρ) leads to the total kinetic energy T_S via

$$\int t\left(\mathbf{r};\rho \right)d\mathbf{r}=T_{\mathrm{S}}$$

(4)

while t(r; ρ) can be obtained from the canonical orbital densities,

$$t\left(\mathbf{r};\rho \right)=\sum _i\frac{1}{8}\frac{\nabla {\rho }_i·\nabla {\rho }_i}{{\rho }_i}-\frac{1}{8}{\nabla }^2\rho$$

(5)

and t_TF(r; ρ) is simply cast in terms of $\rho \left(\mathbf{r}\right)$,

$$t_{\mathrm{T}\mathrm{F}}(\mathbf{r};\rho )=c_{\mathrm{K}}{\rho }^{5/3}\left(\mathbf{r}\right)$$

(6)

Of note, the kinetic-energy density can differ in its form, thus can be used in different contexts. [54,55,56,57,58,59,60,61] But, S_GBP satisfies the maximum-entropy requirement from a mathematical viewpoint. [53]

Moving forward, some ITA quantities have been introduced for chemical reactions. As new reactivity descriptors in conceptual density functional theory (CDFT), [62,63,64,65] one example is relative Rényi entropy [66] of order n

$$R_n^r=\frac{1}{n-1}\mathrm{l}\mathrm{n}\left[\int \frac{{\rho }^n\left(\mathbf{r}\right)}{{\rho }_0^{n-1}\left(\mathbf{r}\right)}d\mathbf{r}\right]$$

(7)

Information gain [67] (also called Kullback−Leibler divergence or relative Shannon entropy) I_G is expressed as follows,

$$I_{\mathrm{G}}=\int \rho \left(\mathbf{r}\right)\mathrm{l}\mathrm{n}\frac{\rho \left(\mathbf{r}\right)}{{\rho }_0\left(\mathbf{r}\right)}d\mathbf{r}$$

(8)

In Equations (7) and (8), ρ₀(r) and ρ(r) satisfy the same normalization condition, and ρ₀(r) denotes the reference-state density.

Recently, [68] one of the present authors proposed three functions G₁, G₂, and G₃, at both atomic and molecular levels. They are defined as below:

$$G_3=\sum _{\mathrm{A}}\int {\rho }_{\mathrm{A}}\left(\mathbf{r}\right){[\nabla \mathrm{l}\mathrm{n}\frac{{\rho }_{\mathrm{A}}\left(\mathbf{r}\right)}{{\rho }_{\mathrm{A}}^0\left(\mathbf{r}\right)}]}^{2}d\mathbf{r}$$

(9)

$$G_1=\sum _{\mathrm{A}}\int {{\nabla }^2\rho }_{\mathrm{A}}\left(\mathbf{r}\right)\frac{{\rho }_{\mathrm{A}}\left(\mathbf{r}\right)}{{\rho }_{\mathrm{A}}^0\left(\mathbf{r}\right)}d\mathbf{r}$$

(10)

$$G_2=\sum _{\mathrm{A}}\int {\rho }_{\mathrm{A}}\left(\mathbf{r}\right)\left[\frac{{{\nabla }^2\rho }_{\mathrm{A}}\left(\mathbf{r}\right)}{{\rho }_{\mathrm{A}}\left(\mathbf{r}\right)}-\frac{{\nabla }^2{\rho }_{\mathrm{A}}^0\left(\mathbf{r}\right)}{{\rho }_{\mathrm{A}}^0\left(\mathbf{r}\right)}\right]d\mathbf{r}$$

(11)

Equations (9)–(11) have been theoretically derived, numerically verified, and have witnessed many applications, as can be found in refs 39, 40, 68, and 69. It is one of our major achievements during the past decade, when we aimed to glue the density functional theory and information theory together in a seamless manner. Because these two theories both can have as electron density as input. Our recent progress along this line can be found in two reviews. [70,71] As another prominent example, we have applied the ITA to appreciate homochirality. [72,73]

Finally, the Hirshfeld’s stockholder approach [43,74–77] is often introduced to partition atoms in a molecule in the literature, as defined in Equation (12),

$${\rho }_{\mathrm{A}}\left(\mathbf{r}\right)={\omega }_{\mathrm{A}}\left(\mathbf{r}\right)\rho \left(\mathbf{r}\right)=\frac{{\rho }_{\mathrm{A}}^0\left(\mathbf{r}\right)(\mathbf{r}-{\mathbf{R}}_{\mathrm{A}})}{\sum _{\mathrm{B}}{\rho }_{\mathrm{B}}^0\left(\mathbf{r}-{\mathbf{R}}_{\mathrm{A}}\right)}\rho \left(\mathbf{r}\right)$$

(12)

Here,${ \rho }_{\mathrm{A}}\left(\mathbf{r}\right)$ is the atomic Hirshfeld density, ${\omega }_{\mathrm{A}}\left(\mathbf{r}\right)$ is a sharing function, ${\rho }_{\mathrm{B}}^0\left(\mathbf{r}-{\mathbf{R}}_{\mathrm{A}}\right)$ is the atomic density of B centered at R_A. The sum over all the free atom densities, typically spherically averaged S₀ atomic densities, is normally termed the promolecular density. The Stockholder approach is natural in the context of ITA because it is also based on information-theoretic arguments. Alternative partitioning schemes include Becke’s fuzzy atom approach [44] and Bader’s zero-flux atoms-in-molecules (AIM) method. [78]

4.2. Computational Details

All density functional theory (DFT) calculations were performed with the Gaussian 16 [79] package. Default options include ultrafine integration grids and tight self-consistent field convergence, which are adopted to eliminate numerical noises. The ground- and excited-state structural relaxation was fully carried out at the CAM-B3LYP/6-311+G(d) [80,81] and TD-CAM-B3LYP/6-311+G(d) [80,81,82–84] level, respectively, for 27 molecular systems, as shown in Figure 1. The optimized atomic Cartesian coordinates are supplied in the Supplementary Materials. Subsequent harmonic vibrational frequency calculations were executed at the same level and no imaginary frequencies were observed by direct visual inspection. The isotropic polarizabilities [α_iso = (α_xx + α_yy + α_zz)/3] and isotropic quadrupole moments [Θ_iso = (Θ_xx + Θ_yy + Θ_zz)/3], molecular volumes (at 0.001 e/Bohr³ contour surface of electronic density), and molecular wavefunctions were obtained at the CAM-B3LYP/6-311+G(d) level. The Multiwfn 3.8 [85] program was utilized to calculate all ITA quantities at S₀ and S₁ by using the checkpoint or wavefunction file as the input. The stockholder Hirshfeld partition scheme of atoms in molecules was employed when atomic contributions were concerned. The reference-state density was the neutral atom calculated at the same level of theory as molecules.