Influence of Pseudo-Jahn–Teller Activity on the Singlet-Triplet Gap of Azaphenalenes

Abstract

We probe the sensitivity of the singlet-triplet energy gap of selected azaphenalenes to symmetry lowering induced by Jahn–Teller interactions. While cyclazine in its characteristic $D_{\rm 3h}$ structure defies Hund’s rule, CCSD(T)-level modeling suggests its structure corresponds to two equivalent minima of $C_{\rm 3h}$ symmetry undergoing rapid automerization. The combined effect of symmetry reduction and high-level corrections indicates a negligible singlet-triplet gap in cyclazine. Notably, pentazine and heptazine prefer symmetric structures exhibiting negative gaps in accord with experiments. Azaphenalenes containing nitrogen atoms at electron-deficient sites exhibit stronger in-plane structural distortion; in their low-symmetry energy minima they adhere to Hund’s rule.

I Introduction

Azaphenalene (1AP in FIG. 1), also known as mono-azaphenalene, cycl[3.3.3]azine, or cyclazine, and its poly-aza analogues, wherein nitrogen (N) atoms replace the peripheral CH groups, display a potential to violate the Hund’s ruleKutzelnigg and Morgan (1996)—their first excited singlet state, S₁, is suggested to be lower in energy than the triplet counterpart, T₁de Silva (2019); Ehrmaier et al. (2019); Pollice et al. (2021); Aizawa et al. (2022); Li et al. (2022); Ricci, Sancho-García, and Olivier (2022); Ghosh and Bhattacharyya (2022); Tučková et al. (2022); Bedogni et al. (2023, 2023); Won, Nakayama, and Aizawa (2023); Drwal et al. (2023); Loos, Lipparini, and Jacquemin (2023); Pollice, Ding, and Aspuru-Guzik (2024); Garner, Blaskovits, and Corminboeuf (2024); Valverde et al. (2024); Pollice, Ding, and Aspuru-Guzik (2024). Three-fold symmetric 1APCunningham et al. (1969) and heptazine (7AP)Hosmane, Rossman, and Leonard (1982); Shahbaz et al. (1984) were among the first APs synthesized, followed by the two-fold symmetric, pentazine (5AP), along with a few other asymmetric APsRossman et al. (1985); see FIG. 1 for the structures of 5AP and 7AP. GimarcGimarc (1983) ascribed the stability of 7AP to the substitution of the CH moiety by N at the sites of 1AP with significant densities of the highest occupied molecular orbital (HOMO), as shown in FIG. 1. This topological charge stabilization effectGimarc (1983) rationalized the synthesis of another three-fold symmetric APLeupin et al. (1986). A similar effect diminishes anti-aromaticity through double-bond localizations in substituted indaceneHeilbronner and Yang (1987) and pentaleneTerence Blaskovits, Garner, and Corminboeuf (2023); Omar et al. (2023); Meiszter et al. (2024) that, along with APs, are among the too few molecular cores identified so far to exhibit negative STGs. A search for other examples in the chemical space of small organic molecules showed no evidence of STG inversionMajumdar and Ramakrishnan (2024).

Figure 1: Representation of APs explored in the present study, and density contours of the frontier MOs of 1AP.

Leupin et al. were the first to suggest the likelihood of a negative S${}_{1}-$T₁ gap (STG) in 7APLeupin and Wirz (1980). Its photochemistry exhibits delayed fluorescence, characteristic of a small, negative STGEhrmaier et al. (2019). Further, fluorescence measurements suggested an STG of $-$0.011 eV for a substituted 7APAizawa et al. (2022). More recently, an absorption spectrum confirmed the STG of 5AP in argon matrix to be $-0.047\pm 0.007$ eVWilson et al. (2024) aligning with independent measurements on dialkylamine-substituted 5APKusakabe et al. (2024). Inverted STG promises the scope for designing the next generation organic light emitting diodes (OLEDs), with efficiencies not limited by spin statisticsAizawa et al. (2022); Li et al. (2022); Won, Nakayama, and Aizawa (2023).

In this study, we investigate structural preferences exhibited by APs shown in FIG. 1 and explore how symmetry lowering, driven by second-order Jahn–Teller (SOJT) interactions, impacts their STGs. We first focus on 1AP, 5AP, and 7AP for which Loos et al. Loos, Lipparini, and Jacquemin (2023) performed detailed calculations of STGs and provided the theoretical best estimates (TBEs) as $-0.131$ eV, $-0.119$ eV, and $-0.219$ eV, respectively. We also examine the three APs—2AP, 3AP, and 4AP—containing N at electron-deficient sites of 1AP (see FIG. 1). Though these molecules lack experimental measurements, they are of theoretical interest. Loos et al. Loos, Lipparini, and Jacquemin (2023) also reported TBEs of their STGs as $-0.071$ eV (2AP), $-0.042$ eV (3AP), and $-0.029$ eV (4AP). In the calculations of TBEsLoos, Lipparini, and Jacquemin (2023), geometries were determined using the accurate method, CCSD(T), by enforcing the maximum symmetry point group for each composition: $D_{\rm 3h}$ for 1AP, 4AP, and 7AP; $C_{\rm 2v}$ for 2AP, 5AP, and 3AP. For this purpose, MP2 minimum energy structures verified through frequency analysis were selected as initial geometries. Other studies have also applied geometries determined with MP2de Silva (2019); Tučková et al. (2022) or using DFT methods such as B3LYPGhosh and Bhattacharyya (2022) and $\omega$B97X-DGarner, Blaskovits, and Corminboeuf (2024).

