Relevance of the DFT method to study expanded porphyrins with different topologies

: Meso-aryl expanded porphyrins present a structural versatility that allows them to achieve different topologies with distinct aromaticities. Several studies appeared in the literature studying these topological switches from an experimental and theoretical point of view. Most of these publications include density functional theory calculations, being the B3LYP the most used methodology. In this work, we show that the selection of the functional has a critical role on the geometric, energetic


1) INTRODUCTION
From the seminal work of Heilbronner, [1] the concept of Möbius aromaticity attracted considerable attention from an experimental and theoretical point of view. [2]Heilbronner, based on Hückel molecular orbital theory, predicted that singlet annulenes with 4n π electrons would be aromatic systems in twisted conformations (with an odd number of half-twists), where the p orbitals lie on the surface of a Möbius strip; whereas singlet [4n+2] Möbius molecules would be antiaromatic.Following the work of Calugareanu, White, and Fuller, [3] Rappaport and Rzepa [4] proposed that the topology of a band can be partitioned into two chiral indices, the writhe and twist, whose sum gives an overall integer linking number, Lk.In this formulation, an even linking number Lk follows the 4n+2 electron rule for the closedshell ground state, while systems with odd values for Lk follows instead the 4n electron rule.
In the last two decades, annulenes and expanded porphyrins have shown a conformational flexibility that allows them to achieve original motifs (for instance Hückel planar, single twisted Möbius, figure-eight, dumbbell, and triply twisted Möbius) with specific aromaticities, magnetic, and electric properties. [5]The meso-aryl expanded porphyrins have provided an effective platform for the synthesis of Möbius topologies using chemical and physical modifications such as temperature control, [6,7] solvent or protonation, [8] metal coordination, [9] phosphorus [10] or silicon insertion, [11] or regioselective fusion of cyclic subunits. [12][15] It is pertinent to note that the largest part of the experimental works published in the literature includes computational results, most of them using the hybrid B3LYP functional as the default methodology.
In a previous work, [14] we benchmarked a variety of computational methods to describe the topological switch between the Hückel planar and the singly twisted Möbius structure for a model of [28]hexaphyrin (1.1.1.1.1.1.1),where the meso-substituents were hydrogen atoms.The calculations were performed with the B3LYP, BH&HLYP, CAM-B3LYP, M05-2X, and MP2 methodologies using the 6-31G and 6-311G(d,p) basis sets.For benchmarking purposes, single-point energies were also calculated, at the optimized structures, employing the CCSD(T)/6-31G level of theory.We concluded that the increment of the size and flexibility of the basis set from 6-31G to 6-311G(d,p) shows a small effect in the relative stabilities and energy barriers for these conformational switches.A very important finding of that work is that the selection of the computational method is crucial for obtaining a correct energetic and geometric description of these systems, and our results showed that CAM-B3LYP and M05-2X are the functionals that provide the most reliable results.Moreover, it is worth noting that B3LYP and MP2 methodologies exaggerate the delocalization of the Möbius systems and overstabilize these topologies with respect to their Hückel analogues. [14]e relevance of these expanded porphyrins as a new promising class of molecules for the creation of topological switches motivated us to perform a more systematic study for two challenging distortions using eleven density functional theory (DFT) functionals and performing local pair natural orbital coupled cluster, DLPNO-CCSD(T), calculations for benchmarking purposes.To this aim, we have studied the figure-eight (E32) → Möbius (M32) → Hückel (H32) switch of [32]-heptaphyrin(1.1.1.1.1.1.1)(see Scheme 1) and the interconversion between the conjugation paths of [26]-hexaphyrin(1.1.1.1.1.1)in the Hückel geometry (see Scheme 2).In the results and discussion section, we show that for these two considered cases the DFT methods employed present: i) discrepancies in the energetic differences between topologies larger than 10 kcal/mol; and ii) differences in the number and/or nature of the optimized stationary points of the PES (i.e.some functionals determine a spurious minimum [16] instead of a transition state).When these important divergences between the results obtained with different DFT methods occur, it results essential to perform high level ab initio calculations to judge the reliability of the PES obtained.
(Insert Schemes 1 and 2 around here)

2) COMPUTATIONAL METHODS
The geometries for all stationary points have been optimized using the 6-31G and 6-311G(d,p) basis sets [17] employing eleven DFT methods (B3LYP, [18] BPE, [19] BP86, [20] M06L, [21] TPSSh, [22] M05-2X, [23] M06-2X, [24] BH&HLYP, [25] CAM-B3LYP, [26] BMK, [27] and WB97XD [28] ).At each stationary point we have carried out harmonic frequency calculations to ensure the nature of the corresponding stationary point (minimum or transition state), and to provide the zero-point vibrational energy (ZPE) and the thermal and entropic contributions to the enthalpy and free energy at T=298 K.In the case of the Hückel, Möbius, and figure-eight topologies of [32]-heptaphyrin with hydrogen atom as meso-substituent, we have used as initial structures the geometries reported by Alonso et.al. [15] Otherwise, in the second illustrative example, the Xray Hückel structure of the meso-hexakis (2,4,6-trifluorophenyl)  [26]-hexaphyrin(1.1.1.1.1.1)has been used as the initial geometry. [6]For five selected functionals (B3LYP, PBE, BP86, CAM-B3LYP, and BMK) the dispersion effects were also included by using the D3-Grimme's dispersion correction [29] with Becke-Johnson damping. [30]For benchmarking purposes, single-point energy calculations at the M05-2X/6-311G(d,p) optimized geometries have been evaluated at DLPNO-CCSD(T) level with VDZ [31] and cc-pVTZ [32] basis sets and using resolution of identity and the corresponding auxiliary basis reference, [33] VDZ/C and cc-pVTZ/C.DLPNO-CCSD(T) is an efficient quantum chemical method developed by Neese and collaborators, [34] which uses the concept of the local pair natural orbitals (LPNO), [35] that allows for coupled cluster calculations including single-, double-, and perturbative triple-excitations on large-scale chemical applications.In order to check the reliability of single-determinant-based methods, we have computed the T1 diagnostic [36] of the DLPNO-CCSD and CCSD wave functions to analyze their multireference character.In all studied systems, the T1 diagnostic values are smaller than 0.015, which confirms the appropriateness of the single-reference wave-functions used in this work.The aromaticity of selected structures has been evaluated using a geometric (harmonic oscillator model of aromaticity, HOMA) [37] and magnetic (nucleus-independent chemical shift, NICS, [38] using the GIAO method) criteria at the B3LYP and M05-2X levels using the 6-311G(d,p) basis set.For some specific structures, the visualization of non-covalent interactions (NCI) has been performed using the NCIPLOT program. [39]All geometry optimizations, frequencies, and energies of the DFT methods were calculated using the Gaussian 09 program package, [40] while the Orca program [41] was used to calculate the DLPNO-CCSD and DLPNO-CCSD(T) single-point energies.

