Performance of Electronic Structure Methods for the Description of Hückel–Möbius Interconversions in Extended π-Systems

Expanded porphyrins provide a versatile route to molecular switching devices due to their ability to shift between several π-conjugation topologies encoding distinct properties. DFT remains the workhorse for modeling such extended macrocycles, when taking into account their size and huge conformational flexibility. Nevertheless, the stability of Hückel and Möbius conformers depends on a complex interplay of different factors, such as hydrogen bonding, π···π stacking, steric effects, ring strain, and electron delocalization. As a consequence, the selection of an exchange–correlation functional for describing the energy profile of topological switches is very difficult. For these reasons, we have examined the performance of a variety of wave function methods and density functionals for describing the thermochemistry and kinetics of topology interconversions across a wide range of macrocycles. Especially for hexa- and heptaphyrins, the Möbius structures have a stronger degree of static correlation than the Hückel and twisted-Hückel structures, and as a result the relative energies of singly twisted structures are a challenging test for electronic structure methods. Comparison of limited orbital space full CI calculations with CCSD(T) calculations within the same active spaces shows that post-CCSD(T) correlation contributions to relative energies are very minor. At the same time, relative energies are weakly sensitive to further basis set expansion, as proven by the minor energy differences between the extrapolated MP2/CBS energies estimated from cc-pV{T,Q}Z, diffuse-augmented heavy-aug-cc-pV{T,Q}Z and explicitly correlated MP2-F12/cc-pVDZ-F12 calculations. Hence, our CCSD(T) reference values are reasonably well-converged in both 1-particle and n-particle spaces. While conventional MP2 and MP3 yield very poor results, SCS-MP2 and particularly SOS-MP2 and SCS-MP3 agree to better than 1 kcal mol–1 with the CCSD(T) relative energies. Regarding DFT methods, the range-separated double hybrids, such as ωB97M(2) and B2GP-PLYP, outperform other functionals with RMSDs of 0.6 and 0.8 kcal mol–1, respectively. While the original DSD-PBEP86 double hybrid performs fairly poorly for these extended π-systems, the errors drop down to 1.9 kcal mol–1 for the revised revDOD-PBEP86-NL, which eliminates the same-spin correlation energy. Minnesota meta-GGA functionals with high fractions of exact exchange (M06-2X and M08-HX) also perform reasonably well, outperforming more robust and significantly less empirically parametrized functionals like SCAN0-D3.


■ INTRODUCTION
Expanded porphyrins have attracted considerable attention over the past few decades in view of their large conformational flexibility, facile redox interconversions, novel metal coordination behaviors, and versatile electronic states. 1 The rich chemistry of these extended macrocycles has led to diverse applications including near-infrared dyes, 2 nonlinear optical materials, 3 magnetic resonance imaging contrast agents, 4 and molecular switches. 5 In contrast to the regular porphyrin, expanded porphyrins are flexible enough to switch between different π-conjugation topologies (Huckel, Mobius, and twisted-Huckel), each with different properties (Figure 1). 6−8 Such a change of topology involves a Huckel−Mobius aromaticity switch in a single molecule and it can be induced by solvent, pH, and metalation, among others. 9,10 These Huckel−Mobius aromaticity switches combine both mechanical and π-electron switching, providing a new route to molecular optoelectronic devices. 11 Indeed, we have recently demonstrated the applicability of these unique aromaticity switches as conductance switching and bithermoelectric devices 12,13 and as efficient nonlinear optical switches. 14 Through extensive density functional theory calculations, 15,16 we have demonstrated that the molecular topology is highly influenced by the number of π-electrons and the size of the macrocycle. While [4n + 2] π-electron expanded porphyrins adopt Huckel conformations that are almost planar and highly aromatic, antiaromatic Huckel and aromatic Mobius conformers coexist in dynamic equilibrium for [4n] π-electrons expanded porphyrins. 17 The larger heptaphyrins and octaphyrins prefer a twisted-Huckel topology (also known as figure-eight conformation) in the neutral state, but the Mobius topology became the most stable in protonated species. 18,19 The stability of these conformers depends on the complex interplay of different factors, such as hydrogen bonding, π···π stacking, steric effects, ring strain, and aromaticity. 19 As a consequence, the selection of an exchange−correlation functional for describing the energy profile of these topological switches remains challenging. An illustrative example is the potential energy curve for the Huckel−Mobius interconversion of neutral [32]heptaphyrin computed with different density functionals ( Figure 2).
In previous studies, we investigated the performance of several exchange−correlation functionals in reproducing the molecular structure of different meso-substituted expanded porphyrins by comparison with the X-ray diffraction data. Our results indicate that the degree of bond-length alternation and π-electron delocalization along the conjugation pathway is highly dependent on the amount of exact Hartree−Fock (HF) exchange. 15,16,18,19 Functionals with a small amount of HF exchange predict more delocalized and bond-length equalized structures than functionals including a larger percentage of exact exchange. 20 Nevertheless, the bond lengths for Mobius and Huckel structures obtained with B3LYP and M06 are in better agreement with the crystallographic data than the bond lengths computed with M06-2X with double the amount of nonlocal exchange. 16,18,19 From the geometrical parameters, it seems that π−π stacking interactions in the twisted-Huckel topology are overemphasized in functionals like M06-2X,  which generally yields reduced inter-ring distances and interplanar angles for the stacked rings compared to the experimental geometries ( Figure S1). These geometrical changes influence the relative energies of twisted-Huckel and Mobius conformations of large expanded porphyrins, resulting in likely overestimation of the stability of twisted-Huckel topologies. As a consequence, the M06-2X functional fails in predicting the Mobius topology as the most stable configuration upon protonation of [32]heptaphyrin and [36]octaphyrins, 18,19 in disagreement with the spectroscopic data and solid-state structures. 21,22 The critical role of the density functional on the energetic, geometric and magnetic results of expanded porphyrins was recently established by Torrent-Sucarrat et al. 23 In this benchmark study, the performance of 11 DFT functionals was assessed with respect to the local pair natural orbital coupled cluster DLPNO-CCSD(T) method 24−26 for topological switches based on [26]hexaphyrin and [32]heptaphyrin. They conclude that CAM-B3LYP, M05-2X, and M06-2X functionals provide a more consistent description of these topological switches. Nevertheless, it is important to underline that the reference values were obtained with the DLPNO-CCSD(T 0 ) approach using default settings, in which off-diagonal Fock matrix elements are neglected in the (T) contribution. 27 In line with previous findings, 28,29 our recent efforts to assess the performance of localized coupled cluster methods for Huckel−Mobius interconversions suggest larger errors for Mobius-type structures for the DLPNO-CCSD(T 0 ) approach, even with TightPNO cutoffs. 30 For all these reasons, we decide to examine the performance of a large variety of density functionals and wave function methods for describing the thermochemistry and kinetics of topological switches involving the interconversion between Huckel, Mobius and twisted-Huckel structures. A number of exchange−correlation functionals, ranging from generalized gradient approximations to double hybrids, has been tested and their performance to describe such unique interconversions has been carefully assessed with respect to canonical CCSD(T)/CBS calculations. To this aim, we have considered a variety of topology interconversions in N-fused [24]penta-, [28]hexa-, and [32]heptaphyrins ( Figure 3). On the basis of our previous conformational analyses, 15,16,18 the more stable conformations for the different expanded porphyrins were selected. In the case of [28]hexaphyrin and [32]heptaphyrin, these conformers correspond to Huckel (H), Mobius (M), and twisted-Huckel topologies (F), also known as figure-eight structures, whereas only Huckel and Mobius structures correspond to energy minima in the smaller N-fused [24]pentaphyrin. Besides the performance on the description of the conformational energies, we considered the key isomerization transition states for these topological interconversions.