Elaborate analysis using wavefunction methods and large basis sets have scrutinized the validity of Hund’s rule violation by APsDreuw and Hoffmann (2023). Other works have shown the limitations of popular density functional theory (DFT) methods in predicting the negative STGs of APsGhosh and Bhattacharyya (2022); Kondo (2022); Tučková et al. (2022). The nature of the potential energy surface (PES) of APs and its impact on their STGs has received less attention. For 1AP, early calculationsLeupin and Wirz (1980) suggested a $C_{\rm 3h}$ structure with alternating bond lengths, contrasting the $D_{\rm 3h}$ structure typical of an aromatic system. However, calculations using DFTde Silva (2019) and MP2Ehrmaier et al. (2019) methods predicted a symmetric structure, which has since been widely adopted. Although more accurate calculations of 1AP have been performed using the CCSD(T) methodLoos, Lipparini, and Jacquemin (2023), vibrational analysis at this level remains challenging when using a triple-zeta quality basis set. Additionally, determining an initial geometry for CCSD(T)-level explorations to assess the relative stability of iso-energetic structures is non-trivial, as approximate methods exclusively prefer symmetric structures. Consequently, the true nature of the PES of 1AP remains elusive.

Given the complexity of the problem, an out-of-the-box approach is necessary to investigate the structural preferences of APs. We examine symmetry-lowering in APs from the highest point group symmetry, permitted by their composition and connectivities, to their maximum-order subgroups. To probe for symmetry lowering and the variation of STG on a lower-dimensional potential energy surface, we conducted a diagnostic procedure, which we denote as Jahn–Teller–Hund’s (JTH) diagnostics. We determine the PES of APs using the CCSD(T) method, which is more accurate than MP2, to probe for changes in qualitative aspects of the PES due to electron correlation effects. Further, we inspect the variation of the STG across the PES.

II Computational details

Geometry optimizations of APs (FIG. 1) were conducted by relaxing the internal coordinates defined through the Z-matrix representation. For JTH diagnostics, two interatomic distances ($r_{1}$ and $r_{2}$ in FIGs. 2 and 4) were selected to sample geometries for constructing two-dimensional PESs. Each point on the PES was obtained using the MP2 method and the cc-pVDZ basis set via constrained geometry optimization with fixed values of $r_{1}\in[1.30\,{\mathrm{\AA}},1.50\,{\mathrm{\AA}}]$ and $r_{2}\in[1.30\,{\mathrm{\AA}},r_{1}]$, in steps of 0.01 Å. For plotting, discrete points were interpolated using cubic functions. Single point calculations of the ground and excited state energies were calculated on the PESs using the CCSD(T) and ADC(2) methods, respectively, in combination with the cc-pVDZ basis set. For each system, the minimum energy structure along the $r_{1}=r_{2}$ line was selected for determining the geometries corresponding to the high-symmetry point groups: $D_{\rm 3h}$ (for 1AP, 4AP, and 7AP) and $C_{\rm 2v}$ (for 2AP, 3AP, and 5AP). For 1AP, 2AP, 3AP, and 4AP, CCSD(T) energies on the PESs identified more stable low-symmetry structures (with $r_{1}\neq r_{2}$), which were selected as initial structures for subsequent geometry refinement. For all molecules, equilibrium geometries were determined using the DFT methods, B3LYPStephens et al. (1994) and $\omega$B97X-DChai and Head-Gordon (2008), along with the wavefunction methods MP2, CCSD and CCSD(T) in combination with the cc-pVTZ basis set. Harmonic frequency analyses were performed with DFT and MP2 methods to confirm if the structures were minima or saddle points on their respective potential energy surfaces.

