Molecular Structure Effects on the Aggregation Motif of Porphyrins: Computational Insights

The latest advancements in semiempirical Hamiltonians have inspired new confidence for the supramolecular computational predictivity. The advanced accuracy of newly developed semiempirical methods is offered for computations of what can provide a valuable database of molecular tuning recommendations for the synthesis of new supramolecular materials. In this work the very versatile and impactful porphyrin is employed for examination of the first basic chemical tuning factors that may drive specific aggregation motifs. The 1D motifs are examined as a function of peripheral substituent steric bulk. Subsequently, a 1D‐wire versus a 3D‐square motif is investigated as a function of the metal–ligand effect. For the first time, an interesting effect of misprediction of semiempirical computations is encountered for a small class of these aggregates and is briefly examined with a conformational search analysis. These findings encourage further in silico work which is greatly required for diminishing the current discovery bottleneck in supramolecular chemistry.


Introduction
The first explorations of noncovalent interactions (NCIs) in 1873 by Johannes van der Waals [1] immediately instigated a plethora of research in a very wide spectrum of research fields, from pure physics, [2] to medicinal sciences. [3] Most notably, the examination of NCIs commenced with biological process applications, as they naturally occur, creating, as a consequence the large sector of research that is supramolecular chemistry. [4] With inspiration from these naturally occurring interactions synthetic materials known as supramolecular aggregates, have been created. These materials are unique and offer new and occasionally rare properties [5] sourced mainly from their weakly bound form and reversible interactions, as well as their desirably controlled DOI: 10.1002/adts.202100050 directionality in space. [6] The versatile nature of NCIs effectively offers a wide range of combinations of interactions between distinct mer units (the repeating noncovalently bound entities), hence an even wider spectrum of supramolecular architectures on assembly. The noncovalently bound mer units assemble in 1D stacks or in 2D arrays, or other more complex 3D architectures (Figure 1). Each structural set up offers unique properties, such as a gel phase for 1D wires, [7] or some beneficial porosity in a 3D network assembly. [8] It has also been found that fine topological interplay in the molecular structure of the mer unit yields dramatic changes in material properties and assembly dynamics, or motif, inspiring even further synthetic investigations. [9,10] During the past decades, the experimental explorations for creating new supramolecular materials and controlling their properties and assembly motif have risen exponentially. One of the most promising building blocks for this research has been porphyrins as mer units. These are tetrapyrrolic macrocyclic biomolecules, with important roles in natural biological processes, and thus they are occasionally utilised to mimic such processes in synthetic materials or nanodevices. [11] Moreover, their extended organic conjugated core and their ready availability for peripheral substitution and metallization makes them ideal candidates for versatile combinations of noncovalent linkages in self-assemblies, creating numerous architectures in all three dimensions. The collection of controlling factors that produce the desired material have been explored in multiple projects, entailing structural control, [9] thermodynamic or solvent conditions. [12,13] Thus far, these investigations have yielded a multitude of materials with notable photoelectronic device applications [14] and very impactful use in medicinal materials such as in photodynamic therapy. [15] One of the major bottlenecks hindering the advancement of this field has been the costly and time-consuming nature of synthetic explorations of these materials. Many synthetic processes can be lengthy. [16] Therefore, testing the nature and effects of chemical tuning or assembly dynamics is currently a slow process with little knowledge of any global categorisation of these effects for synthetic recommendations. Although encouraging, only a very few databases have been published linking either crystallographic [17] or synthetic data to assembly motif. [18] However, the initial mer structural design for new supramolecular materials remains intuitive and explorative. Recently, theoretical chemistry approximations were able to offer reasonably accurate and tractable predictions for supramolecular systems through improved semiempirical methods, thus offering fertile ground for systematic synthetic recommendations that may allow this field of material discovery to finally bloom. [19] It is a common approach to validate supramolecular aggregation data with quantum chemical semiempirical methods. [10,20] Very few studies have offered a predominantly in silico study of aggregates thus far. [21] Yet recent advancements and improvements, such as the 2019 development of GFN2-xTB [22] and the D4 dispersion correction, [23] have offered improved reliability for the modelling of large noncovalently bound species and an increasing number of more theoretical focused studies have been seen. [18] Herein, we report a theoretical investigation of exclusively porphyrinic assemblies in the 1D and 3D design. Inspired by published experimental work on porphyrinic self-assemblies, the role of steric and electronic effects of peripheral substitution, and the central metal held in the porphyrin is explored for controlling the assembly motif. This theoretical path highlights the strengths as well as obstacles currently occurring in computational predictions and aims to produce some quantifiable synthetic recommendations for a successful structure-property design. The aim of this study was to investigate in silico the chemical tuning factors that promote different aggregation motifs in supramolecular porphyrins. This was achieved by incrementally varying the nature of peripheral substituents in experimentally studied aggregates in terms of steric bulk or electronic effect, as well as the central metal when present. The study undertook 1D wires, which assemble in H (face to face) or J (head to tail) conformations, initially. Subsequently, the preference for 1D or 3D assembly was investigated with a second set of aggregates in a systematic set of computations. Through this path of consecutive semiempirical computations, the notable advantages as well as obstacles of the currently developed Hamiltonians are highlighted in order to open the dialogue for further advancements in the field.

The Measurement of Intermolecular Slippage
The threshold between a J and an H aggregation is described in the exciton coupling theory of Kasha. [24] The strength of head to tail (in J) interactions produce a varied spectroscopic effect as compared to a face-to-face (in H). The threshold angle that is formed by those planes is suggested at a value of 55.5 o . [25] To deduce this angle from the optimised computed structures, the monomer units were set as planes (each plane considered the planar porphyrin core). The axis of the stack is defined by two central (metal or dummy) atoms. The angle of the projection of the axis on the below porphyrin plane is the angle defining a J or H mode (Figure 2). It is shown that defines the degree of slippage between adjacent mer units in a supramolecular polymer stack.
The angle according to Figure 2 is then defined in the optimised coordinates as the angle formed between two successive central metal/dummy atoms and the succeeding connected nitrogen forming the narrowest angle of the porphyrin plane, given as numbered (1-2-3) = ο . The of the stack is given as the average of all consecutively measured angles.

Computational Methods
The weak nature of NCIs dominant in supramolecular aggregates demands a good treatment of the dispersion effects present for an accurate prediction of their geometry and Binding Energy (BE). Traditional computational mean field approaches have been notably poor in this description. [26] This is associated with the nature of the dispersive forces and the correlation of instantaneous fluctuation in the electron clouds of the adjacent interacting regions. This has been initially presented as a computationally demanding description especially for molecules the size of porphyrins. Several approaches to tackle this issue have been presented with empirical corrections, the most successful approach to this challenge has been a dispersion correction term Scheme 1. The chemical structures of the zinc porphyrin mer unit whole as in ref. [12] above and fragmented with a range of additional sterically ranged substituents 1A to 1H below.
to the calculated mean field energy. For all methods Grimme's D3 [27] dispersion correction was used, except for GFN2-xTB (using the updated D4) with the Becke-Johnson's dampening function, which is applied to avoid incorrect description below what is considered a common van der Waals distance. The size and complexity of the aggregates concerning this work suggests semiempirical computations as optimal for these predictions. A further consideration is that solvent and thermochemistry calculations have been considered to be nontractable for this set of coordinates. Therefore, all computations were carried out in vacuum with no constraints. Several published benchmarks and previous supramolecular computations in the literature suggest PM6-D3H4X, [28] PM7, [29] and GFN-xTB [30] and GFN2-xTB for optimal cost-accuracy balance for this load of calculations. These literature suggestions have been further investigated and supported in a benchmark in Tables S1a,b and S2a,b in the Supporting Information.
Geometries of structures 1S, 1M, and 1L were drawn in Avogadro 1.02 based on the Ogi et al. [12] structures and average known bond and intermolecular distances and angles, and the remaining structures of the first section (Scheme 1) are fragmented (via Avogadro 1.02) coordinates of those. All structures were optimised in the gas phase with no constraints with GFN2-xTB, except aggregates 1Z and 1H that were also optimised with PM6-D3H4X and PM7 for comparison (Table 3). Though frequency calculations were not tractable for this set of aggregates to assure a minimum in the potential energy surface, it was confirmed that all optimisation criteria were reached for all aggregates. The BE of the aggregates were calculated as Where n is the number of mer units in the aggregate, E n is the ground state energy of an aggregate with n mer units and E 1 the energy of the monomer. These are reported in kcal mol −1 .
All GFN2-xTB calculations were carried out via the xTB program. [22] The PM6-D3H4X and PM7 calculations were carried out with MOPAC2016. [31] The default convergence criteria were used in all cases. A conformational search analysis was used with CREST [32] via the xTB program for structure 1Z. All molecular illustrations were produced with Avogadro 1.02 [33] and all graphs were plotted with MS Excel.

Porphyrin Aggregations in 1D, J versus H
This section aims to evaluate the structural mer unit factors affecting the 1D aggregation motif of supramolecular porphyrins. These aggregation motifs have been categorised in two main forms the H and J (Figure 3). The findings of Ogi et al. [12] were utilized as a representative study for J versus H porphyrin aggregation, which reports thermodynamic control of the intermolecular slippage angle in three porphyrin stacks 1S, 1M, and 1L (Scheme 1).
These aggregates differ only on the two lateral substituents which are very small for 1S (a proton), 1M has a larger peripheral group (4-methoxy phenyl) and 1L, even larger, the chiral 4-[(R)-(+)-2-phenyl-1-propoxy]phenyl groups, implying an increment of the steric effect on the aggregation motif. The distinction between H and J conformers is described in the Experimental Section and is defined according to the angle . The investigation commenced by initially replicating the whole original structures as in 12 ] with GFN2-xTB. The predicted geometries were in good agreement with the published data, showing a pronounced that ranges around the threshold of distinction ≈57 o (Table S4, Supporting Information). For computational economy, the structures were then fragmented by strategically removing the long side chains, thus producing structures 1A, 1Z, and 1H (Scheme 1) as detailed in Table 1.  The original aggregates are assigned in Ogi [12] as H for 1S and 1M and J for 1L under thermodynamic control. The semiempirical computations herein agree in distinction for 1S and 1M being H (57 o ), yet an unexpected higher is produced for 1L. One of the initial considerations is the value of = 58 o being very close to the threshold of distinction (55.5 o ). The method of measurement could compensate for this discrepancy, yet it does not reflect a low enough value to account for the spectroscopic effects seen for a true J aggregate. This raises a question on the reliability of the predictive method. Consequently, when the peripheral groups are altered for the remaining aggregates 1A (hydrogen), 1B (methyl), 1 (isopropyl), and 1 (tert-butyl) the issue within the methodology was revealed. The aggregates have incrementally increased peripheral steric bulk and showed a clear decrease in the angle as the steric bulk increased ( Table 2).
The obvious difference between the latter set and the 1L (1H) misprediction is the existence of a benzene ring in the side substituents. To examine if the issue is present in the other relevant semiempirical Hamiltonians, was re-evaluated with PM6-D3H4X and PM7 for tetramers 1Z and 1H (that contain the benzene part in the side groups) shown in Table 3. The repredicted values only further confirmed this discrepancy, showing even higher (more H) values than the xTB method ≈59 o , further contrasting experimental findings. It should be further considered that ab initio or DFT methods that have  been reported as more precise for NCIs would not be efficiently applied for this type of study. The reason being that such a study would only be tractable for monomers or fragments of dimers of these entities, and as we have shown (Table S4, Supporting Information) is inconclusive for dimers of these aggregations. Consequently, the hypothesis for the misprediction in benzene-containing substituted porphyrins, may initially be rationalised by the rigidity present in the computational optimisation formed by these aromatic rings interacting and between adjacent units (Figure 4).
Computational optimisations are designed to locate the lowest energy structure in a potential energy surface. However, a real system is dynamic in vacuum or solvent with continuous rotations and vibrations. Therefore, to examine this hypothesis, a conformational analysis was carried out for 1Z (or 1M fragmented). The aim being to examine the number of H-type and Jtype conformers present in the ensemble of structures produced.
The range of angles for the 102 conformers from the CREST conformational analysis are shown in Figure 5. The majority of the conformers (72.6%) within the ensemble being a J-type aggregation, which is in contrast to the H-type isomer predicted based on the single lowest energy conformer. It is additionally considered that this set provides the ensemble of conformers for the default energy window of 6 kcal mol −1 . This may reflect too large  a window for the conformers corresponding to the experimental conditions. This conclusion can be reached by noting a majority of J conformer within this energy window, while an H conformer is experimentally recorded.
Given this result, an analysis of the ensemble set of the CREST conformers for a 3 kcal mol −1 energy window was performed. This revealed a contrasting higher H-conformer population (66.7%), in agreement with the experimental findings. [12] This adjustment of the energy window shows the importance of the choice of the ensemble of structures to correctly reflect the experimental data for porphyrin systems. This can be attributed to the systems used to develop the CREST algorithm which are mainly small organic molecules, aromatic fragments, and some metalized coordinated species. These do not reflect the aryl substituted porphyrins that are the focus of the work. However, in Figure 5 the angle for the 102 conformers of the CREST evaluation, where the yellow line indicated the threshold between a J (below the line) and H in (a) and (b), within the default 6 kcal mol −1 energy window of the CREST analysis, and the equivalent (c) and (d) for the adjusted 3 kcal mol −1 energy window for 16 structures.
The results of this analysis indicate that CREST calculations have the potential to correctly predict the J versus H aggregation preference for aryl substituted porphyrin systems, depending on the selection of an appropriate energy window. This is an impor-tant finding given that the default set GFN2-xTB did not reach the precise J versus H preference in these systems. Finally, it is worth noting that this work only studied a single system, thus additional systems need to be investigated to identify the optimal energy window that may be required for the correct prediction of J and H structures in porphyrin systems.
Beyond the computational predictive precision, this examination for 1D aggregates clearly confirms the trends that were originally observed by Ogi et al. In this work two lateral peripheral substituents of the mer unit were placed in increasing steric bulk 1S, 1M, and 1L respectively. The same principle was obeyed for the optimized structures presented herein 1A, 1B, 1 , and 1 yielding an increasingly pronounced (towards J), effectively agreeing with the initial hypothesis. The encountered misprediction of aryl substituted porphyrins has been presented for the first time here and it raises further examination of semiempirical predictions for this class of porphyrinic entities.

Examination of Factors Promoting 3D Assembly
The following data aim to evaluate the intermolecular interactions promoting 3D aggregation in porphyrin self-assemblies. This is achieved by a systematic examination of the effect of the central metal and the electronic effect of the peripheral groups on a model porphyrinic system, and by comparing the energetic preference for a 1D-wire or 3D-square motif. The dominant recorded strategies for synthesizing 3D supramolecular aggregates are the inclusion of a metal in the mer unit structure that coordinates with an electronegative group of an adjacent mer unit. Other strategies such as columnar assemblies on a surface have been successful. However, in this study we attempt to examine and replicate the attractive interactions that successfully produced the 3D-square assembly (Figure 6) from the study carried out by Oliveras-Gonzalez et al. [20] The decomposition of the forces promoting this cooperative assembly show a favourable Zn-pyridine coordination, that leads to a T-shaped stack of adjacent porphyrins and consequently promotes assembly in the third dimension. The original work additionally exploits a peripheral group size and hydrogen bonding for fine changes in spatial motif. These are important for illustrating the wide versatility of substituted porphyrins for forming different interactions in different topologies and of variable bonding strength. However, for a systematic investigation into the computational predictivity the focus is kept on the metal to linker interaction in order to show which is most favourable for promoting square (3D) or wire (1D) assembly.
For a systematic study the central metal held in the porphyrin was varied down and across the periodic table of elements, by selecting only appropriately divalent metals for keeping a neutral species, and thus omitting the additional complexity of charged interactions (these vary in computational predictivity). Furthermore, the peripheral groups were varied symmetrically on the porphyrin mer unit to differently hybridized groups. These are the cyano sp hybridized, and the carboxy sp 2 hybridized groups. The sp hybridized pyridine group, presented in the original work of Oliveras-Gonzalez [20] was also included for comparison, yet presented previously encountered issues as detailed in the next section. The sum of these combinations produced the aggregates listed in Figure 7.
The 3D aggregation motif is driven by different forces than a 1D stack. In a 1D-wire aggregate the -interactions between two consecutive conjugated stacked porphyrins, the metalmetal interactions and any -, vdW or hydrogen bonding (when present) of the peripheral substituents are the NCIs driving the spatial directionality. In contrast, the 3D square motif is dominated by metal-ligand coordination (Figure 8). The hypothesis herein is that differently hybridized electron withdrawing groups interact differently with various central metals, therefore promoting alternative aggregation motifs.   The BEs and aggregation motif of the optimized tetramers of this set were tabulated and compared in Tables 4 and 5. A first assessment, that simplifies the remaining analysis is seen for computations of all C-M aggregates (pyridine substituted). Interestingly, the computational data for C-Zn predict a wire aggregation motif in contrast to experimental findings. A reflection on these results leads back to findings of the previous section, where all ring (benzene or pyridine) substituted porphyrins promote an H aggregation with GFN2-xTB computations. As this prediction is repeated here for a pyridine substituted porphyrin, this further highlights the difficulty encountered in computations of these systems where lateralinteractions are found.
With consideration to the above findings, where the pyridine remains an obstacle for a reliable predictivity of these systems, alternative substituents were chosen which do not form the prob- lematic -interactions. The sp hybridized lone pair of the cyano group and the sp 2 hybridized lone pair of the carbonyl bind differently with the adjacent coordinating metal. The first forms a synergistic interaction entailing electron donation from ligand to an empty orbital of the metal and electron donation from the filled d orbitals of the metal to the * antibonding orbital of the ligand (Figure 9). The sp 2 carbonyl ligand is only able to bond to the metal via -donation, making the bond weaker. This establishes the hypothesis that the stronger bonding of ligands capable of synergistic bonding to the metal (such as cyano) should promote the formation of the 3D-square assembly. In Table 5 the BEs for the 1D-wire and 3D-square assemblies of the tetramers are collected and their preference for a dimensionality motif is assigned. The cyano substituted aggregates show an even distribution between 1D-wire and 3D-square assembly with the cyano containing aggregates. However, there appears to be a preference for first row metals being predicted as 3D-square when substituted with cyano groups. Upon examination of the trend, when descending the group 10 metals there is no clear pattern, with Ni preferring consistently the square, Pd the wire and Pt a wire with the carbonyl group and a square with the cyano. In this case the strong metal-metal interactions (Figure 8) which promote the 1D-wire stack formation, might be at play.
Furthermore, on reviewing the carbonyl substituted porphyrins a clear trend is seen with all but the Ni-and Cd-porphyrin aggregates predicted as a 1D-wire. The collection of these findings can confirm a preference for sp 2 hybridised ligands for 1D-wire formation. Alternatively, sp ligands favour more 3Dsquare assemblies, particularly for the first-row metals. The sp hybridised ligand predictions possess, though, some complication as the metal-metal bonding is at play, therefore possibly distorting the hypothesized trends.

Conclusion
The first part of this work focused on the 1D-wire porphyrin stacks, that assemble in H or J motifs, with respective distinct properties. It was confirmed that the xTB method yields reliable predictions for geometries and binding modes of the original aggregates, with the optimized 1L, 1M, and 1S, being generally in good agreement with the literature description (as aggregation interactions are concerned) and showing cooperative endergonic BEs. It was found that the molecular structures can be reliably studied once fragmented, keeping consistency in predictions with the whole. A clear trend was revealed where an increase in steric bulk of the peripheral substituents of the porphyrin leads to J motif of aggregation.
The computational mispredictions found within this data set revealed, for the first time, a short-coming of the semiempirical methods, where aromatic -interactions of the peripheral substituents strongly promote an H motif. A CREST analysis provided an attempt to reflect the dynamic free aggregate. This computational analysis opens the dialogue for more creative investigations into this issue. A larger set of porphyrin assemblies, and more in-depth examination of the conformational energy window could provide a better insight or further improvement into the semiempirical approaches for supramolecular studies.
The key factors known to promote 3D assembly namely metalligand coordination were examined in the second section of this work. It was noted that the previously seen peripheral -interaction issue was observed anew, this time persistently yielding 1D-wire assemblies instead of the experimentally reported 3Dsquare, with the trend disappearing where nonaromatic peripheral groups were introduced. Overall, it was recorded that the metal to ligand interaction plays an important role in the motif preference. This section suggests that further parametrization of methods might be needed to yield possibly consistent trends through metalation and substitution that were not seen herein.
It is important that more in depth computational investigations on the subject are undertaken for promoting a further advancement in supramolecular material discovery and their computational predictivity. The demand for such materials has risen exponentially yet their discovery has still to catch up to this necessity. We aspire that this work opens the dialogue for further methodology assessments and new improvements in chemical predictivity for supramolecular aggregates.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.