■ COMPUTATIONAL DETAILS
Calculations. The geometries of all stationary points involved in the topology interconversions of N-fused [24]pentaphyrin, [28]hexaphyrin, and [32]heptaphyrin were fully optimized and characterized by harmonic vibrational frequency calculations at the B3LYP 31 /6-311G(d,p) level of theory. The performance of the B3LYP hybrid functional on the geometries of expanded porphyrins was assessed in our previous benchmarks toward crystallographic data. 16,18,19 The nature of the stationary points was ascertained from the appropriate number of negative eigenvalues of the Hessian matrix from the frequency calculations. Minima are characterized by positive eigenvalues, whereas transition states exhibit one negative eigenvalue corresponding to the rotation of one internal dihedral angle. Starting from those optimized geometries, single-point energy calculations were performed using various DFT functionals and wave function methods. Most density functional calculations were carried out using either the Gaussian 09/Gaussian16 software package 32 or ORCA versions 4.0 through 4.2, 33 whereas Molpro 2018 34 was used to perform several Møller−Plesset (MP2) and coupled cluster calculations. MP2.X and MP2.5 calculations were conducted with the Turbomole program system. 35 As byproducts of the MP2 and MP2-F12 calculations, we also obtained spin-component-scaled varieties such as SCS-MP2-F12 36 and S2-MP2. 37 The spin-component scaled MP2 theory treats the correlation effects of electron pairs with opposite spin and same spin differently by means of different scaling parameters. 38,36 Next to SCS-MP2, MP2.5 and MP2.X methods were also considered. MP2.5 is constructed as the average between MP2 and MP3 energies and produces accurate energies at a computational cost significantly lower than that of CCSD(T). 39 MP2.X is a generalized MP2.5 method of the form 40 where C is a basis set-specific empirical scaling parameter for the third-order correction term (C = 0.72 for cc-pVDZ).
MP2.X is designed to reduce the errors exhibited by MP2.5 when using smaller basis sets. Spin-component scaled coupled-cluster singles and doubles, such as SCS-CCSD 41 and SCS(MI)-CCSD, 42 were also considered, as they are a zero-cost byproduct of the CCSD results. SCS(MI)-CCSD is a reparametrized version of the SCS-CCSD method, in which the same-and opposite-spin scaling parameters were optimized to minimize the error toward CCSD(T)/CBS interaction energies for various noncovalent interaction dimers. Furthermore, the performance of the distinguishable cluster approach DCSD 43 and its spincomponent scaled version SCS-DCSD 44 was assessed.
In our benchmark, the following set of DFT functionals was considered, classified according to Jacob's ladder: 45 • on the second rung corresponding to generalized gradient approximation functionals (GGA): BLYP, 46 PBE, 47 and revPBE 48 • on the third rung, we consider the following metageneralized gradient approximation functionals (mGGA): TM, 49 78 ), and a random-phaseapproximation dRPA-based double hybrid (dRPA75 79 ) We have also considered Grimme's DFT-D3 atom-pairwise dispersion corrections 80 with the Becke−Johnson (BJ) damping function 81 for most of the DFT energies. For the M06-2X functional, the zero-damping version D3(0) 82 was used since the parameter optimization for the M06-2X-D3(BJ) diverge. For the revDSD-PBEP86 and revDOD-PBEP86 functionals, the newly developed atomic-charge dependent London dispersion correction DFT-D4 83,84 dispersion correction was also considered, aside from the Vydrov−Van Voorhis "nonlocal" (NL) dispersion functional, 85 in which an a posteriori correction is obtained from the electron density.
For orbital-based ab initio calculations, we mostly employed correlation-consistent polarized basis sets. 86,87 For the canonical CCSD(T) calculations, the cc-pVDZ with no p-type polarization functions on hydrogen atoms was employed due to the systems' sizes. For the explicitly correlated methods MP2-F12, the cc-pVDZ-F12 basis set 88 together with the associated complementary auxiliary basis sets (CABS) 89 was used. The latter accelerate the basis set convergence of explicitly correlated methods. 90 For most DFT calculations, the Weigend−Ahlrichs basis set def2-TZVP 91 was employed using the corresponding auxiliary basis sets for simultaneously fitting Coulomb and exchange. 92 For several functionals, we compare the def2-TZVP basis set with the 6-311G(d,p) basis set, which is commonly used in the modeling of expanded porphyrins. 15,19,23 The small differences in the RMSDs ( Figure  S3) indicate that basis set sensitivity for the DFT energies of Huckel and Mobius expanded porphyrins is relatively weak. In addition, to assess the influence of the expansion of the basis set on the DFT energies, we explore the basis set convergence of several functionals using the def2-QZVP basis set (Table  S2), which approaches the CBS limit for the rung 1−4 functionals. 93,94 It is encouraging that the differences between the def2-TZVP and def2-QZVP are remarkably small, with RMSD differences below 0.3 kcal mol −1 . These results suggest that a triple-ζ basis set is sufficient for more efficient yet acceptably accurate DFT calculations for the considered systems.
Choice of the Reference Method. As a reference, we employed fully canonical CCSD(T) calculations, widely regarded as the gold-standard in quantum chemistry. 95 More specifically, relative energies evaluated at the CCSD(T)/cc-pVDZ and CCSD(T)/CBS levels serve as a reference to benchmark the relative energies obtained with wave function and DFT methods, respectively, for each family of macrocycles shown in Figure 3. The accuracy of the different computational methods was assessed through the mean unsigned error (MUE) and the root-mean-square deviation (RMSD) relative to the fully canonical CCSD(T) energies: It is noteworthy that the computational resources for performing canonical CCSD(T)/cc-pVDZ on these extended systems were enormous in terms of CPU and disk space, precluding the use of a larger basis set. To assess the effect of  88 and associated auxiliary basis sets 89 for all species. We note in passing that for the largest systems, even cc-pVDZ-F12 already required about 10 TB of scratch space per run. As the largest nodes at Weizmann only had about 5 TB of SSD scratch space available, we achieved this feat by cross-mounting scratch file systems from other nodes via NFS over InfiniBand. Typically, the explicitly correlated calculations with appropriate basis set gain about two angular momenta in basis set convergence. 97 Hence, the MP2-F12/cc-pVDZ-F12 energies ought to be comparable or superior to MP2/cc-pVQZ in terms of convergence. This statement was verified by carrying out additional RI-MP2/cc-pVTZ and cc-pVQZ calculations and extrapolating to the complete basis set limit using the Helgaker formula. 98 Interestingly, the resulting MP2/ cc-pV{T,Q}Z relative energies deviate from their MP2-F12/cc-pVDZ-F12 counterparts by just 0.1 kcal mol −1 RMSD. The basis set extension effect itself, from MP2/cc-pVDZ, is just 0.9   The CCSD(T)/CBS energies were obtained from extrapolation of the MP2 energies at the complete basis set limit from the (heavy-aug-)cc-pVTZ and (heavy-aug-)cc-pVQZ basis sets. Alternatively, the CCSD(T)/cc-pVDZ-F12 were estimated from explicitly correlated MP2-F12/cc-pVDZ-F12 calculations.
The Journal of Physical Chemistry A pubs.acs.org/JPCA Article kcal mol −1 RMS in both cases. As shown in Figure 4, the MP2 relative energies computed with the cc-pVDZ and the cc-pVDZ-F12 basis sets are almost perfectly linearly related (coefficient of determination R 2 = 0.998) and the MUE is 0.74 kcal mol −1 for the complete set of structures. We might hence safely assume that the coupling term C in the equation below is negligible: And thus, we can apply the familiar "high-level correction" (HLC) approximation: In order to further assess the influence of the diffuse functions, we performed additional MP2 calculations with the heavy-augcc-pV{T,Q}Z, in which only the heavy (non-hydrogen) atoms are augmented with diffuse functions. Similarly, the extrapolated MP2/CBS energies estimated from cc-pV{T,Q}Z and diffuse-augmented heavy-aug-cc-pV{T,Q}Z are rather similar. A perfect linear correlation is obtained for the two sets of MP2/CBS energies (R 2 = 1.000) and the MUE is only 0.14 kcal mol −1 for the complete set of structures ( Figure S2). Overall, these results reinforce the hypothesis that the relative energies of topology interconversions are comparatively insensitive to basis set expansion and that our estimated CCSD(T)/CBS should be close to the basis set limit. These results are in line with those previously obtained which showed that the expansion of the Pople basis set from 6-31G to 6-311G(d,p) hardly influences the relative stabilities and energy  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article barriers for topological switches based on hexaphyrins. 99 Our best estimates for the relative energies of the expanded porphyrin database are collected in Table 1.
Although the extrapolation to the complete basis set limit does not affect significantly the energetic description of Huckel−Mobius interconversions, energy deviations larger than 1 kcal mol −1 are found for several conformations. Accordingly, we employed the CCSD(T)/CBS energies as a reference to assess the performance of the DFT methods, whereas the CCSD(T)/cc-pVDZ energies were used to benchmark the lower-level electron correlation methods.