DFT-level geometry optimizations were conducted using Gaussian (version 16 C.01)Frisch et al. (2016). MP2, CCSD, and CCSD(T) geometry optimizations were performed using Molpro (version 2015.1)Werner et al. (2015). For determining STGs, S₁ and T₁ energies were calculated using the excited states method, ADC(2), employing the resolution-of-identity (RI) approximationVahtras, Almlöf, and Feyereisen (1993); Kendall and Früchtl (1997) as implemented in QChem (version 6.0.2)Krylov and Gill (2013). All calculations based on the correlated wavefunction methods MP2, CCSD, CCSD(T), and ADC(2) were performed with the frozen-core approximation. Complete basis set (CBS) extrapolation was performed using the two-point formula that uses CCSD(T) energies determined with the cc-pVTZ and cc-pVQZ basis sets. Hartree–Fock (HF) energy from the cc-pVQZ basis set is taken as the reference for the CBS estimate, while the correlation energy is extrapolated using the formula $E_{\rm corr}^{n}=E_{\rm corr}^{\rm CBS}+\alpha n^{-3}$, where $n$ is the cardinal number of the basis set. Normal mode and nucleus-independent chemical shifts (NICS)Schleyer et al. (1996) analyses were conducted with the $\omega$B97X-D3 method and the cc-pVTZ basis set using Orca (version 5.0.4)Neese (2012, 2018). For both analyses, minimum energy geometries calculated using $\omega$B97X-D/cc-pVTZ were utilized. The NICS value was calculated as the negative of the magnetic shielding of a ghost atom placed 1 Å above the centroid of each ring to determine isotropic and out-of-plane NICS, denoted NICS(1)_iso and NICS(1)_zz, respectively.

Figure 2: Contour plots of two-dimensional potential energy surfaces and corresponding S₁-T₁ gaps for 1AP, 5AP, and 7AP. MP2 (relaxed) indicates that each point on the surface was calculated through constrained geometry optimization at the MP2-level for fixed values of $r_{1}\in[1.30\,{\mathrm{\AA}},1.50\,{\mathrm{\AA}}]$ and $r_{2}\in[1.30\,{\mathrm{\AA}},r_{1}]$, in steps of 0.01 Å. Discrete points were interpolated as a surface using cubic functions. The point group symmetries were constrained to $C_{\rm 3h}$  for 1AP and 7AP, and $C_{\rm s}$  for 5AP; along $r_{1}=r_{2}$, the symmetries are $D_{\rm 3h}$ and $C_{\rm 2v}$ , respectively. CCSD(T)@MP2 indicates the ground state energy, while ADC(2)@MP2 indicates S₁-T₁ gap calculated on geometries optimized with MP2 at each point. The ground state energy (in kJ/mol) and the S₁-T₁ gap (in eV) are color-coded by rainbow and red-to-blue color bars, respectively. Energies above 20 kJ/mol are not shown. In all calculations, cc-pVDZ basis set was used.

Figure 3: Equilibrium geometry parameters of 1AP in $D_{\rm 3h}$ (saddle point) and $C_{\rm 3h}$ (minimum) configurations calculated with the CCSD(T)/cc-pVTZ method. Dashed lines show two-fold rotation axes, perpendicular to the principal (three-fold rotation) axis, present in $D_{\rm 3h}$ . MP2 describes the PES of 1AP as a single-well with the $D_{\rm 3h}$ structure to be the minimum, while CCSD and CCSD(T) suggest a double-well behavior with the lower symmetric $C_{\rm 3h}$ structure as the minimum. At the CCSD(T)-level, the barrier to automerization is vanishingly small, indicating the vibrational ground state’s wavefunction is strongly delocalized over the equivalent $C_{\rm 3h}$ minima.

Figure 4: Contour plots of two-dimensional potential energy surfaces and corresponding S₁-T₁ gaps for 2AP, 3AP, and 4AP. The point group symmetries were constrained to $C_{\rm 3h}$  for 4AP, $C_{\rm s}$  for 2AP and 3AP; along $r_{1}=r_{2}$, the symmetries are $D_{\rm 3h}$ and $C_{\rm 2v}$ , respectively. See the caption of FIG. 2 for further details.

III Results and Discussion

We begin our analysis with the two-dimensional PESs of 1AP, 5AP, and 7AP, obtained through constrained geometry optimizations by fixing two interatomic distances related by two-fold rotation, see FIG. 2. Using the Z-matrix representation, we constrained the molecules to be in the $C_{\rm 3h}$ (for 1AP and 7AP) and $C_{\rm s}$ (for 5AP) subspaces. MP2-level ground state electronic energy surface shows prominent single-well behavior for all three systems with 1AP, 5AP, and 7AP preferring $D_{\rm 3h}$ , $C_{\rm 2v}$ , and $D_{\rm 3h}$ structures with $r_{1}=r_{2}=1.42$ Å, 1.41 Å, and 1.40 Å, respectively. Subsequent MP2/cc-pVTZ-level geometry optimizations and harmonic vibrational frequency analysis described these symmetric structures as minima on the PESs.

CCSD(T) energies calculated in a single-point fashion identify two equivalent minima of 1AP corresponding to $C_{\rm 3h}$ symmetry at $\{r_{1}=1.44,\,r_{2}=1.41\}$ and $\{r_{1}=1.41,\,r_{2}=1.44\}$, see FIG. 2. The $D_{\rm 3h}$ structure predicted as a minimum by MP2 shows characteristics of a saddle point at the CCSD(T)-level. It is important to note that the transition state structure is more sensitive to electron correlation effects and MP2 has a tendency to overestimate the importance of double excitations for transition state structuresScuseria and Schaefer (1990). Another classic example is cyclobutyne, which MP2 predicts as a minimum but found to be a first-order saddle point by the more accurate CCSD(T) and CCSDT methodsSun and Schaefer III (2019). As the HOMO-LUMO gap increases upon substitution of N in the electron-rich sites of 1AP—with values of about $6/8/10$ eV for 1AP/5AP/7AP, respectively—electron correlation effects become less relevant. Hence, the PESs of 5AP and 7AP remain similar at both MP2 and CCSD(T) levels. PES plots calculated with the DFT methods, B3LYP and $\omega$B97X-D, are on display in Figure S1 in the Supplementary Information (SI). Both methods forecast single-well type PES for 1AP, 5AP, and 7AP.

In general, for probing the inverted nature of STGs of APs, correlated excited state methods such as ADC(2) or CC2 have provided good estimatesTučková et al. (2022); Loos, Lipparini, and Jacquemin (2023). STGs calculated on the PES with the ADC(2) method and cc-pVDZ basis set identifies regions where S₁ is lower in energy than T₁, see FIG. 2. For all three molecules, symmetric structures exhibit inverted STG. Individual plots of S₁ and T₁ relative to S₀ are on display in Figure S2, where one notices that the S₁ energy drops more rapidly when deviating from the diagonal line (corresponding to symmetric structures) compared to the T₁ energy. Overall, the JTH diagnostics on a low-dimensional PES hints at a symmetry-lowered structure for 1AP with an STG of diminished magnitude while suggesting symmetric structures for 5AP and 7AP with negative STG in agreement with experimental measurementsEhrmaier et al. (2019); Wilson et al. (2024).

Table 1 presents STGs of the selected APs shown in FIG. 1 calculated with ADC(2)/aug-cc-pVTZ using equilibrium geometries determined using the accurate CCSD(T)/cc-pVTZ method. Results based on geometries from other methods are available in the SI (see Tables S1, S2, and S3 for 1AP, 5AP, and 7AP, respectively; Table S4 for 2AP, 3AP, and 4AP). FIG. 3 shows the corresponding equilibrium geometry parameters for both the structures at the CCSD(T)/cc-pVTZ-level. For the $D_{\rm 3h}$ structure of 1AP, our ADC(2) value of the STG=$-$0.142 eV, deviates slightly from the previously obtainedLoos, Lipparini, and Jacquemin (2023) TBE of $-$0.131 eV. Using the basis set limit values of CCSD(T) energies, we find the $C_{\rm 3h}$ structure of 1AP to be lower in energy than the $D_{\rm 3h}$ structure by only 0.3 kJ/mol.

Had the symmetry lowering been predicted by MP2 or DFT methods but not at the more accurate CCSD(T)-level, then the distortion should be considered artifactualDavidson and Borden (1983). According to the accurate CCSD(T)-level description, the actual PES of 1AP should be thought of as a very shallow double well permitting strong carbon tunnelingZuev et al. (2003), which according to quantum mechanics has to be interpreted as two equivalent $C_{\rm 3h}$ structures undergoing rapid automerization. Furthermore, as pointed out earlier by Leupin et al. the nuclear magnetic resonance (NMR) spectrum of 1AP does not rule out the possibility of a rapid isodynamic equilibrium between two $C_{\rm 3h}$ minimaLeupin and Wirz (1980), which aligns with its very low barrier reported in this study.

For 1AP in the $C_{\rm 3h}$ structure, the ADC(2)/aug-cc-pVTZ value of STG is $-0.06$ eV. Post-ADC(2), high-level corrections can be estimated as $+0.009$ eV by comparing the STG of the $D_{\rm 3h}$ structure from ADC(2) ($-$0.142 eV) with the TBE ($-$0.131 eV). Furthermore, for the $D_{\rm 3h}$ structure of 1AP, the vertical STG underestimates the adiabatic transition energy (0-0), which accounts for geometric effects in S₁ and T₁ along with the zero-point vibrational energy, by +0.057 eV at the ADC(2)-levelTučková et al. (2022), see Table 1. Using these corrections, one can estimate the 0-0 STG of $C_{\rm 3h}$ structure to be $+0.006$ eV explaining the experimental value of +0.04 eVLeupin and Wirz (1980). On the other hand, 0-0 correction (+0.057 eV) when added to the TBE value of the $D_{\rm 3h}$ structure results in a too negative value, $-$0.074 eV. Given the vanishing barrier of 0.3 kJ/mol for automerization, one can expect the wavefunction of the vibrational ground state to span double well broadly. Hence, when the 0-0 STG of 1AP is weighted by the square of this wavefunction, one can expect STG to be vanishingly small, as in the $C_{\rm 3h}$ structure.

We performed separate calculations for high-symmetry and low-symmetry geometries of 5AP and 7AP, and found no evidence for symmetry lowering via in-plane distortions. For 5AP and 7AP, STGs calculated with ADC(2)/aug-cc-pVTZ based on geometries determined with CCSD(T)/cc-pVTZ show deviations of about 0.03 eV from TBEs. The structural stability of 5AP and 7AP, confirmed by the CCSD(T)-level analysis, indicates that both molecules violate Hund’s rule in their minimum energy geometry with high-symmetry. From previously reportedTučková et al. (2022) ADC(2)-level vertical and 0-0 values of STG, one can correct the TBE of the vertical STG of 5AP from Ref. Loos, Lipparini, and Jacquemin, 2023 to be $-$0.07 eV, agreeing with the experimental value of $-0.047$ eVWilson et al. (2024). As the precise value of the STG of 7AP is unknown, the TBE values corrected for 0-0 excitation amounts to $-0.083$ eV, which is lower than $-$0.011 eV from the fluorescence measurementsAizawa et al. (2022) of a substituted 7AP.

Table 1: For the six azaphenalenes shown in FIG. 1, vertical excitation energies of the S₁ and T₁ states with respect to the ground state, S₀, along with the singlet-triplet gap (S₁-T₁) are given in eV. The excited state energies were calculated using the ADC(2)/aug-cc-pVTZ method using geometries determined with the CCSD(T)/cc-pVTZ method. The barrier for automerization, $E^{\ddagger}$ (in kJ/mol), determined using two-point CBS extrapolation of CCSD(T) energies with cc-pVTZ and cc-pVQZ basis sets, is stated for the symmetric saddle point. Results from other studies are included for comparison; 0-0 indicates adiabatic transition energies accounting for the energy of the vibrational ground state.

System	$E^{\ddagger}$	S1	T1	S1-T1	Source
1AP ($D_{\rm 3h}$ )	0.3	0.979	1.121	$-$0.142	this work
1AP ($C_{\rm 3h}$ )		1.130	1.190	$-$0.060	this work
1AP		1.047	1.180	$-$0.133	ADC(2)Tučková et al. (2022)
1AP		0.992	1.068	$-$0.076	0-0, ADC(2)Tučková et al. (2022)
1AP ($D_{\rm 3h}$ )		0.979	1.110	$-$0.131	TBELoos, Lipparini, and Jacquemin (2023)
1AP		0.97	0.93	+0.04	exp.Leupin and Wirz (1980)a
5AP ($C_{\rm 2v}$ )		2.128	2.274	$-$0.147	this work
5AP		2.231	2.365	$-$0.134	ADC(2)Tučková et al. (2022)
5AP		1.971	2.056	$-$0.085	0-0, ADC(2)Tučková et al. (2022)
5AP ($C_{\rm 2v}$ )		2.177	2.296	$-$0.119	TBELoos, Lipparini, and Jacquemin (2023)
5AP		1.957	2.003	$-$0.047	0-0, exp.Wilson et al. (2024)
7AP ($D_{\rm 3h}$ )		2.636	2.889	$-$0.253	this work
7AP		2.756	2.998	$-$0.242	ADC(2)Tučková et al. (2022)
7AP		2.512	2.618	$-$0.106	0-0, ADC(2)Tučková et al. (2022)
7AP ($D_{\rm 3h}$ )		2.717	2.936	$-$0.219	TBELoos, Lipparini, and Jacquemin (2023)
7AP				$\quad<$0	exp.Ehrmaier et al. (2019)
2AP ($C_{\rm 2v}$ )	2.1	0.838	0.941	$-$0.102	this work
2AP ($C_{\rm s}$ )		1.246	1.121	$+$0.125	this work
2AP ($C_{\rm 2v}$ )		0.833	0.904	$-$0.071	TBELoos, Lipparini, and Jacquemin (2023)
3AP ($C_{\rm 2v}$ )	5.3	0.689	0.768	$-$0.079	this work
3AP ($C_{\rm s}$ )		1.335	1.061	$+$0.274	this work
3AP ($C_{\rm 2v}$ )		0.693	0.735	$-$0.042	TBELoos, Lipparini, and Jacquemin (2023)
4AP ($D_{\rm 3h}$ )	11.1	0.550	0.623	$-$0.073	this work
4AP ($C_{\rm 3h}$ )		1.409	1.017	$+$0.392	this work
4AP ($D_{\rm 3h}$ )		0.554	0.583	$-$0.029	TBELoos, Lipparini, and Jacquemin (2023)

• ^a Using S₁=0.78 $\mu$m^-1 and T₁=0.75 $\mu$m^-1 from Ref. Leupin and Wirz, 1980 multiplied by 1.2398 eV/$\mu$m^-1.

