Atreyee
Majumdar
and
Raghunathan
Ramakrishnan
*
Tata Institute of Fundamental Research, Hyderabad 500046, India. E-mail: ramakrishnan@tifrh.res.in
First published on 7th May 2024
We embark on a quest to identify small molecules in the chemical space that can potentially violate Hund's rule. Utilizing twelve TDDFT approximations and the ADC(2) many-body method, we report the energies of S1 and T1 excited states of 12880 closed-shell organic molecules within the bigQM7ω dataset with up to 7 CONF atoms. In this comprehensive dataset, none of the molecules, in their minimum energy geometry, exhibit a negative S1–T1 energy gap at the ADC(2) level while several molecules display values <0.1 eV. The spin-component-scaled double-hybrid method, SCS-PBE-QIDH, demonstrates the best agreement with ADC(2). Yet, at this level, a few molecules with a strained sp3-N center turn out as false-positives with the S1 state lower in energy than T1. We investigate a prototypical cage molecule with an energy gap <−0.2 eV, which a closer examination revealed as another false positive. We conclude that in the chemical space of small closed-shell organic molecules, it is possible to identify geometric and electronic structural features giving rise to S1–T1 degeneracy; still, there is no evidence of a negative gap. We share the dataset generated for this study as a module, to facilitate seamless molecular discovery through data mining.
In 2019, two independent studies confirmed STGs < 0 in the prototypical cycl[3.3.3]azines—cyclazine10 and heptazine11—using time-dependent density functional theory (TDDFT) approximations, as well as a few correlated wavefunction methods. Since then, there has been a renewed interest in exploring the historically significant inverse-STG candidates: N-containing triangular molecules12–20 and non-alternant hydrocarbons.21–23 Besides these classes of molecules, Bedogni et al.18 showed the possibility of STG < 0 in CnHnNn aza-rings.
Despite mounting computational evidence supporting the likelihood of Hund's rule violation, the credibility of negative STGs still attracts criticisms.24 This skepticism arises from the failure to account for the experimental conditions in computational modeling and the challenges posed by the large molecules for accurate ab initio calculations. More recently, Loos et al.25 confirmed negative STGs in triangulene systems through composite excited states modeling and provided theoretical best estimates (TBEs).
The present study aims to report STGs calculated using a many-body method and double-hybrid density functional theory (dh-DFT) models for 12880 small organic molecules with systematically varying compositions and structures. Using this data, we verify the possibility of Hund's rule violation in the chemical space. Our workflow for data generation and the design of the module pymoldis for querying the reported data is illustrated in Fig. 1. The rest of this article discusses the qualitative aspects of the data, their analysis, and the technical details of the calculations.
Fig. 1 Workflow outlining the design of the pymoldis module26 for data mining S1 and T1 energies across 12880 molecules with up to 7 CONF atoms in the bigQM7ω dataset.27 ADC(2) and TDDFT calculations were performed as a part of the present work. See (ESI†) for screenshots of example queries. |
Fig. 2 Cumulative distribution of S1–T1 energy gaps of 12880 bigQM7ω molecules calculated with various ab initio methods. The inset shows the distribution in the range of −0.2 to 0.2 eV. |
Fig. 2(a) shows the cumulative count of STGs of 12880 molecules predicted with ADC(2) and selected hybrid-DFT, long-range-corrected-DFT, and long-range-corrected-hybrid-DFT approximations. All the theoretical models featured in this plot predict positive STGs. Past studies17,28,29 have shown that explicit incorporation of electron correlation, for example, at the MP2-level (second-order many-body perturbation theory) as in dh-DFT, is a requirement to predict STG < 0. However, it has come as a surprise that at the ADC(2) level, which is the excited-state equivalent to MP2, none of the 12880 molecules show a negative STG. Similar results from dh-DFT methods, along with their spin-component-scaled (SCS) and opposite-spin-component-scaled (SOS), are shown in Fig. 2(b). The zoomed-in inset shows the distribution of values from RSX-QIDH to shift towards the positive domain compared to ADC(2) values. Upon SCS/SOS corrections,30 the distribution is shifted slightly to the negative domain. At the dh-DFT level, SCS-PBE-QIDH, an accurate method for modeling STGs of triangulenes, a few molecules exhibit STG < 0 eV; we give a detailed discussion of individual values later.
Fig. 3 illustrates the effect of SCS corrections to PBE-QIDH predicted values of STG in the form of probability densities of the shift in S1 and T1 energies with the inclusion of SCS. Overall, the SCS corrections lower the S1 energies while raising the T1 values, illustrating why SCS-PBE-QIDH favors smaller STGs than the unscaled method, PBE-QIDH. Fig. S1 of the ESI† displays similar plots for SOS-PBE-QIDH, SCS-RSX-QIDH, and SOS-RSX-QIDH. While the SCS/SOS variants of RSX-QIDH shift the S1 and T1 energies further apart compared to the PBE-QIDH variants, previous benchmark studies25,31 have shown SCS-PBE-QIDH to be more accurate for modeling molecules with negative STGs. Hence, one can conclude scaled-RSX-QIDH to result in more false-positive predictions (i.e. spurious predictions of STG < 0 eV) than scaled-PBE-QIDH. Notably, the SCS/SOS corrections applied to the PBE-QIDH and RSX-QIDH dh-DFT methods are specifically tailored for the TD-DFT framework.30 Hence, STGs predicted using the ΔSCF approach will be similar in unscaled and scaled DFT methods.
Fig. 3 Probability density of the shift in S1 and T1 energies (in eV) for 12880 molecules upon the inclusion of spin-component-scaling (SCS) in PBE-QIDH. |
The scatterplot in Fig. 4 offers a detailed view of the distribution of STGs around 0 eV. We find that the predominant entries shown in this plot are in the blue region denoted “true negatives in SCS-PBE-QIDH”, implying that the TDDFT method agrees with ADC(2), both predicting these molecules as Hund's rule obeying systems with positive STGs. Molecules shown in the red region of Fig. 4 are “false positives in SCS-PBE-QIDH” as these systems show STG > 0.04 eV according to ADC(2), while in the TDDFT formalism they show STG < 0 eV. As already highlighted in Fig. 2, there are no “true positives” or “false negatives” as none of the 12880 molecules in the bigQM7ω dataset exhibit a negative STG at the ADC(2) level.
Fig. 4 Scatterplot of S1–T1 energy gap ≤0.4 eV from SCS-PBE-QIDH and ADC(2) for the bigQM7ω molecules. |
The boundary separating positives and negatives in Fig. 4 is not sharp as both ADC(2) and SCS-PBE-QIDH have uncertainties >0.1 eV in their predictions (see benchmarks in the ESI†). Fig. S2 (ESI†) shows scatterplots of the joint distributions of STG with the S1 and T1 energies at both the SCS-PBE-QIDH and ADC(2) methods. At both levels of theories, one finds the small STG systems to have S1 and T1 energies in the range of 6–8 eV. Additionally, we find a molecule with S1 and T1 energies in the 3–4 eV range to have a small STG at the SCS-PBE-QIDH level (see Fig. S2c and d, ESI†). By querying the dataset using the pymoldis module (see Fig. S3–S12, ESI†), we found the corresponding molecule to be 2,6-dihydro-1H-pyridin-3-one (SMILES: OC1CNCCC1), which is cyclohexenone with an N atom at the δ-position.
To identify geometrical aspects common to the small STG molecules, we have queried the pymoldis database for ten molecules with the smallest STG according to TDA/SCS-PBE-QIDH. Fig. 5 shows a screenshot of the query. A geometric moiety common to all N-containing systems is a substantial deviation of the bonding environment of the sp3-N center from the ideal geometry. For the ten benchmark triangular systems, the S1 and T1 excitation energies are <3 eV due to the possibility of a low-energy n → π* transition. However, in the systems shown in Fig. 5, the S1/T1 excitations are of the n → σ* type with excitation energies >6 eV. Nine of the ten molecules shown in Fig. 5 contain a 3-membered heterocyclic ring, while one is a fluorinated cyclopropane derivative; the S1 and T1 excitation energies of the latter are >8 eV, at both TDDFT and ADC(2) levels while the corresponding STG = −0.01 eV and 0.13 eV at TDDFT and ADC(2) levels, respectively. The molecule corresponding to SMILES, CC1C2CCN1C2, in Fig. 5 contains a propellane-type cage with STG = −0.01 eV at the TDDFT level and an STG of 0.06 eV at the ADC(2) level. The bigQM7ω dataset comprises several such cage systems with small STGs.
Fig. 5 Data-mining the bigQM7ω dataset in pymoldis26 to identify molecules with small S1–T1 gaps. For the ten molecules with the lowest gaps according to TDA-SCS-PBE-QIDH, SMILES strings and S1/T1 energies are displayed. The corresponding ADC(2) values are also given alongside. Molecular structures are displayed in cartoon format below. See the ESI† for more examples of data-mining exercises. |
Starting with a cage-type molecule with a strained N, we have explored the possibility to design a molecule with the character of the lowest excitation as n → π*. We started with a symmetric cage system, quinuclidine, a [2.2.2]propellane with an axial CH group replaced by an N atom, structure 1 in Fig. 6. This molecule has an STG of 0.1 eV at the TDA-SCS-PBE-QIDH/def2-TZVP level. This small gap is because the excitation is primarily from non-bonding MO of N (highest occupied MO, HOMO) to C–H σ* of the cage (lowest unoccupied MO, LUMO); their corresponding densities show poor overlap leading to a small value of exchange integral, Kar. Further, to arrive at a local-geometric environment of the N atom as in cycl[3.3.3]azines, we have introduced an additional cage to constrain the N-center to a plane (structure 2 in Fig. 6). This structure comprises perfectly co-planar C–N bonds resulting in degenerate S1–T1 levels. We have also modified quinuclidine by attaching three ethylene groups (structure 3). In this structure, the S0 → S1 excitation has the character n → π* (MO indices, n: 52, which is the HOMO, and π*: 53–55) while the S0 → T1 excitation has the character π → π* (MO indices, π: 49–51) with a large STG of 1.5 eV. Finally, we combined structural modifications introduced in structure 2 and structure 3 to arrive at structure 4 with a planar N interacting with π moieties through space. Interestingly, this system resulted in an STG of −0.21 eV at the DFT level. To verify this prediction of a negative STG, we have also performed ADC(2)/def2-TZVP calculations. For both structure 1 and structure 4, we find the magnitudes of the ADC(2) excitation energies to be lower than the DFT values. While both energies drop by a similar magnitude in the former, in structure 4, the energy of T1 drops more than the S1 energy, giving rise to a nearly zero STG at the ADC(2) level. This case study indicates that for molecules such as structure 4, even some of the best double-hybrid DFT methods can spuriously predict a negative STG; hence, one may consider many-body methods such as ADC(2) as a baseline theory.
While structure 4 is a minimum on the potential energy surface as verified through vibrational frequency analysis, we do not expect the system to be relevant to the thermally activated delayed fluorescence (TADF) applications.10,13,32–34 On the other hand, it is a compelling computational chemistry exercise to modulate a molecule's STG by chemical modifications. Hence, even though structure 4 seems to be yet another false positive in the search for a Hund's rule-violating molecule, we have examined it further. We inspected the shape of the MOs involved in the S1 and T1 excitations and found the excitations to be primarily HOMO→LUMO type. These MOs are on display in Fig. 7, from which we visually conclude that the densities of HOMOs and LUMOs do not overlap. For a more quantitative analysis, we calculated the Λ-index35 defined as using Multiwfn,36 and obtained the values: 0.37 and 0.40 for the S1 and the T1 states, respectively. The Λ-index quantifies the degree of overlap between hole and electron in S0 → S1 and S0 → T1 excitations. In comparison, for both cyclazine and heptazine, the values of Λ for the S1 and T1 states are nearly 0.50.
Fig. 7 Plots of HOMO and LUMO of the ethylene substituted double cage derivative of quinuclidine, structure 4, shown in Fig. 6. |
For all 12880 molecules in bigQM7ω, we performed single-point vertical excited state calculations of the S1 and T1 energies using 12 DFT methods: PBE0,43 B3LYP,44 CAM-B3LYP,45 ωB97X-D3,46 LC-BLYP,47 LC-PBE,48 PBE-QIDH,49 SCS-PBE-QIDH, SOS-PBE-QIDH,30 RSX-QIDH,50 SCS-RSX-QIDH30 and SOS-RSX-QIDH.30 We also calculated S1 and T1 energies using the correlated excited state method: second-order algebraic diagrammatic construction, ADC(2). The accuracy of ADC(2) and SCS-PBE-QIDH in combination with other settings is evaluated in the ESI† using previously reported25 TBEs of 10 triangular systems as references. For these systems and the cage-type molecules, we performed geometry optimization – using the ωB97X-D3 DFT method with tightscf and tightopt keywords in combination with the def2-TZVP basis set. Minimum energy structures and S1/T1 energies of the bigQM7ω molecules can be queried using the pymoldis module presented in this study (see Fig. S8 and S12 in the ESI†).
We performed ADC(2) calculations using QChem 6.0.2 and DFT calculations using ORCA 5.0.4.51,52 With in TDDFT, we calculated twelve energy eigenvalues—six singlets and six triplets—which we sorted separately to extract S1 and T1 excitation energies. In all calculations, we used the resolution-of-identity (RI) approximation.53,54 In DFT calculations, we used the ‘chain-of-spheres’ (COS) algorithm for exchange integrals (RIJCOSX). In dh-DFT calculations, we used the universal fitting auxiliary basis sets by Weigend55 (denoted def2/J) along with the def2-TZVP/C and aug-cc-pVTZ/C basis sets for the orbital basis sets def2-TZVP, and aug-cc-pVTZ, respectively.
The data presented in this study is importable in Python code for data mining endeavors. Using this infrastructure, we identified molecules with vanishing STGs, some of which have negative values at the SCS-PBE-QIDH level while not violating Hund's rule as per ADC(2) predictions. A common geometric feature of these molecules was a substantial deviation of an N-atom from the typical sp3 environment with both singlet and triplet excitations showing the n → σ* character and potentially vanishing exchange interaction integral between the MOs involved in excitation.
We have selected a cage structure and attached ethylene groups to mimic the environment of the N atom as in the well-known cases of cycl[3.3.3]azines. The corresponding MOs involved in S1/T1 excitations exhibit characteristics seen in previously studied triangular negative-STG systems. Upon further scrutiny, we showed this molecule to obey Hund's rule. Yet, introducing a polarizable environment in this system through donor–acceptor groups may selectively stabilize S1 over T1.10,23 In this study, we did not investigate the practical utility of the small molecules studied here in the context of TADF. Such exploration necessitates meticulous consideration of adiabatic effects on a case-by-case basis, a task that exceeds the scope of our present investigation. The present study demonstrates that a data-driven approach allows for gaining insight into the molecular structural factors that can quench the singlet–triplet energy gap. We offered evidence that the chemical space of small closed-shell organic molecules lacks geometric and electronic structural factors that may be necessary for a negative S1–T1 energy gap.
Theoretical studies have identified only a few molecular fingerprints to favor negative STGs. Dynamic spin polarization, attributed to double excitation effects involving frontier orbitals, has emerged as a potential mechanism to induce a negative STG.1,57,58 While a quantitative relationship exists between molecular structural features and zero STG, a corresponding structure–property relation for negative STG remains elusive. Introducing functional groups is one promising avenue for designing large synthetically tractable molecules with negative STGs. For theoretical explorations, our research highlights the limitations of using DFT methods, which can result in false positives. Consequently, there is a pressing need for efficient strategies to accelerate predictions using correlated wavefunction methods. For large-scale investigations, data-driven approaches such as machine learning can complement first-principles modeling, especially when combined with inverse-design strategies, such as genetic algorithms.59,60
Footnote |
† Electronic supplementary information (ESI) available: Contains the following: assessment of S1–T1 gaps from ADC(2) and SCS-PBE-QIDH for 10 triangular molecules. Table S1 compares S1/T1 energetics calculated using the ADC(2) method with theoretical best estimates. Table S2 compares the S1/T1 energetics predicted by TDDFT and TDA. Table S3 contains ADC(2) and TDA energies calculated using DFT-level geometries. Table S4 provides various error metrics for ADC(2) and SCS-PBE-QIDH predicted S1/T1 energetics. Fig. S1 illustrates the shifts in S1 and T1 due to SCS/SOS corrections to PBE-QIDH and RSX-QIDH methods. Fig. S2 shows a scatterplot of S1 and T1 energies with S1–T1 gaps. Fig. S3–S12 offer screenshots of data mining exercises. Minimum energy geometries of structures 1–4 in Fig. 6 are also listed. Sample Python notebooks and further details are available at https://github.com/moldis-group/pymoldis.26 See DOI: https://doi.org/10.1039/d4cp00886c |
This journal is © the Owner Societies 2024 |