■ RESULTS AND DISCUSSION
Database of Huckel and Mobius Expanded Porphyrins. In the present work, we introduce a representative database of topological switches based on expanded porphyrins with varying ring size. The database covers a broad spectrum of structures with Huckel (H), Mobius (M), and twisted-Huckel (F) topologies for representative [4n] π-electron expanded porphyrins ( Figure 3). Besides the more stable configurations, we also considered the key isomerization transition states for the topological interconversions. Here, it is important to note that the topology switching is achieved via internal rotations, 11 without "dissecting" the ring. This interesting feature is displayed in Figure 2, which shows that the interconversion between the twisted-Huckel and the Mobius conformers is achieved only by variation of one torsional angle (φ 1 ). As a result of the rotation of φ 1 , a ct-aligned pyrrole ring is transformed into a tt-aligned subunit, leading to an odd number of trans bonds in the Mobius topology. 17 The macrocycle is preserved during the switching process, but the π system is temporarily broken when one torsion angle becomes close to 90°. The latter can be better visualized in Figure 5. For these reasons, these Huckel−Mobius inter-conversions can be regarded as challenging isomerization processes for future benchmark studies.
In the case of N-fused [24]pentaphyrin, the low-energy pathway for the rotation of an imine-type pyrrole ring corresponds to a two-step mechanism, as shown in Figure  S4. Accordingly, two transition structures are considered for the Huckel−Mobius interconversions in the pentapyrrolic system. 16 In the case of [28]hexaphyrin, two different pathways for the interconversion between the Huckel planar (28H) and the Mobius structure (28M) were considered, on the basis of the exhaustive study of the reaction mechanism reported by Torrent-Sucarrat et al. 99 The two pathways differ in the rotating carbon−carbon bond leading to the intermediate Mobius structures (28M 1A and 28M 2A ), as indicated in Figure 6. Mechanism A involves the rotation of the φ 1 angle (in blue), whereas mechanism B involves the variation of the dihedral angle φ 2 (in red). The second step in both mechanisms corresponds to the proton transfer between both nitrogen atoms, leading to the most stable tautomer of the Mobius topology (28M). The rate-determining step can be either the bond rotation or the proton transfer, depending on the rotating bond and the meso-substituents. 99,100 For this particular switch, several density functionals were benchmarked against CCSD(T)/6-31G data, concluding that both CAM-B3LYP and M05-2X provide the most reliable results. 99 Besides the Huckel−Mobius interconversion, we also considered the switching from 28M to 28F with the associated transition state (28TS 3 ), as shown in Figure S5.
For the [32]heptaphyrin, we have considered the interconversions between twisted-Huckel (32F), Mobius (32M), and Huckel (32H) topologies with the associated transition states (32TS 1 and 32TS 2 ). The schematic energetic profile for the three-level topology switching in [32]heptaphyrin is displayed in Figure S6. In total, 21 structures Multireference Character of Expanded Porphyrins. In order to assess the reliability of single reference quantum methods for this set of expanded porphyrins, several diagnostics for static correlation have been examined (Table  2). 101 The single-excitation amplitude vector t 1 from the CCSD calculations offers up two diagnostics. The first, the T 1 diagnostic of Lee and Taylor, 102 is the vector norm of t 1 normalized by the number of correlated electrons, ||t 1 ||/N 1/2 = (t 1 T ·t 1 /N) 1/2 . The second, the D 1 diagnostic of Janssen and Nielsen, 103,104 is the square of the largest eigenvalue of the matrix product t 1 ·t 1 T . For the didactic example of a complex between ozone (strong static correlation) and octadecane (purely dynamical correlation) at very long distance, it is obvious from the structure of t 1 T ·t 1 and of t 1 ·t 1 T that D 1 will approach its value for the worst monomer (in this case, ozone), while T 1 will approach its value for the largest monomer (octadecane). For our set of expanded porphyrins, we see that T 1 varies in a rather narrow band ( Table 2) and deceptively stays below the allegedly "safe" limit of 0.02, while D 1 shows a more pronounced variation, ranging from 0.08 to 0.13. D 1 is largest for the Mobius conformers and transition states with Mobius topology, and smallest for the Huckel and twisted-Huckel structures. In a recent benchmark study of explicitly correlated coupled cluster methods for the W4-17 thermochemical benchmark, 105 we obtained a set of diagnostics for a broad spectrum of small molecules as a byproduct. Even the lowest D 1 diagnostics in the porphyrinoid set are already in the same range as O 3 (0.077), FOOF (0.087), and C 2 ( 1 Σ + g ) (0.086)all cases with strong static correlation. Despite the fairly large correlation (R 2 = 0.93) between the T 1 and D 1 diagnostics ( Figure S7), these results suggest that the T 1 diagnostic is not informative in terms of the degree of multireference character in Huckel and Mobius porphyrins.
Another diagnostic, nearly as old as CASSCF and CI calculations, that can shed light about the multireference character is 1 − C 0 2 , where C 0 is the coefficient of the reference determinant in a CASSCF or CI calculation. 106 1 − C 0 2 is intuitively understood as the weight of determinants other than the reference in the energy. Generally, the system is regarded to possess significant multireference character if the coefficient C 0 is smaller than 0.95 or the weight C 0 2 is smaller than 0.90. 107 The specific values given in Table 2 were obtained from both iterative full CI (ICE-CI) calculations with 30 electrons in 30 orbitals and CASSCF (12,12) wave functions. The 1 − C 0 2 values obtained from CASSCF (12,12) and ICE-CI(30,30) methods provide very similar information with R 2 = 0.95. Similar to D 1 , we see a much more pronounced variation in static correlation among the different π-conjugation topologies than for T 1 .
Yet a fourth is Truhlar's M diagnostic, 108 which for closedshell systems is effectively the average of the deviations from 2 and 0, respectively, of the HOMO and LUMO natural orbital occupations; therefore, one must always have 0 ≤ M ≤ 1 for closed-shell species. 1 − C 0 2 and M are statistically very similar for the problem at hand, with R 2 = 0.965 (Figure 7a). Interestingly, the relationship M vs 1 − C 0 2 shows the separation of the porphyrinoid structures into two groups that exhibit different degree of static correlation: (a) the Huckel and twisted-Huckel topologies as well as the Mobius pentaphyrins and (b) the Mobius structures of hexa-and heptaphyrins.
A fifth diagnostic to measure the importance of static correlation is I ND of Matito and co-workers, 109 as defined in eq 20 of that paper and obtained from the 30-in-30 FCI natural orbital occupations. We found it to be statistically most similar to 1 − C 0 2 for the same wave function (R 2 = 0.97), less so to M (R 2 = 0.84), and still less to T 1 (R 2 = 0.81) and D 1 (R 2 = 0.79).
Finally, there is the von Neumann correlation entropy, 110 where the n i are the natural orbital occupation numbers: S 2 is not size intensive. In the present work, we use just the HOMO−1, HOMO, LUMO, and LUMO+1 orbitals out of the 30-in-30 FCI to mitigate size disintensivity. Thus, the S 2 obtained contains basically the same information as M (R 2 = 0.995, Figure S7) and, accordingly, we have hence not considered it further. The M diagnostic, in particular, turns out to be a fair predictor for the energy difference between canonical CCSD(T) and MP2 calculations, with R 2 = 0.83 ( Figure 7b); CASSCF 1 − C 0 2 has a more modest but still respectable R 2 = 0.79, while remaining diagnostics (S 2 aside, see above) fare more poorly (R 2 = 0.65 or less). The CCSD(T) − MP2 energy differences can be seen as a pragmatic measure of the importance of higher-order correlation effects.
In addition, we have computed the fractional occupation number weighted density (FOD) as a real-space measure for The Journal of Physical Chemistry A pubs.acs.org/JPCA Article static electron correlation effects. 111 This density is obtained by performing finite-temperature DFT calculations and allows for the localization of strongly correlated electrons in molecules with a complicated electronic structure. 112 A representative example of the FOD plots for the different topologies of [28]hexaphyrin (28F, 28M 1A , and 28H) is shown in Figure S8. For the three topologies, the FOD is significantly delocalized over the entire macrocyclic ring and not localized on any particular moiety of the molecules, which would indicate that the whole system is strongly correlated. Unexpectedly, no significant differences are observed in the FOD plots for the Huckel and Mobius topologies, in contrast to other diagnostics computed from the ICE-FCI (30,30) wave functions which clearly show that the 28M 1A exhibits stronger multireference character than 28H and 28F. This prompts the question whether we can trust CCSD(T) at all for such extended π-systems. Unfortunately, higher-order correlated methods with the full complement of orbitals are not a practical option. We can, however, carry out iterative full CI calculations in a small orbital "active space" by means of the ICE-CI method, 113 and we have done so for 12-in-12, 18-in-18, 24-in-24, and 30-in-30 active spaces and then carried out CCSD(T) in the same restricted spaces for comparison. The results are collected in Table S3 and summarized in Table 3. It is clear that, at least for the relative energies of topology interconversions in expanded porphyrins, CCSD(T) closely tracks FCI for all orbital spaces (see also Figure S9). The fact that post-CCSD(T) contributions are surprisingly small might arise from a fortunate compensation of errors between the repulsive contribution of higher-order iterative triple excitations [CCSDT−CCSD(T)] and the attractive contribution of connected quadruple excitations [CCSDTQ-CCSDT]. Similar cancellations have been already observed for small molecules with strong multireference character, such as the C 2 ( 1 Σ + g ) molecule. 114,115 Then, what sets the Mobius structures's CI wave functions apart from those for the twisted-Huckel and Huckel structures? Table 4 shows the coefficients of the leading configurations in a 30-in-30 full CI expansion for the different structures involved in the topology interconversions of [28]hexaphyrin. From Table 4, it is clear that the Mobius structures have qualitatively different multireference character than the twisted-Huckel and Huckel-untwisted structures, with a prominent second component that effectively corresponds to exciting one electron each from the two quasidegenerate HOMO and HOMO−1 orbitals to the quasidegenerate LUMO and LUMO+1 orbitals. There are six unique determinants that fit this pattern, which generate three singlet and three M S = 0 triplet configuration state functions. The same pattern is seen for the TS 2A and TS 2B Mobius transition states. Effectively, a single Slater determinant is a fairly poor zero-order reference for these systems, and a 4-in-4 CASSCF would be required. This would also embrace most of the other prominent determinants in the expansions.
In contrast, despite having significant contributions from excited determinants, the HF reference determinant of the twisted-Huckel and Huckel structures as well as the transition structures (28TS 1A , 28TS 1B , and 28TS 3 ) has a coefficient of 0.895 or more. For them, a 6-in-6 CASSCF would be preferable over 4-in-4, but even a single HF determinant is still a reasonable reference unlike for the Mobius structures.
Performance of Lower-Level Electron Correlation Methods. Having assessed the reliability of the canonical CCSD(T) method for this set of Huckel and Mobius expanded  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article porphyrins, the next step involves the performance of lowerlevel wave function methods. Tables S4 and S5 collect the relative energies computed with different CCSD-based approaches and MP-based procedures, respectively, together with the mean unsigned error (MUE) and root-mean-square deviation (RMSD). For the purpose of assessing the performance of CCSD and MP-based methods against canonical CCSD(T) results, we resort to the CCSD(T)/cc-pVDZ reference energies, since the same basis set is employed in all wave function approaches. The RMSD deviations for each wave function approach with respect to canonical CCSD(T)/cc-pVDZ energies are visualized in Figure 8.
First of all, the importance of connected triple excitations is illustrated by the RMSD of 5.1 kcal mol −1 for CCSD. The breakdown of the CCSD/cc-pVDZ reference energies into the SCF, CCSD, and (T) components (Table S6) shows that the inclusion of triple excitations has a major effect on the relative energies of the Mobius topologies, which become greatly stabilized upon inclusion of higher excitations into the coupled cluster calculations. The statistical error drops to 3.9 kcal mol −1 for SCS-CCSD and 2.8 kcal mol −1 for SCS(MI)-CCSD, 42,116 or to 3.5 kcal mol −1 for the distinguishable cluster singles and doubles (DCSD) method, 43 which costs the same as CCSD. Adding spin-component scaling, i.e., SCS-DCSD, 44 improves the distinguishable cluster approach to an RMSD of  The Journal of Physical Chemistry A pubs.acs.org/JPCA Article 2.7 kcal mol −1 . As expected, the larger statistical errors of all CCSD-based approaches are found for the Mobius structures of [28]hexaphyrin and [32]heptaphyrin, which exhibit a more pronounced multireference character. These results confirm that the DCSD is less sensitive to static correlation than CCSD and that the addition of the SCS correction improves the accuracy of coupled-cluster-based approaches. Figure 9 shows the box-and-whisker plots for different CCSD-based approaches. It becomes clear that the error distribution becomes progressively narrower for the series CCSD, SCS-CCSD, and SCS-DCSD. The RMSD decreases to 3.9 and 2.7 kcal mol −1 for SCS-CCSD and SCS-DCSD, showing that spin-component scaled approaches improves the accuracy of coupled-clusterbased approaches. Among spin-component-scaled approaches, the distinguishable SCS-DCSD approach results in a narrower error distribution as compared to SCS-CCSD for the expanded porphyrin database. Since the localized orbital coupled cluster theory has recently emerged as an efficient quantum chemical method for coupled cluster calculations on large systems, 26 we have also investigated the accuracy of the DLPNO-CCSD(T 1 ) approach for these extended π-systems using two different sets of truncation thresholds, termed "NormalPNO" and "TightP-NO". 29 The energy differences with respect to canonical CCSD(T)/cc-pVDZ relative energies are collected in Table S6 and the box-and-whisker plots of the energy deviations for the different settings are shown in Figure 9. With the TightPNO setup, DLPNO-CCSD(T 1 )/cc-pVDZ delivers relative energies with an RMSD of 1.3 kcal mol −1 with respect to canonical CCSD(T)/cc-pVDZ calculations. Loosening of the thresholds results in a slightly larger error distribution and the RMSD increases to 1.7 kcal mol −1 for NormalPNO. Although the statistical errors increased with the default PNO settings, the wall time is reduced by 75−80%.
Overall, similar trends are provided by the local pair natural orbital coupled-cluster theory, but the errors for the Mobius structures reach up to 4.7 kcal mol −1 due to their more pronounced multireference character. Indeed, the M diagnostic turns out to be a good predictor for the energy difference between canonical CCSD(T) and localized DLPNO-CCSD-(T) calculations, with R 2 = 0.89 ( Figure S10). For these reasons, we conclude that canonical coupled-cluster theory is needed for benchmarking purposes in the case of Huckel− Mobius interconversions in expanded porphyrins. The assessment of the performance of a large variety of localized coupled cluster methods for these challenging structures is currently underway. 30 With HF being a poor zero-order wave function for the Mobius species, it is not surprising that the statistics for normal MP2 are also not good (RMSD = 9.5 kcal mol −1 ). In contrast to the CCSD approach, MP2 overstabilizes the Mobius-type structures to a huge extent, leading to energy deviations up to 20 kcal mol −1 relative to canonical CCSD(T) method (Table  S5). Orbital-optimized Møller−Plesset perturbation (OO-MP2) theory 117 in fact performs worse than conventional MP2, leading to an RMSD of 12.0 kcal mol −1 . By contrast, spin-component scaling 36 greatly mitigates the problem, with an RMSD = 3.0 kcal mol −1 for SCS-MP2. Intriguingly, Fink's S2-MP2, 37 despite its sounder theoretical foundation, actually performs worse at RMSD = 7.9 kcal mol −1 . A possible reason The Journal of Physical Chemistry A pubs.acs.org/JPCA Article for the bad performance of S2-MP2 is the high coefficient for the same-spin MP2 correlation energy (c ss = 0.75), which seems to have a negative impact on the performance of spincomponent-scaled approaches for the description of these topology interconversions. SOS-MP2, 118 however, which completely neglects same-spin correlation, achieves a stunning RMSD = 0.8 kcal mol −1 . This behavior is not specific to this subset of systems. For the MOBH35 organometallic transition state benchmark, 119 standard MP2 has an RMSD = 6.0 kcal mol −1 , compared to 2.8 and 2.5 kcal mol −1 , respectively, for SCS-MP2 and SOS-MP2 ( Table 2 in ref 74). This is not the case for problems less driven by static correlation, as shown for the very extensive GMTKN55 (General Main-group Thermochemistry, Kinetics, and Noncovalent interactions, 55 problem types) benchmark of Grimme, Goerigk, and co-workers, 120 in which SOS-MP2 is clearly inferior to standard MP2 and SCS-MP2. 74 Inclusion of damped third-order contributions, however, appears to improve on the shortcomings of same-spin MP2 for strongly correlated systems, as seen from the RMSDs of SCS-MP3 121 (0.9 kcal mol −1 ) and to MP2.X 40 (0.8 kcal mol −1 ).
Performance of Density Functional Methods Belonging to Different Rungs of Jacob's Ladder. Since DFT remains the workhorse for modeling expanded porphyrins, we have investigated the performance of a variety of exchange− correlation functionals from different rungs of Jacob's ladder. 45,93 Tables S8−S15 collect all the relative energies for our expanded porphyrin database computed with all the DFT methods. RMSDs for each functional relative to canonical CCSD(T)/CBS energies are visualized in Figure 10.
Overall, Figure 10 reveals a poor performance for GGAs, meta-GGAs, GGA hybrids, range-separated hybrids, and double hybrids for describing the thermochemistry and kinetics of topological interconversions in expanded porphyrins. Among all the functionals tested in this work, the rangeseparated double hybrids, such as ωB97M(2) 71 and ωB2GP-PLYP, 73 are the most reliable approaches for describing the relative energies of Huckel−Mobius interconversions, providing RMSDs within chemical accuracy. In contrast, rangeseparated hybrids perform badly for these topological interconversions, including the promising ωB97MV functional with nonlocal correlation. 65 In a recent benchmark study of 200 density functionals, 93 ωB97MV emerges as the best choice for nearly all interaction types with the exception of systems highly sensitive to the self-interaction/delocalization errors. For our set of highly delocalized π-systems, 20 the Minnesota meta-GGA functionals with high fractions of exact exchange (M06-2X and M08-HX) perform reasonably well, outperforming more robust and significantly less empirically parametrized functionals within their class (SCAN0-D3 and PW6B95-D3). It is important to note that the conclusions regarding the performance of the different type of functionals remain the same when the relative errors per functional is analyzed separately for reaction energies and barrier heights ( Figure  S11).
GGA functionals, like BLYP, PBE, or revPBE, exhibit very poor performance, similarly to the meta-GGA functionals. Global hybrids outperform GGA and meta-GGA functionals, but still they are not accurate enough to describe these Huckel−Mobius topological interconversions. Dispersion corrections significantly improve the performance of the standard functionals, in contrast to previous findings ( Figure  11). 18,23 For nearly all the tested functionals across the rungs of Jacob's Ladder, we found a significant improvement for describing the relative energies of these challenging systems upon addition of atom-pairwise dispersion corrections. The improvement of performance upon addition of the D3(BJ) dispersion can be linked to the variety of noncovalent interactions present in the different conformations, ranging from hydrogen bonding to π···π stacking interactions ( Figure  S13). For the M06-2X functional, we found that the D3(0) yields only 0.2 kcal mol −1 improvement with respect to the  (Table S16). These are the prefactors s 6 and s 8 for the R −6 and R −8 terms, respectively, and the damping function turnover parameters a 1 and a 2 (eq 7).
where the C 6 AB and C 8 AB are dispersion coefficients specific to the atom pair and°= R C C / AB 8 AB 6 AB . The larger a 1 and a 2 become, the further out to the periphery is the dispersion correction pushed. Conversely, the larger s 8 becomes (for typical values of a 1 and a 2 ), the more prominent mediumdistance noncovalent interactions become, to the extent that they may interfere with static correlation effects. For typical GGAs and hybrid GGAs, the large values of s 8 , combined with typical values of a 1 and a 2 , ensure large D3BJ contributions. The role of the dispersion is perhaps best understood as favoring the spatially more compact twisted-Huckel and Mobius structures over the spatially more spreadout untwisted Huckel topology. For SCAN, SCAN0, and M06-2X, s 8 = 0, so the effect of D3BJ is mitigated.
An alternative way to improve the accuracy of GGA hybrids is introducing additional ingredients into the density functional form, as in meta-GGAs hybrids. Within this group, we have considered several Minnesota density functionals with variable fraction of HF exchange, highly parametrized functionals, and extreme sensitivity to the integration grid, 122 but showing good performance for main group thermochemistry and kinetics. 93,123,124 However, none of them are state-of-the art for noncovalent interaction data sets. 120,124 The fraction of exact exchange across these functionals varies from 27% (M06) to 44% (MN15) to 52% (M08-HX) and 54% (M06-2X). M06-2X and M08-HX exhibit the best performance with an RMSD about 1.7 kcal mol −1 , followed by MN15 (RMSD = 2.1 kcal mol −1 ). The lower performance of MN15 with respect to M06-2X for expanded porphyrins is rather unexpected considering the partial multireference character of Mobius structures. MN15 was specifically developed for preserving the performance of M06-2X for single-reference main-group thermochemistry, while improving its description of multireference cases. 67 The fact that more exact exchange reduces the errors reinforce the idea that the delocalization error represents an additional issue in Huckel and Mobius expanded porphyrins. 125 Range-separated hybrid functionals show significantly poorer performance compared to the meta-GGA hybrids, except for CAM-B3LYP-D3 (RMSD = 1.8 kcal mol −1 ), as pointed out in previous benchmark studies on expanded porphyrins. 23,99 We recently proved that M06-2X and CAM-B3LYP exhibit very similar degrees of π-electron delocalization in expanded porphyrins. 20 Disappointingly, the promising range-separated hybrid ωB97M-V with nonlocal correlation, 65 designed using a combinatorial approach and containing only 12 parameters, clearly underperforms the older ωB97XD for topology interconversions in expanded porphyrins. Similarly, the range-separated GGA hybrid with 10 parameters ωB97XV exhibits larger errors than the original ωB97XD.
Regarding DSD double hybrids, one might expect these functionals to suffer from the same issues as the MP2 method. Indeed, the original DSD-PBEP86 performs poorly at RMSD = 6.4 kcal mol −1 , but the error drops to 4.7 kcal mol −1 upon addition of the D3(BJ)-dispersion correction. However, the entire DSD family of functionals has recently been reparametrized 126 against the much larger and more chemically diverse GMTKN55 database, 120 which led to a substantial change in parameters and improved performance across the board. Our revDSD-PBEP86-NL functional indeed reduces the RMSD to 2.6 kcal mol −1 , presumably due to the much smaller coefficient for same-spin MP2-like correlation in revDSD The Journal of Physical Chemistry A pubs.acs.org/JPCA Article compared to DSD. As we have seen earlier for transition metals reaction energies, 127 the revDSD parametrizations are not only more accurate than the original DSD but more robust as well. Furthermore, the comparison of spin-component (DSD) and spin-opposite (DOD) scaled double hybrids ( Figure S12) again shows that elimination of the same-spin correlation in DOD-type functionals results in a significant improvement over DSD-type functionals, except for the DSD-SCAN x -D3 functionals. Consequently, the best performers, among the spin-component scaled dispersion-corrected double hybrids, are the revDOD-PBEB95-D3 (RMSD = 1.8 kcal mol −1 ), revDOD-PBEP86-NL (RMSD = 1.9 kcal mol −1 ), and especially the nonempirical SOS0-PBE0-2-D3(BJ) 78 (RMSD = 0.9 kcal mol −1 ). Finally, the RMSDs of selected hybrids and double hybrid functionals for each family of expanded porphyrin are carefully analyzed. Figure 12 summarizes the errors for the B3LYP-D3(BJ) (GGA hybrid), M06-2X-D3(0) (meta-GGA hybrid), CAM-B3LYP-D3(BJ) (range-separated GGA hybrid), re-vDOD-PBEP86-NL (double hybrid), and ωB97M(2) (rangeseparated double hybrid) for N-fused [24]pentaphyrin, [28]hexaphyrin, and [32]heptaphyrin. The performance of the functionals is clearly different for each type of macrocycle, and only ωB97M(2) shows errors within 1 kcal mol −1 for all three expanded porphyrins considered here. Not surprisingly, most of the functionals perform worse for the relative energies of [28]hexaphyrin, which contain a larger number of Mobius structures with a more challenging electronic structure. By contrast, the relative energies of Huckel and Mobius structures of the pentapyrrolic macrocycle are relatively well described by the different functionals except B3LYP-D3, which can be expected from the minor degree of static correlation in the smallest macrocycle. In this case, the five functionals predict the same global minimum (24H b ) and similar activation barriers for the Huckel and Mobius interconversion, in good agreement with experimental studies. 128 In the case of [28]hexaphyrin, only ωB97M(2) provides reasonably accurate energies relative to CCSD(T)/CBS, followed by the M06-2X-D3(0) functional. According to canonical CCSD(T)/CBS calculations, unsubstituted [28]hexaphyrin can adopt three different topologies with small energy differences (28H, 28M, and 28F), explaining why mesoaryl substituted [28]hexaphyrin exists in solution as an equilibrium of several equivalent twisted Mobius conformations and a Huckel rectangular conformation. 129 According to CCSD(T)/CBS energies, the Mobius topology 28M corresponds to the thermodynamically more stable conformation, in line with the experimental findings from spectroscopic investigations. 130 Here, it is important to note that M06-2X-D3(0) and CAM-B3LYP-D3(0) both predict the Huckel untwisted topology 28H to be the global minimum.
For the flexible [32]heptaphyrin, the canonical CCSD(T) calculations reveal the predominance of the twisted-Huckel conformer (32F), in agreement with the experimental findings. 21,22 The five functionals predict the twisted-Huckel as most stable, but important energy differences between functionals appear for the relative energies of the Mobius topologies (32M a and 32M b ), belonging to the group of structures with larger multireference character. The difficulty to describe the relative energies of the interconversion between the twisted-Huckel and Mobius topologies in [32]heptaphyrin arise from the partial multireference character of the Mobius topologies as well as the delocalization error. In this regard, we observe that the statistical errors for B3LYP arises mainly from the overstabilization of the Mobius topologies with respect to the Huckel ones. This overstabilization can be traced to the delocalization error, since B3LYP exaggerates the degree of πelectron delocalization in aromatic macrocycles. 20 This overstabilization is mitigated to a great extent by means of the D3(BJ) dispersion. Overall, these results highlight the critical role of the exchange−correlation functional for studying the conformational preferences and topology interconversions in expanded porphyrins.