3) RESULTS AND DISCUSSIONS
In a previous work Alonso et.al. [15] reported the relative energies of the figure-eight and Möbius conformers of three different meso-substituted [32]-heptaphyrins using B3LYP, B3LYP-D, and M06-2X functionals.The meso-substituents studied were H, C6F5, and C6H3Cl2.They concluded that B3LYP methodology shows the best overall performance for describing the geometries and thermochemistry of these topological switches.Moreover, the relative energies obtained by Alonso et.al. [15] presented important discrepancies among the different DFT functionals, e.g. 10 kcal•mol -1 between M06-2X and B3LYP methodologies, which motivated us to consider this system as a good candidate to test the performance of different theoretical approaches.
Table 1 gathers the relative electronic energies, electronic energies plus zero-point vibrational energies, enthalpies, and free energies for the figure-eight (E32) → Möbius (M32) → Hückel (H32) switch of hydrogen meso-substituent [32]-heptaphyrin (see Scheme 1 and Figure 1) calculated using the 6-311G(d,p) basis set and eleven DFT functionals and DLPNO-CCSD and DLPNO-CCSD(T) methodologies with the cc-pVTZ basis set (see Table S1 of the supporting information for the results at DFT/6-31G and CCSD(T)/6-31G levels).The results obtained at the DLPNO-CCSD(T)/cc-pVTZ//M05-2X/6-311G(d,p) level indicate that the figure-eight conformation is the most stable and the Möbius and Hückel topologies lie 14.2 and 32.2 kcal•mol -1 , respectively, above E32 (14.7 and 34.2 kcal•mol -1 at the CCSD(T)/6-31G//M05-2X/6-311G(d,p) level).This expanded porphyrin shows a 4n (n=8) π number of electrons and E32, M32, and H32 conformations present twisted π surface, half-twisted π surface, and planar topologies, respectively (notice that in the case of hydrogen as meso-substituent the figure eight structure has almost no twist).According to the Hückel-Heilbronner rules, [1] the E32 and H32 (M32) structures must show antiaromatic (aromatic) character (in the particular case of E32, it presents almost no twist and the applicability of these rules is questionable).The fact that the figure-eight structure shows the most stable conformation can be rationalized in terms of its more effective network of intramolecular NH•••N hydrogen bond interactions between imine-type and amine-type pyrroles. [15]th respect to DLPNO-CCSD(T)/cc-pVTZ//M05-2X/6-311G(d,p) benchmark values, the methodologies that show the most important discrepancies (in increasing order of the maximum error) are TPSSh, PBE, B3LYP, and BP86.For instance, B3LYP (BP86) predicts that M32 and H32 structures are, with respect to E32, 9.4 (11.7) and 5.8 (4.8) kcal•mol -1 , respectively, more stable than the DLPNO-CCSD(T) results.It is important to remark that an energetic difference of 9.4 (11.7) kcal•mol -1 represents a decrease of ca.66% (82%) on the DLPNO-CCSD(T) relative stability of M32 with respect to E32, pointing out the key relevance of the computational method for a correct energetic description of these topological switches.An intermediate behavior with respect to the DLPNO-CCSD(T) calculations is obtained from CAM-B3LYP, M06-2X, BH&HLYP, and M06L functionals (in increasing order of the maximum error), whose differences with respect to DLPNO-CCSD(T) lie in the range of 2.3-5.7 kcal•mol -1 .The best performance is obtained for the BMK and M05-2X functionals, with energetic differences with respect to DLPNO-CCSD(T) method smaller than 3.1 and 2.5 kcal•mol -1 , respectively.
A similar scenario has been observed in the case of the annulenes, for which it has been reported that some DFT and MP2 methodologies exaggerate the delocalization of conjugated systems and overstabilize their structures, while the HF level underestimates this delocalization in aromatic molecules. [42]The present results (and the obtained in the following studied case, vide infra) are very relevant, because B3LYP functional is the most popular methodology used in the computational characterization of these expanded porphyrins and consequently, the reliability of some of the theoretical data reported in the literature could be questioned.The WB97XD functional has not been considered, because the optimization process of the Hückel structure yields a Möbius topology (indicated as M232 in Table 1), resulting in a different Möbius conformation than M32.The selection of the M05-2X geometries for the benchmark DLPNO-CCSD(T) single-point energy calculations is based on previous studies of annulenes [42] and our previous work of expanded porphyrins. [14]nsert Table 1 and Figure 1 around here) We have also analyzed the effects of the basis set and geometry selection.The DFT functionals show differences between 6-311G(d,p) and 6-31G basis sets of around 2 (4) kcal•mol -1 for the relative stability of M32 (H32) with respect to E32, see Tables 1 and S1.These values represent ca.10-20% of their relative stability values, which is a not negligible contribution, although it is small in comparison to the relevance of the computational method (vide supra).Moreover, the effect of the geometry has been evaluated doing single-point energy calculations at all DFT methodologies using the M05-2X/6-31G optimized geometries (see Table S1).The differences of the relative stability values between the optimized DFT/6-31G and DFT/6-31G//M05-2X/6-31G geometries are smaller than 2 kcal•mol -1 .These results indicate that the effect of the DFT methodology in predicting the optimized geometry is nonnegligible, although less relevant than its effect in the electronic energy.
The good performance of the M05-2X functional (and in a minor degree M06-2X), which is one of the functionals that that better includes the dispersion, motivates us to explicitly analyze the role of the dispersion correction.For this reason, in Table 2 it is collected the results of the dispersion-corrected DFT methods (it has been only considered the functionals for which D3-Grimme's dispersion correction with Becke-Johnson damping has been implemented in the Gaussian 09 program package).The inclusion of dispersion correction induces a strong increment of the energy difference of M32 and H32 with respect to E32; e.g. at BP86 (BP86-GD3BJ) level M32 and H32 lie 2.5 (14.1) and 27.4 (45.2) kcal•mol -1 , respectively, above the figure eight structure.The important energetic differences between dispersion-corrected and -uncorrected DFT methods vary from 5.2 to 17.9 kcal•mol -1 .With respect to DLPNO-CCSD(T)/cc-pVTZ//M05-2X/6-311G(d,p) benchmark values, the inclusion of the dispersion correction does not lead to an overall improvement of the relative energies.The root-mean-square deviations (RMSD) of DFT and DFT-D3BJ values with respect to DLPNO-CCSD(T) results are 6.7 and 6.8 kcal/mol, respectively.The best performance is observed for B3LYP-GD3BJ (RMSD = 5.2 kcal/mol), which shows a slight better agreement with respect to DLPNO-CCSD(T) results than B3LYP (RMSD = 7.8 kcal/mol).On the other hand, the good behavior of the BMK functional (RMSD = 2.3 kcal•mol -1 ) is lost with the consideration of the dispersion correction (RMSD = 8.5 kcal•mol -1 ).The NCI isosurfaces of M32, H32, and E32 are reported in the Figure S1 and they indicate that as one can expect E32 structure shows the most relevant π-π stacking interactions.Lastly, these results allow concluding that the good performance of M05-2X and BMK functionals (and in minor degree CAM-B3LYP, M06-2X, BH&HLYP, and M06L) is probably due to their large percentage of exact Hartree-Fock exchange.

(Insert Table 2 around here)
Another aspect that can generate misleading conclusions is the comparison between the DFToptimized geometries and X-ray structures, which is a frequently used argument in the literature to validate the DFT results for these expanded porphyrins.To give more insight about this point, we will follow the analysis performed by Schleyer and co-workers [43] for the structure of [18]-annulene.The Xray structure of [18]-annulene [44] was determined to have a D6h symmetry, which was also confirmed by MP2 and various DFT calculations. [45]urthermore, it was found that only the DFT C2 geometry gave acceptable computed shifts. [43]The explanation of this discrepancy is that the real minimum structure of [18]-annulene is a C2 geometry, and the D6h structure is a low energy transition state joining the two equivalent C2 minima.The CCSD(T) calculations support also this conclusion. [43]On the contrary, the degree of delocalization of the D6h geometries is overestimated by various DFT and MP2 methodologies over-stabilizing such D6h isomers.Only the hybrid density functional methods with a larger HF contribution lead to the correct C2 minima.Schleyer and co-workers [43] explained the presumed disagreement between CCSD(T) results and X-ray structures in terms of the possible existence of static and dynamic disorders of the X-ray structures.Static disorders are a superposition of lower symmetry structures leading to apparent higher symmetry.A dynamic disorder is given when several isoenergetic conformations are linked by almost barrierless transition states which provokes that in the X-ray investigation, the C-C distances are averaged over time and over a large number of unit cells during the diffraction process.This phenomenon results in an incorrect experimental assignation of a D6h symmetry for the geometrical structure of [18]-annulene. [43,47] he link between the [18]-annulene and expanded porphyrins is direct, indicating that a cautious comparison between DFT and the X-ray structures is required, especially in conjugated systems with almost isoenergetic conformations.Nevertheless, it should be pointed out that a careful analysis should be performed in the comparison between theoretical and experimental X-ray structures.For instance, the crystal packing effects can exhibit a strong influence on the topology and the more planar Hückel structures usually are favored.
An additional example of malfunction of the B3LYP methodology predicting Möbius topologies is the case of the cyclic 8π electrons [9] annulene cation, C9H9 + .In 1971, it was reported its first experimental evidence such a short-lived intermediate. [48]Two decades after, Mauksch et al. [49] carried out a theoretical study of C9H9 + and they concluded that the aromatic Möbius conformation is the most stable structure.The geometry optimizations and single point calculations were performed at B3LYP/6-311+G(d,p) and CCSD(T)/DZP levels, respectively.In 2009, Herges and co-workers [50] presented highlevel coupled-cluster calculations including corrections for solvent and vibrational contributions and they predict that the Hückel and the Möbius isomer of C9H9 + are quasi-degenerate (the energy difference is only 0.04 kcal•mol -1 ).Moreover, this work also contains the UV/Vis spectra of the [9]  annulene cation measured in laser flash photolysis, which concludes that the Hückel structure is the most stable conformation since it is the unique conformer detected, although it was considered that the Möbius structure could exist in small concentrations below the detection limit.
Another peculiar case for which the selection of the DFT method used for the calculations of these Hückel-to-Möbius topological switches plays a critical role is the Hückel structure of [26]-hexaphyrin.
Initially, we have studied the [26]-hexaphyrin with hydrogen atoms as meso-substituents, see Scheme 2 and Table 3.For this model, the B3LYP, PBE, BP86, M06L, and TPSSh functionals predict that the Hückel structure presents a spurious Ci minimum geometry, H26(Ci), which shows an important bond equalized structure.On the other hand, the remaining six DFT methods (M05-2X, M06-2X, BH&HLYP, CAM-B3LYP, BMK, and WB97XD) predict that the stationary point with Ci symmetry is a transition state TS26(Ci), that joins two equivalent C1 minima, which display a longer bond alternating structure and represent the correct Hückel structures, H26.In Scheme 2, it is pertinent to note that H26 and TS26(Ci) are not planar geometries.Our DLPNO-CCSD(T)/cc-pVTZ//M05-2X/6-311G(d,p) benchmark calculations validate this description and TS26(Ci) conformation lies ca. 5 kcal•mol -1 above H26.Then, we can conclude that some DFT methods exaggerate the bond equalization of the Hückel structure of [26]-hexaphyrin and overstabilize its structure.This overstabilization converts the conjugated transition state geometry to a spurious minimum with the corresponding disappearance in the PES of the correct Hückel minima structures.
The computed relative energies of H26(Ci), H26, and TS26(Ci) using the eleven DFT methods and DLPNO-CCSD and DLPNO-CCSD(T) methodologies are collected in Table 3.For the sake of completeness, we have also included the energetic results for the Möbius structure (M26), which according to the DLPNO-CCSD(T)/cc-pVTZ//M05-2X/6-311G(d,p) method, lies 8.6 kcal•mol -1 above the Hückel conformation.From the six DFT functionals that correctly describe the bond alternating Hückel structure, the CAM-B3LYP and WB97XD are the methodologies that show the closest results to the CCSD(T) benchmarked values, with energetic differences respect to CCSD(T) ca. 1 kcal•mol -1 .The performance of the BH&HLYP, M05-2X, and M06-2X methods is also quite acceptable, with energetic differences respect to DLPNO-CCSD(T) ca. 2 kcal•mol -1 .Finally, BMK functional shows an incorrect description of the TS26(Ci) energetic barrier, which is 4.8 kcal•mol -1 lower than the DLPNO-CCSD(T) barrier.Moreover, it is worth noting that for these geometrical configurations the London dispersion GD3BJ correction has a small effect on the energies and the dispersion -corrected and -uncorrected methods lead to the same spurious, H26(Ci), or correct, H26, minima (see Tables 3 and 4).This small effect of the dispersion is also corroborated by the H26 NCI isosurface (see Figure S2).
(Insert Tables 3 and 4 around here) In a second step, we have considered a [26]-hexaphyrin model with 2,4,6-trifluorophenyl as mesosubstituent (see Figure 2), which Hückel X-ray structure is available. [6]For this chemical species we have calculated their Hückel, H26F or H26F(Ci-like), Möbius, M26F, and the transition state that joins the two Hückel structures transition state, TS26F(Ci-like), structures (the subscript like is to point out that these structures are equivalent to the Ci geometries with the hydrogen as meso-substituent, although the 2,4,6trifluorophenyl meso-substituents break the Ci symmetry).To reduce the computational effort, we have focused the analysis to only six different functionals (dispersion-corrected and -uncorrected B3LYP, PBE, BP86, M05-2X, BH&HLYP, and CAM-B3LYP, see Tables 5 and 6) and the DLPNO-CCSD(T) calculations were performed using a smaller basis set (VDZ).In analogy to the results obtained with the hydrogen as meso-substituent, B3LYP, PBE, and BP86 methodologies locate the spurious H26F(Ci-like) minima instead of the correct TS26F(Ci-like) structure.Furthermore, there are important energetic differences of ca.6 kcal•mol -1 between the Hückel-Möbius relative energies obtained with these three functionals and their counterparts computed with the M05-2X, BH&HLYP, and CAM-B3LYP methods.The former (later) group of functionals estimates that M26F lies in the range of 11.9-14.1 (6.1-8.6)kcal•mol -1 above the Hückel conformation energy.Moreover, the DPLNO-CCSD(T)/VDZ//M05-2X/6-311G(d,p) results are again in a better agreement with the M05-2X, BH&HLYP, and CAM-B3LYP methodologies than B3LYP, PBE, and BP86.The wrong description of these topological switches performed by some DFT functionals has important implications in the analysis and rationalization of molecules that have been synthesized and characterized experimentally.In the [26]-hexaphyrin model with 2,4,6-trifluorophenyl as meso-substituent, one could expect more relevant London dispersion interactions between the meso-substituents and between the porphyrin ring and the 2,4,6-trifluorophenyl, which are confirmed by the NCI isosurfaces (see Figure S3).However, the results displayed in Tables 5 and 6 do not show large differences between dispersioncorrected and -uncorrected results.For this chemical system, the best performance is given by the CAM-B3LYP-GD3BJ methodology.
(Insert Tables 5 and 6 and Figure 2 around here) In accordance with previous systems and the [18]-annulene, the explanation of these discrepancies can be found through a prompt analysis of the H26F, H26F(Ci-like), and TS26F(Ci-like) bond distances of the expanded porphyrin ring (see Figures S4 and S5 for more details of the bond distances).
The H26F(Ci-like) determined with B3LYP, PBE, or BP86 methodologies has an important bond equalized structure, e.g. the bond length alternation (BLA) shows a value of 0.033 Å at B3LYP level, see Table 4. On the other hand, with the M05-2X, BH&HLYP, and CAM-B3LYP methods, H26F Hückel geometries present a more bond alternating structure, e.g. a BLA value of 0.065 Å is obtained at M05-2X level.Then, we can conclude that the B3LYP, PBE, and BP86 functionals overestimate the bond equalization of H26F(Ci-like) structures, which must be considered as spurious minima, and a more correct description is provided by the M05-2X, BH&HLYP, and CAM-B3LYP methodologies.Analyzing in detail the bond distances of H26F and H26F(Ci-like) structures, one can see that the C-C bond distances adjacent to the meso-substituent are the main responsible for these BLA discrepancies.For instance, at B3LYP level the differences between these adjacent C-C bond distances in the H26F(Ci-like) structure are not larger than 0.01 Å.On the other hand, at M05-2X level in the H26F structure, these differences increase to 0.09 Å.Then, a very simple geometrical descriptor of delocalization (and potentially aromaticity) of these expanded porphyrins can be just the average difference between C-C bond distances adjacent to the meso-substituent.
In Table 7, the influence of the DFT functional used to evaluate the aromaticity using geometric and magnetic descriptors can be appreciated.The difference of the HOMA and NICS between H26F or H26F(Ci-like) and M26F structures is larger in the B3LYP method (0.33 Å and 24.5 ppm) than the M05-2X one (0.15 Å and 11 ppm).Moreover, the H26F X-ray and B3LYP and TS26F(Ci-like) M05-2X geometries show very similar aromatic measures (ca.0.88 Å and -17 ppm values for the HOMA and NICS, respectively) and predict a more aromatic character than the H26F M05-2X structure (0.69 Å and -9 ppm values for the HOMA and NICS, respectively).This large and "spurious" aromatic character of B3LYP H26F(Ci-like) is in agreement with their more bond equalized structure.We have selected the HOMA and NICS (evaluated at the geometrical center of the expanded porphyrin ring) as aromaticity descriptors because they are extensively used in the literature and they are quite simple to evaluate.The aim of these results is just to illustrate again the essential importance of the selection of the adequate functional for the evaluation of the geometry, energy, and magnetic properties of these topological switches.
Nevertheless, it is pertinent to note that they are not the ideal aromaticity descriptors for some of the studied structures due their non-planar geometry.Thus, the complex electronic structure of these systems probably require more elaborated methodologies for an accurate description of their aromatic, e.g. the AV1245 index developed by Matito, [51] and an additional work on this point is in progress in our laboratory.
(Insert Table 7 around here)

4) CONCLUSIONS
With the aim of bringing a rigorous framework of the DFT functional methods that can be used to obtain an accurate energetic, geometric, and magnetic description for the topological switches based on expanded porphyrins, a systematic study using eleven dispersion-corrected and -uncorrected DFT functionals has been performed.Two different molecular distortions (the figure-eight → Möbius → Hückel switch of [32]-heptaphyrin and the interconversion between the conjugation paths of [26]hexaphyrin in the Hückel conformation) have been studied as illustrative examples.Benchmarking calculations were performed for selected structures at the DLPNO-CCSD(T)/cc-pVTZ level of theory.
The obtained results highlight that the selection of a suitable DFT method for the study of the topologies of the expanded porphyrins plays a critical role.Some functionals exaggerate the bond equalization of the Hückel aromatic structures and add a fictitious stabilization to their structures.This overstabilization can represent an energetic value of 10 kcal•mol -1 , which account for ca.80% of unreal relative stability of the Hückel topologies with respect to their figure-eight counterparts.This artificial additional stabilization can also provoke the formation of spurious minima and the disappearance on the PES of the correct transition state and minima structures.The methodologies that present the worst description of these topological distortions are the B3LYP, PBE, BP86, WB97XD, and TPSSh functionals.An intermediate performance is obtained from the BMK, M06L, and BH&HLYP methods.
The methodologies that show a more consistent behavior with respect to the DLPNO-CCSD(T) results are the CAM-B3LYP, M05-2X, and M06-2X functionals.The London dispersion effects can show important effects, although overall, the corrected-dispersion functionals do not imply an improvement with respected to uncorrected-dispersion functionals.Moreover, the influence of the used methodology has also important implications in the optimized geometries and the measure of the aromaticity descriptors.
In summary, PBE, BP86, WB97XD, TPSSh, and particularly B3LYP, provide very inaccurate geometric, energetic, and magnetic descriptions of these topological switches.On the contrary, the CAM-B3LYP, M05-2X, and M06-2X are a suitable choice to perform a reliable theoretical study of the topologies of the expanded porphyrins.ACKNOWLEDGMENTS: The calculations described in this work were carried out at the Consorci de Serveis Universitaris de Catalunya (CSUC), SGI/IZO-SGIker, and DIPC.Financial support was provided by the Ministerio de Economía y Competitividad (MINECO) of Spain and FEDER (projects CTQ2016-80375-P, CTQ2014-52525-P, and Red de Excelencia Consolider CTQ2014-51912-REDC), the UPV/EHU (UFI11/22 QOSYC), the Basque Government (GV/EJ, grant IT673-13), and the Generalitat de Catalunya (Grant 2014SGR139).We would also like to thank the reviewers for their comments and helpful suggestions.

Supporting information:
Table S1 contains relative energies, energies plus zero point energies, enthalpies, and free energies of the E32 → M32 → H32 switch calculated using the 6-31G basis set.Figures S1, S2, and S3 display the NCI isosurfaces for some selected structures.Figures S4 and S5 show the H26F, H26F(Ci-like), TS26F(Ci-like), and M26F structures calculated at the B3LYP and M05-2X levels and the X-ray structure of H26F with the bond distances of the expanded porphyrins rings.Figure S3 displays the paths used to calculate the HOMA values of the H26F, H26F(Ci-like), TS26F(Ci-like), and M26F structures.Finally, Table S2 contains the Cartesian coordinates of some selected stationary points.a To see the paths used to calculate the HOMA and BLA values of the H26F, H26F(Ci-like), M26F, and TS26F(Ci-like) structures, see Figure S6 (Supporting information).b NICS evaluated at the geometrical ring center of the expanded porphyrin ring.c X-ray Hückel structure obtained from Ref. [6].

Table 5 .
Relative energies, energies plus zero-point energies, enthalpies, and free energies for the
a Imaginary vibrational frequencies of the transition states (cm -1 ) are given in square brackets.