FIG 4 presents the JTH diagnostic plots for 2AP, 3AP, and 4AP. MP2 describes these systems as prominent single-well structures of maximum symmetry ($C_{\rm 2v}$ for 2AP and 3AP; $D_{\rm 3h}$ for 4AP). MP2/cc-pVTZ geometry refinements and frequency analysis indicated all three symmetric structures to be minima. CCSD(T) energies show a trend contrasting with MP2-level description for all three molecules; in all cases, the high-symmetry structure turns out to be the saddle point connecting low-symmetry structures belonging to the largest subgroup ($C_{\rm s}$ for 2AP and 3AP; $C_{\rm 3h}$ for 4AP ). As seen in FIG 4, the barrier increases with the number of N atoms on charge-destabilizing positions (i.e., LUMO sites of 1AP, see FIG.1). STG plots of all three molecules indicate a profile similar to the one noted for 1AP in FIG. 2. The STG is negative along the $r_{1}=r_{2}$ line containing the high-symmetry structures, while for the distorted structures, the plot suggests positive STGs. CCSD(T)/cc-pVTZ geometries were determined using the coordinates of the distorted structures inferred from FIG 4 as the starting point. Subsequent ADC(2)/aug-cc-pVTZ excited state calculations for 2AP, 3AP, and 4AP resulted in STGs of +0.125 eV, $+$0.274 eV, and +0.392 eV, respectively (see Table 1).

Figure S1 shows the PESs of 2AP, 3AP, and 4AP generated using the DFT methods B3LYP and $\omega$B97X-D. While B3LYP predicts 2AP and 3AP to be single-wells, for 4AP, the method suggests a slight distortion with $r_{1}$ and $r_{2}$ differing by 0.03 Å. Frequency analysis confirmed the $D_{\rm 3h}$ and $C_{\rm 2v}$ structures of 2AP and 3AP to be minima at the B3LYP/cc-pVTZ-level, while the symmetric structure of 4AP showed one imaginary frequency (494$i$ cm^-1, $A_{2}^{\prime}$). The long-range corrected DFT method, $\omega$B97X-D, has been shown to outperform B3LYP in determining minimum energy structures of unusual covalent bond connectivitiesSenthil, Chakraborty, and Ramakrishnan (2021). For 2AP, $\omega$B97X-D shows a weak distortion in agreement with CCSD(T). At the same time, the method indicates the symmetric geometries of 2AP, 3AP, and 4AP to be saddle points, suggesting the possibility of locating the low-symmetry minima with less expensive DFT modeling (see Figure S1). We found the NICS(1)_iso values of 5AP and 7AP to have very small magnitudes, indicating these two systems to be non-aromatic with localized charge density. The other four systems in their symmetric structure exhibit large positive NICS(1)_iso values 1 Å above the centroid of their rings indicating anti-aromaticityKaradakov (2008), in the order: 1AP $<$ 2AP $<$ 3AP $<$ 4AP (see Figure S4).

We rationalize the symmetry lowering in APs as a vibronically driven process on the basis of the Jahn–Teller (JT) theoremJahn and Teller (1937); Bader (1962); Pearson (1975) . Specifically, we verify if the symmetry of the vibrational normal mode corresponding to the distortion, as seen in DFT results, agrees with the JT theorem. The electronic ground state energy, $E_{0}(q)$, is expanded as a truncated Taylor series in the distortion coordinate, $q$:

$$E_{0}(q) = E_{0}+\langle 0|\hat{H}^{\prime}|0\rangle q+\frac{1}{2}kq^{2},$$ (1)

where $E_{0}$ corresponds to the ground state energy in the undistorted equilibrium geometry (at $q=0$), and $\hat{H}^{\prime}$ is the first derivative of the electronic Hamiltonian with respect to $q$. At $q=0$, the coefficient in the first-order JT contribution, $\langle 0|\hat{H}^{\prime}|0\rangle$, vanishes for APs; however, this term can be non-zero for systems with degenerate electronic states. For pseudo-degenerate (i.e., near degenerate) electronic states, the second-order term in $q^{2}$ can be non-zero, where $k$ is the total force constant that includes contributions from vibronic couplingNakajima, Toyota, and Kataoka (1982)

$$k = \langle 0|\hat{H}^{\prime\prime}|0\rangle-2\sum_{n\neq 0}\frac{|\langle 0|\hat{H}^{\prime}|n\rangle|^{2}}{E_{n}-E_{0}}.$$ (2)

Here $\hat{H}^{\prime\prime}$ is second derivative of the electronic Hamiltonian with respect to $q$, and $E_{n}$ is the energy of the $n$-th electronic state at $q=0$.

The second-order term gives rise to the pseudo-JT contribution capable of stabilizing the system upon distortion from $q=0$, provided there is a low-lying state, $n$, of correct symmetry such that $\langle 0|\hat{H}^{\prime}|n\rangle\neq 0$. For this, the product of the corresponding irreducible representations, $\Gamma_{0}\times\Gamma_{q}\times\Gamma_{n}$, should contain the totally symmetric representation of the point group corresponding to the symmetric structure. If this group is Abelian, then the condition is $\Gamma_{q}=\Gamma_{0}\times\Gamma_{n}$. For 1AP and 4AP, the totally symmetric irreducible representation is $A_{1}^{\prime}$ of $D_{\rm 3h}$ , and for 2AP and 3AP, it is $A_{1}$ of $C_{\rm 2v}$ . We write the irreducible representations in small case alphabets when denoting the MOs and use large cases for the states and vibrational modes.

Assuming that there is only one low-lying excited state, $|1\rangle$, coupling with the ground state, we arrive at

$$k = \langle 0|\hat{H}^{\prime\prime}|0\rangle-2\frac{|\langle 0|\hat{H}^{\prime}|1\rangle|^{2}}{E_{1}-E_{0}}=k_{0}-k_{0,1}.$$ (3)

Here, $k_{0}\geq 0$ is the force constant of the symmetric structure in the absence of vibronic coupling, while the condition for distortion is $k<0$. The second term of Eq. 3, $k_{0,1}$, cannot be negative, as it is a ratio of two positive numbers. Hence, symmetry-lowering is prominent when $k_{0}<k_{0,1}$, which requires a significantly large $|\langle 0|\hat{H}^{\prime}|1\rangle|^{2}$, small $E_{1}-E_{0}$ or both. This explains why symmetry lowering is not observed in 5AP and 7AP with fairly large S₀$\rightarrow$S₁ excitation energies of about 2.1, and 2.6 eV (see TABLE 1). On the other hand, the excitation energy of S₁ is 1.0 eV in 1AP, which decreases to 0.8, 0.7, and 0.6 eV in 2AP, 3AP, and 4AP, respectively. Overall, our analysis suggests that inverted-STG candidates with small values of ${\rm S}_{0}\rightarrow{\rm S}_{1}$ may exhibit strong structural distortions warranting more detailed analysis such as the JTH diagnostics presented in this study.

The symmetry of the ground electronic state of the APs, in their high-symmetry geometry, correspond to the totally symmetric irreducible representations: $A_{1}$ of $C_{\rm 2v}$ and $A_{1}^{\prime}$ of $D_{\rm 3h}$ . If the irreducible representation of $|1\rangle$ is non-degenerate, the product of the irreducible representations of $\hat{H}^{\prime}$ and $|1\rangle$ should be totally symmetric, which is possible only when $q$ and $|1\rangle$ belong to the same irreducible representation. The character of the HOMO and LUMO in all APs are similar, see Figure S3. For molecules in $D_{\rm 3h}$ point group, the symmetries are $a_{1}^{\prime\prime}$ and $a_{2}^{\prime\prime}$, respectively, while for those in $C_{\rm 2v}$ , the symmetries are $a_{2}$ and $b_{1}$, respectively. The symmetry of the excited state, $|1\rangle$, when considering a Slater-determinant form, is the product, $\Gamma_{\rm HOMO}\times\Gamma_{\rm LUMO}$. Hence, the symmetry of the excited state under concern is $A_{2}^{\prime}$ for 1AP and 4AP (in $D_{\rm 3h}$ ) and $B_{2}$ for 2AP and 3AP (in $C_{\rm 2v}$ ).

Figure 5: Pseudo-Jahn–Teller distortion in 4AP according to $\omega$B97X-D/cc-pVTZ modeling: a) Potential energy profile of 4AP along the Jahn–Teller active dimensionless normal coordinate, $q$, of $A_{2}^{\prime}$ symmetry. The force constant matrix was computed at the $D_{\rm 3h}$ geometry ($q=0$), a first-order saddle point with an imaginary frequency (1113$i$ cm^-1); the inset shows the normal mode. The DFT ground electronic state energy is denoted S₀ (pink) and the symmetry lowers to $C_{\rm 3h}$ when $q\neq 0$. The excited state energies, S₁ (blue) and T₁ (red), are determined using the ADC(2)/cc-pVTZ method. The singlet-triplet gap, S₁-T₁, (green) is negative within the shaded region $-0.3\leq q\leq 0.3$. b) Frontier MOs of 4AP in $D_{\rm 3h}$ and $C_{\rm 3h}$ structures showing an increase in the HOMO-LUMO gap due to better symmetry-allowed atomic orbital interactions in the $C_{\rm 3h}$ structure.

In our $\omega$B97X-D/cc-pVTZ calculations, the wavenumbers of the distortion modes of 2AP, 3AP, and 4AP are 340$i$ cm^-1 ($B_{2}$), 773$i$ cm^-1 ($B_{2}$), and 1113$i$ cm^-1 ($A_{2}^{\prime}$), respectively; see Table S4 in the SI. The corresponding $k$ (in Eq. 3) are $-$0.64, $-$3.82, and $-$8.82 mdyne/Å, respectively. This trend shows how the distortion becomes stronger with increasing N in the electron-deficient sites of 1AP. The Jahn–Teller active normal mode corresponds to in-plane stretching, with symmetries $A_{2}^{\prime}$ (for $D_{\rm 3h}$ $\rightarrow$ $C_{\rm 3h}$ transition) and $B_{2}$ (for $C_{\rm 2v}$ $\rightarrow$ $C_{\rm s}$ transition). Upon symmetry lowering, the irreducible representation is $A^{\prime}$, which is the totally symmetric representation in $C_{\rm 3h}$ and $C_{\rm s}$ . This observation follows the epikernel principle: preferred distortions of JT unstable molecules are directed towards the largest subgroup such that the JT active mode in the subgroup is totally symmetric Ceulemans, Beyens, and Vanquickenborne (1984).

FIG. 5 shows the ground and excited state energies of 4AP along the JT active mode of the $D_{\rm 3h}$ configuration. The STG remains negative for slight distortions from the equilibrium geometry. At $q\thickapprox$1.8 (corresponding to the $C_{\rm 3h}$ minimum), the STG increases to +0.344 eV (see Table S4). The main driving force for the pseudo-JT effect is the improved orbital interactions permitted by symmetry lowering. Analysis of the HOMO and LUMO plots of 4AP in $D_{\rm 3h}$ and $C_{\rm 3h}$ configurations reveals that the $a_{1}^{\prime\prime}$ and $a_{2}^{\prime\prime}$ MOs (HOMO and LUMO) transform as $a^{\prime\prime}$ upon symmetry lowering, see FIG. 5. A mild splitting of the HOMO-LUMO gap is noted upon distortion, imparting more stability to the ground state. Subsequently, symmetry lowering enhances the exchange interaction between HOMO and LUMO in the $C_{\rm 3h}$ structure, resulting in a more stable T₁ state compared to S₁.

IV Conclusions

In conclusion, we conducted diagnostic calculations to determine if the symmetric structures of APs are true minima on the PES and to investigate the dependence of their STG on structural distortions. Out of six APs studied here, the $D_{\rm 3h}$ structure of the prototypical molecule, cyclazine (1AP), is predicted to be a saddle point by CCSD(T) with the correct minimum identified to have the $C_{\rm 3h}$ symmetry. The mildly anti-aromatic nature of the ring-current in the symmetric form suggests the possibility of symmetry-lowering, as in the case of cyclobutadiene in the $D_{\rm 4h}$ structureKaradakov (2008); while the latter exhibits first-order JT effect, in the closed-shell system 1AP, the effect stems from second-order interactions. When including high-level and adiabatic corrections, our estimation of the STG of 1AP in the $C_{\rm 3h}$ structure explains the experimental value of 0.04 eVLeupin and Wirz (1980) better than the more negative STG of the $D_{\rm 3h}$ structure. Two APs (5AP and 7AP) with N atoms at electron-rich sites do not show symmetry lowering at the CCSD(T)-level and exhibit negative STG in ADC(2) calculations in line with the experimental resultsEhrmaier et al. (2019); Wilson et al. (2024) and TBEsLoos, Lipparini, and Jacquemin (2023). Three APs (2AP, 3AP, and 4AP) with N at the electron-deficient sites show stronger tendencies for symmetry lowering, even when modeled with a long-range corrected DFT method. At the correct minima determined with CCSD(T), they show positive STGs. Popular methods, MP2 and B3LYP, prefer symmetric structures leading to negative STGs in subsequent excited state modeling, thereby suggesting a violation of Hund’s rule. More rigorous studies should be conducted on the sensitivity of STGs to symmetry constraints of inverted-STG candidates not explored in this study. It is essential to ensure that novel properties capable of a paradigm shift in the theoretical understanding of molecules withstand the artifacts of quantum chemistry approximations.

V Supplementary Information

i) Tables S1–S4 compares STGs from our calculations to results from past works; ii) Figure S1 presents two-dimensional surface plots of S₁ and T₁; iii) Figure S2 presents two-dimensional surface plots of ground state energies calculated with DFT methods; iv) Figure S3 presents plots of HOMO and LUMO; v) Figure S4 presents NICS(1)_iso and NICS(1)_zz values; vi) Frozen-core CCSD(T)/cc-pVTZ equilibrium coordinates.

VI Data Availability

The data that support the findings of this study are within the article and its supplementary material.

VII Acknowledgments

We acknowledge the support of the Department of Atomic Energy, Government of India, under Project Identification No. RTI 4007. All calculations have been performed using the Helios computer cluster, which is an integral part of the MolDis Big Data facility, TIFR Hyderabad (http://moldis.tifrh.res.in).

VIII Author Declarations

VIII.1 Author contributions

AM: Conceptualization (main); Analysis (equal); Data collection (equal); Writing (main). KJ: Analysis (equal); Data collection (equal). SD: Analysis (equal); Data collection (equal). RR: Conceptualization (main); Analysis (equal); Data collection (equal); Funding acquisition; Project administration and supervision; Resources; Writing (main).

VIII.2 Conflicts of Interest

The authors have no conflicts of interest to disclose.

IX References