■ CONCLUSIONS
The ultimate goal of this study is to identify affordable wave function methods and density functionals for the prediction of accurate relative energies for Huckel-Mobius topology interconversions in expanded porphyrins. Such topology interconversions can be exploited in the development of molecular switches for multiple applications; computational chemistry represent a powerful tool toward the design of efficient switches. DFT remains the workhorse for modeling such extended macrocycles, but these macrocycles are challenging from an electronic structure point of view. In this work, we demonstrate that the degree of static correlation largely varies between the different topologies, as revealed by several diagnostics for static correlation. Although the T 1 diagnostic values stay deceptively below 0.02 for all the investigated systems, D 1 , M, I ND , S 2 , and 1 − C 0 2 diagnostics show a more pronounced variation between Huckel and Mobius π-systems. According to these tests, the selected porphyrinoid structures can be split into two groups as a function of the degree of static correlation, with the Mobius structures of hexa-and heptaphyrins being more "troublesome" for single-reference methods.
Accordingly, we assess the performance of a variety of wave function methods and density functionals for thermochemistry and kinetics of topology interconversions across a wide range of macrocycles. As a reference, fully canonical CCSD(T) energies were used; post-CCSD(T) correlation effects were estimated by comparing CCSD(T) with approximate full configuration interaction (ICE-CI) in an expanding sequence of active orbital spaces. According to our results, the contributions of post-CCSD(T) correlation and basis set expansion to relative energies are minor. Hence, our CCSD(T) reference values are reasonably well-converged in both 1particle and n-particle spaces.
Regarding lower-cost wave function methods, we found that conventional CCSD, MP2, and MP3 perform very poorly, while spin-component scaling approaches like SOS-MP2 or SCS-MP3 reproduce the CCSD(T) relative energies within chemical accuracy. Third-order contributions become important to mitigate the shortcomings of same-spin MP2 for strongly correlated systems; lower-cost MP2.X approaches CCSD(T) quality for the description of topology interconversions.
The accurate description of the relative energies of our expanded porphyrin database becomes difficult for most of the density functionals tested in this work, in particular for the largest macrocycles. The largest errors are found for the Mobius topologies, which exhibit a stronger multireference character. Inclusion of atom-pairwise dispersion corrections and a large percentage of exact HF exchange improve the performance of the exchange−correlation functionals. Overall, the safest choice to describe the energetic profile of Huckel− The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript.