Origins of Offset-Stacking in Porous Frameworks

: Parallel-displaced π -stacking in the benzene dimer and larger polycyclic aromatic hydrocarbons is driven by competition between dispersion and exchange-repulsion interactions. The present work examines whether the same is true in porous frameworks that exhibit stacking interactions, including the [18]annulene dimer, porphyrin dimer, and several models of the covalent organic framework known as COF-1. Interaction energies and their components are computed using extended symmetry-adapted perturbation theory along two-dimensional scans representing slip-stacking. As in the polycyclic aromatic hydrocarbons studied previously, we find that the van der Waals interaction potential (defined as the sum of dispersion and Pauli repulsion) drives the system into a slip-stacked geometry. Electrostatics is a relatively small component of the total interaction energy. In the case of COF-1, the van der Waals potential drives the conformational preference whether or not a solvent molecule intercalates into the framework, although the presence of the guest (mesitylene) molecule substantially limits the low-energy slip-stacking configurations that are available. Even when the COF-1 pore is empty, a modest lateral offset of ≲ 1.5 Å is preferred, which is small compared to the pore size.


INTRODUCTION
The benzene dimer has long been studied by quantum chemists as an example of a π−π stacking interaction, 1−11 although its archetypal status in that regard has been questioned. 10Complexes of larger polycyclic aromatic hydrocarbons, starting with naphthalene dimer, demonstrate a much stronger preference for cofacial stacking as compared to the "T-shaped" or C−H•••π geometry that is similar in energy in the case of (benzene) 2 but is much less stable for larger acene dimers. 10All of the acene dimers exhibit offset π-stacking, in which the center-to-center coordinate of the aromatic rings is parallel-displaced, e.g., by 1.8 Å in the case of benzene dimer.Although long attributed to a competition between dispersion interactions and quadrupolar electrostatics, the origin of slipstacked π−π geometries (in benzene dimer, larger acene dimers, and in benzene on the surface of graphene) has more recently been shown to arise from a competition between dispersion and Pauli repulsion, with electrostatics playing a negligible role in driving offset stacking. 9−12 (In larger acene dimers, electrostatics does play an important role in stabilizing cofacial geometries relative to T-shaped orientations. 10) The aim of the present work is to examine whether these same conclusions hold for macrocycles that assemble into porous frameworks in the solid state.For example, the [18]annulene molecule (Figure 1b) packs in such a way that the π system of one monomer lies atop the pore of the neighboring monomer, in a form of offset-stacking. 13,14−18,23−28 The preference for eclipsedcofacial ("AA") versus fully staggered ("AB") stacking may be controllable via functionalization. 29n the present work, we examine stacking interactions in homodimers of the four monomers shown in Figure 1, with the goal of elucidating the intermolecular forces that drive offsetstacking.−12 For small noncovalent dimers where high-level calculations are feasible, XSAPT +MBD is comparable in accuracy (to within ≲1 kcal/mol) with the best ab initio benchmarks, 31,33,34 yet for large systems its monomer-based nature makes it more affordable than density functional theory (DFT). 30,35,36Crucially, XSAPT +MBD comes equipped with a physically motivated energy decomposition analysis, 11,12 which allows us to investigate the forces responsible for offset-stacking.The XSAPT+MBD energy decomposition is faithful component-by-component with respect to high-level ab initio benchmarks, again at the level of ≲1 kcal/mol. 34,37

METHODS
Systems investigated here include the benzene dimer, the [18]annulene dimer that we will call (annulene) 2 , the porphyrin dimer, and a COF-1 dimer that we will call (COF1) 2 .Monomer units for these four systems are shown in Figure 1.−18 The COF-1 monomer unit was optimized at the TPSS+D3/ def2-SV(P) level of theory, 38−40 as the TPSS+D3 functional reproduces the 3.4 Å face-to-face separation of the paralleldisplaced benzene dimer that is determined in high-accuracy calculations. 5The eclipsed-cofacial (COF1) 2 structure was then optimized, subject to a constraint that both monomers remain planar, resulting in a face-to-face separation of R = 3.28 Å that is used for subsequent two-dimensional potential energy scans using XSAPT+MBD.For (COF1) 2 (mesitylene), the same procedures results in an interlayer separation R = 3.4 Å that is used for the two-dimensional scans in that system.For (porphyrin) 2 and (annulene) 2 we use separations of R = 3.7 Å (TPSS+D3/def2-TZVP level) and R = 3.2 Å (TPSS+D3/def2-TZVPD level), respectively.All calculations were performed using Q-Chem (v.5.4), 41 and all structures are available in the Supporting Information.
Interaction energies were computed using XSAPT +MBD. 31,32Recent overviews and illustrative examples of this methodology can be found elsewhere. 11,12Briefly, XSAPT +MBD is a hybrid model in which the dispersion and exchange-dispersion terms of conventional second-order SAPT are replaced with a modified version 31,32 of the MBD model developed by Tkatchenko and co-workers. 42A chargeembedding scheme is also employed, 32 based on Charge Model 5 for the atomic charges, 43 although induction energies are small for the nonpolar monomers that are considered here, and the charge embedding changes the interaction energies only by ∼0.5 kcal/mol (see Table S1).All together, the result is a decomposition of the total intermolecular interaction energy (E int ) into physically meaningful components: 12 As in conventional SAPT, 44 these energy components are defined directly and not by difference (in contrast to supramolecular calculations of E int ), thus there is no basis-set superposition error and no need for counterpoise correction.
Components include electrostatics (E elst , defined as Coulomb interactions between isolated-monomer charge densities), exchange or Pauli repulsion (E exch , arising from the antisymmetry requirement and responsible for steric repulsion), induction (E ind , which contains both polarization and charge-transfer interactions), 45 and dispersion (E disp , for which we use MBD).In contrast to other recent XSAPT+MBD calculations, 9,10,33,34,46,47 we do not add the so-called "δHF" correction to E ind . 44This term accounts for induction effects beyond second order in perturbation theory but requires a supramolecular Hartree−Fock (HF) calculation that is expensive for the larger systems considered here, such as (COF1) 2 .Furthermore, we will demonstrate that induction energies are generally small in the systems considered here, and do not have any qualitative effect on the conformational energy landscape.
Previous work has demonstrated that XSAPT+MBD calculations with triple-ζ Karlsruhe basis sets afford convergence errors of ∼0.4 kcal/mol for small dimers, although convergence errors are larger when using def2-SVPD (∼1.4 kcal/mol) or def2-SVP (∼2.7 kcal/mol). 34For the larger coronene dimer, (C 24 H 12 ) 2 , triple-ζ basis sets again afford subkcal/mol convergence errors, although results in double-ζ basis The Journal of Physical Chemistry C sets were ≈10 kcal/mol from the basis-set limit. 34The doubleζ interaction energies are therefore not expected to be quantitative but we will see that a very similar qualitative picture arises regarding the origins of offset-stacking.For some of the COF-1 structures, we will report additional single-point calculations using the minimally augmented def2-ma-SVP basis set. 34

RESULTS AND DISCUSSION
We will first consider some one-dimensional interaction potential surfaces along a cofacial sliding coordinate that keeps the intermolecular separation fixed, corresponding to offset-stacking at fixed separation.
3.1.One-Dimensional Surfaces.3.1.1.Benzene Dimer. Figure 2 shows potential surfaces for (benzene) 2 as one monomer slides along the other at a fixed intermolecular separation of 3.4 Å, corresponding to the parallel-displaced minimum-energy structure. 9Previous work 9−11 suggests that the offset-stacking phenomenon in this system (as well as larger acene dimers) can be understood as a competition between Pauli repulsion (E exch ) and dispersion (E disp ).The latter favors placing atoms in close proximity, due to its rapid R −6 decay with atom−atom distance, but exchange-repulsion introduces a significant penalty for eclipsed-cofacial stacking because this aligns the density maxima that occur at the positions of the nuclei.The competition between these two forces is encapsulated in a "van der Waals" (vdW) energy, so named because vdW forces (as codified in the eponymous equation of state) involve both long-range attraction (E disp ) and short-range repulsion (E exch ).It can be seen from Figure 2b that these are the only two energy components that are significant in (benzene) 2 , as both E ind and E elst are relatively small and also relatively flat along the cofacial sliding coordinate that we consider.(The electrostatic energy is as large as |E elst | = 2.5 kcal/mol in the eclipsed-cofacial geometry, but the induction energy is truly negligible, |E ind | ≤ 0.2 kcal/ mol.)As documented previously using several different energy decomposition analyses (EDAs), 9 including the XSAPT+MBD method that is used here, the vdW surface E vdW exhibits the same conformational preferences as E int .In particular, the eclipsed-cofacial geometry is a saddle point that connects to a slip-stacked local minimum. 5The molecular physics is that Pauli repulsion drives displacement of the nuclei on one monomer relative to those on the other.Displacement by 1.8 Å relieves enough of the repulsion associated with the nuclear cusps in the electron density, while preserving favorable dispersion interactions that are reduced as a function of parallel-displacement. 9n contrast to this physical picture, electrostatic interactions are still widely invoked to explain offset-stacking, 14 often under the moniker of the "Hunter−Sanders model". 51This interpretation has been questioned in numerous computational 9−11,52−59 and experimental 60−63 studies, however.−11 Although E elst does change as the intermolecular distance is decreased from 3.8 Å (corresponding to the eclipsed-cofacial saddle point) to 3.4 Å (parallel-displaced local minimum), 64 Figure 2b shows conclusively that the gradient of E elst points toward the eclipsed-cofacial geometry.Electrostatics may set the depth of the π−π interaction potential, 10 but it provides no driving force for offset-stacking and in fact exerts a mild pull toward the sandwich geometry, although this effect is small in the case of benzene dimer and E elst has little effect on the energy landscape.
Asymptotically, electrostatics does have a role to play as it ultimately becomes repulsive, although it does so beyond 4 Å in the sliding coordinate (Figure S1b).Beyond that, charge penetration effects cease to be important, and E elst assumes its leading-order multipolar form, which is quadrupolar repulsion. 12,65This compensates for the second local minimum in E vdW that can be observed beyond 4 Å in the sliding coordinate, which does not exist in E int , and its presence in E vdW originates in the fact that exchange-repulsion decays faster than dispersion.This suggests that E vdW provides a meaningful The Journal of Physical Chemistry C explanation for short-range conformational preferences of benzene dimer but should not be used for asymptotic analysis.

Porous π-Stacked
Dimers.With this analysis of (benzene) 2 serving as a backdrop and a baseline, we next consider the porous systems: (annulene) 2 , (porphyrin) 2 , and (COF1) 2 .Figure 3 presents one-dimension energy component profiles for the [18]annulene dimer, in a similar manner to what was presented for the benzene dimer, starting in an eclipsed-cofacial configuration and sliding one annulene monomer across the other at 3.2 Å face-to-face separation, corresponding to the minimum-energy separation in the offsetstacked structure.At this separation, a local minimum in E int is observed at a parallel displacement of 3.5 Å relative to the eclipsed-cofacial structure.
One-dimensional profiles of E int and E vdW in Figure 3a bear much qualitative resemblance to the corresponding quantities in benzene dimer, except that the magnitudes are larger because the overall π-stacking well depth is considerably larger for (annulene) 2 .Also similar to the case of benzene dimer is the fact that the induction energy remains negligible at ≤1 kcal/mol, or <5% of the total interaction energy.What is different in this case is that E elst is no longer negligible and varies from −22 to −1 kcal/mol (Figure 3b).Nevertheless, E vdW mimics the overall profile of E int quite well, and no combination of electrostatics with another energy component can do so; see Figure S2.As in (benzene) 2 , the slope of E elst along the sliding coordinate points toward eclipsed-cofacial stacking rather than offset-stacking, such that electrostatics provides a driving force that opposes offset-stacking, favoring instead the close-approach of the monomers that enhances charge penetration, which is the origin of short-range electrostatic attraction. 12he analogous energy profiles for porphyrin dimer are provided in Figure 4.Because of the larger separation between the porphyrin monomers (3.7 Å) as compared to that in [18]annulene dimer (3.2 Å), magnitudes of the various energy components are smaller for porphyrin dimer.The induction energy is negligible (|E ind | < 1 kcal/mol), as it was in other systems consider thus far, but despite the similar size of (annulene) 2 and (porphyrin) 2, the electrostatic energy is much  From these results, a consistent picture emerges for the porous π-stacked systems.The vdW energy profile closely resembles the total interaction energy profile, even more so for the porous π-stacked dimers than for benzene dimer.As compared to the latter system, the larger (porous) dimers can relieve Pauli repulsion with a small parallel-displacement that leaves most atoms in close proximity, so that not too much favorable dispersion is lost.In sharp contrast to the conventional Hunter−Sanders (quadrupolar electrostatic) picture, 51 electrostatics opposes (rather than drives) offsetstacking.
3.2.Two-Dimensional Energy Profiles.In order to gain more insight about the role of the vdW contribution to the interaction energy, we next examine two-dimensional potentials representing potential scans at fixed face-to-face separation.(These scans represent rigid translations of fixed monomer geometries, with fixed orientation as well.)Only the vdW interaction energy (E vdW ) is considered in detail due to its qualitative resemblance to E int that was documented above, although two-dimensional scans for some other energy components can be found in Figures S7−S10.
3.2.1.Benzene, Annulene, and Porphyrin Dimers.The two-dimensional potential energy surface for benzene dimer (Figure 6a) has a similar topography to that reported in previous theoretical studies. 5,9,10It is repulsive in the eclipsedcofacial configuration (which is the coordinate origin in Figure 6), but remains so only for displacements <0.25 Å. Symmetryequivalent minima (with well depths of −2.4 kcal/mol) occur at about 2.0 Å displacement along directions that parallel the C−H bonds.(This corresponds to directions of 30°and 90°in Figure 6a; note that the one-dimensional cuts that are plotted in Figure 2 do not correspond to this direction but rather to

The Journal of Physical Chemistry C
horizontal displacement along the 0°direction.)The maximum well depth is slightly smaller than that reported in previous XSAPT+MBD studies where the δHF correction was included. 9,11he vdW potential surface (Figure 6b) has similar characteristics, although the orientation of the potential wells is rotated by 30°so that these occur along the 0°and 60°d irections.However, there is no other energy component (or combination of two energy components) that generates potential minima that can be identified with those occurring on the E int surface; see Figure S7 for examples.
For (annulene) 2 , similar two-dimensional contour plots of E int and E vdW are provided in Figure 7.As compared to benzene dimer, the repulsive region extends to larger displacements, consistent with a larger offset: 2.5 Å in (annulene) 2 as compared to 1.8 Å in (benzene) 2 .A potential well of −13.1 kcal/mol appears at 3.2 Å displacement along the vertical direction, which is the displacement direction of the one-dimensional scans in Figure 3, and another local minimum of −14.4 kcal/mol can be found at a displacement of 2.6 Å in an oblique direction.In this case, the vdW surface predicts local minima along similar radial directions but with an additional offset of ≈0.5 Å.These local minima in E vdW are much shallower (at −3.5 kcal/mol) than those in E int surface.These differences are ascribable to the more significant magnitude of E elst in the case of (annulene) 2 .For example, E elst = −22 kcal/mol at the coordinate origin in Figure 7, and can also be understood as the difference between E int (Figure 7a) and E vdW (Figure 7b).However, a two-dimensional scan of E elst (Figure S8b) does not reveal any local minima.Instead, E elst simply becomes less attractive in all directions emanating from the eclipsed-cofacial coordinate origin, again demonstrating that electrostatics opposes rather than drives offsetstacking.The only other energy component to exhibit any local minima is E ind (see Figure S8), but these minima have well depths of less than 1 kcal/mol.

COF-1 Dimer
Model.Two-dimensional energy profiles for E int and E vdW are plotted in Figure 8 for the (COF1) 2 model.In order to include large displacements representative of AB stacking, these scans extend to 18 Å displacements and use a coarse spacing of 1.0 Å to generate contours.Close-up views of the region up to 3 Å displace-ments, generated using a denser mesh, can be found in Figures S13 and S14.
At first glance, these (COF1) 2 surfaces are quite different from the ones considered above, but in fact they share many of the same features.Most important is that a slip-stacked conformation is preferred, with three local minima that are evident in Figure 8a and perhaps more evident in Figure S14a.The deepest of these (E int = −115.5kcal/mol) occurs at displacements Δx = 1.25 Å and Δy = 0.75 Å.In Figure S14a, one can observe a low-energy region corresponding to slipstacking by ≈1.25 Å in any direction.This is the magnitude of the lateral offset that our calculations predict, relative to the eclipsed-cofacial or AA-stacking geometry, and it is similar to offsets of 1.4−2.8Å that have been computed for this and other COFs, 16−18,24−28 depending somewhat on the material and on the computational method.Lateral offsets of this magnitude are much smaller than what one would expect for AB-stacking, in which the vertices of one COF layer are centered within the voids of the adjacent layer.In comparison to the 14.8 Å pore size of the COF-1 monomer, such small displacements are not much different from eclipsed-cofacial stacking.
Two other local minima, with E int = −42.7 kcal/mol and E int = −42.4kcal/mol, are evident in Figure 8a at coordinates (10.0, 6.0 Å) and (0.0, 12.0 Å), respectively, corresponding to radial displacement along polar angles of 30°and 90°.The vdW potential (Figure 8b) picks up all of these features, with semiquantitative energetics albeit not as attractive as the E int surface.As in the case of the [18]annulene dimer, the difference is electrostatics, and a two-dimensional E elst surface is shown in Figure S10b, with a close-up view of small lateral displacements in Figure S13b.The E elst surface does display some structure that is not unlike E int (or E vdW ), although the potential wells are much shallower and E vdW is a much better match to E int .For (COF1) 2 , we thus conclude that dispersion provides the dominant stacking interaction, albeit significantly enhanced by electrostatics, and that the competition between dispersion and Pauli repulsion (i.e., the vdW interaction) is sufficient to explain offset-stacking.(See Figures S13f and S14b for close-up views of E vdW .) In contrast to this rather consistent picture, a previous investigation of stacking interactions in a semi-infinite two- The Journal of Physical Chemistry C layer model of COF-1 concluded that offset-stacking is driven by electrostatics and dispersion, with no role for Pauli repulsion. 18This analysis was based on results from a periodic EDA 66 applied to periodic PBE+D3/TZ2P calculations, which favor an offset of Δx = 1.8 Å and Δy = 3.3 Å at an interlayer separation of 3.1 Å. 18 These offsets are at least roughly consistent with (if somewhat larger than) those predicted in the present work and in other DFT calculations, 25 and we confirmed that PBE+D3/def2-ma-SVP calculations for the finite (COF1) 2 model that is considered here do predict a larger offset.(A two-dimensional E int surface at the counterpoise-corrected PBE+D3/def2-ma-SVP level of theory is shown in Figure S16a.)This indicates that the different offset as compared to the present work mostly reflects the level of theory rather than the use of a semi-infinite model in ref 18.
Whereas the overall slip-stacking behavior obtained at the PBE+D3 level is therefore similar to that obtained using XSAPT+MBD, several aspects of the analysis presented in ref 18 seem questionable.First, it is not clear how large the basisset superposition error (BSSE) might be, or where it appears in the various components of E int when the periodic EDA developed in ref 66 is applied.Although BSSE is often a minor concern at the DFT/triple-ζ level, it does grow with system size and we have documented examples where the BSSE in a DFT/def2-TZVP calculation can be as large as 7−8 kcal/mol for systems with ∼300 atoms. 67For a layered material, BSSE should be size-extensive and will therefore increase in larger models of COF-1.
As noted in ref 18 and confirmed above (using both XSAPT +MBD and PBE+D3), electrostatics alone does favor offsetstacking in COF-1, although the authors of ref 18 rationalize this using a common trope about repulsive nuclear−nuclear and electron−electron interactions that has been thoroughly debunked. 12The EDA used in ref 18 is essentially a periodic version 66 of the Kitaura−Morokuma scheme, 68 refined for use with DFT, 69 and it employs a quasi-classical electrostatic analysis based on isolated-monomer charge densities.This is consistent with the definition of E elst in XSAPT, 12 and this definition of electrostatics is often attractive at vdW contact distances (meaning E elst < 0), even in cases where leadingorder multipoles of the monomer charge densities suggest that electrostatics will be repulsive when the monomers are well separated. 12,58,65Indeed, the original authors of the periodic EDA scheme note explicitly that the electrostatics term is usually attractive. 66We investigated this by performing SAPT +MBD calculation (with no charge embedding) using the PBE functional to compute the monomer charge densities.Although this type of "SAPT(KS)" calculation is strongly discouraged when combined with density functionals (such as PBE) that have incorrect asymptotic behavior, 37 such a calculation does allow us to obtain quasi-classical electrostatics based on a PBE description of the monomers.Twodimensional scans of the energy components computed at the SAPT(PBE)+MBD level for (COF1) 2 are shown in Figure S15, which should be compared to XSAPT+MBD results in Figure S10b.The electrostatic surface is in semiquantitative agreement at both levels of theory, so indeed the PBE electrostatics predicts an offset.However, the well depth is significantly shallower than E int .The vdW surface from the same SAPT(PBE)+MBD calculation (Figure S10f) is a better overall match to E int , with E elst as the last bit of the interaction energy that enhances offset-stacking.The SAPT(PBE)+MBD calculations predict a lateral offset that is consistent with XSAPT+MBD, and smaller than that predicted at the PBE+D3 level.
A final peculiarity of the PBE+D3 calculations reported in ref 18 is the authors' decision to include the correlation energy within the dispersion energy.While it is true that the D3 correction alone is not a pure dispersion energy, 70−73 because it is an empirical correction that is contaminated by other contributions (in the vdW contact region where the monomer charge densities overlap), certainly not all correlation effects are dispersion.This definition has the curious effect that dispersion alone favors offset-stacking in the analysis of ref 18, The Journal of Physical Chemistry C by about 5 kcal/mol.Physically speaking, E disp usually favors geometries that maximize the number of atoms in close spatial proximity, due to its steep falloff with distance, therefore dispersion often competes directly with Pauli repulsion.In porous systems such as COF-1, the dispersion energy should therefore be most attractive in the eclipsed-cofacial geometry, because any sizable offset places the atoms of one layer atop voids in the adjacent layer(s).Indeed, we observe this behavior using the MBD model for E disp , whether it is used with our preferred LRC-ωPBE description of the monomers (Figure S10e) or in conjunction with PBE (Figure S15e).Even the D3 model is most attractive in the eclipsed-cofacial configuration (Figure S16b).As such, the conclusion in ref 18 that dispersion alone favors offset-stacking appears to result from a misattribution of correlation effects.
3.2.3.COF-1 with Intercalated Mesitylene.We next investigate the effect of an intercalated mesitylene molecule in (COF1) 2 , as it is reported that the structure of COF-1 is sensitive to the presence or absence of residual mesitylene from the solvent. 16Two-dimensional E int and E vdW surfaces for (COF1) 2 (mesitylene) are presented in Figure 9 and twodimensional scans for other energy components can be found in Figure S11.This system is not symmetric so it is necessary to consider both positive and negative displacements Δx and Δy, and due to the large number of displacements that are required to construct the potential surfaces, we used the SVP basis set whereas def2-SVPD was used for (COF1) 2 in Section 3.2.2.To validate and calibrate this choice, we recomputed the (COF1) 2 results with def2-SVP, with results shown in Figures S5 and S12.Although the global minimum value of E int is less attractive by about 20 kcal/mol in the smaller basis set (which is a matter of basis-set convergence), 34 there are no qualitative changes in the nature of the energy surfaces.
Examining the E int surface for (COF1) 2 (mesitylene) in Figure 9a, we note that the global minimum (with E int = −128 kcal/mol) resides not at the coordinate origin but rather at x = −1.25 Å and y = 0. Two other minima with nearly identical interaction energies can be found at x = 0 and y = ± 1.25 Å.
For displacements of ≳2 Å along either x or y, however, the total interaction quickly becomes repulsive as one of the COF1 monomers encounters the intercalant.This behavior is modeled rather well by the vdW potential (Figure 9b), which exhibits three minima in roughly the same locations as those in the E int surface.None of the other energy components show this type of behavior.The electrostatic interaction, in particular, exhibits a local maximum at the coordinate origin and becomes increasingly attractive upon lateral displacements, including displacements greater than 2 Å, while the exchange interaction is repulsive everywhere in this example.The dispersion energy is most attractive at the coordinate origin, in contrast to results in ref 18 that were discussed above.Unsurprisingly, steric interactions prevent further slip-stacking in this system.This analysis suggests that the structure is more likely to exhibit "AA" stacking (albeit possibly with an offset) when the pores are occupied by solvent, but that staggered ("AB") configurations might occur when the pores are empty, especially in an amorphous or polycrystalline material where individual microdomains may be kinetically trapped in geometries other than the free energy minimum, because the solvent-free material exhibits broad swaths of configuration space where the potential is strongly attractive.In any case, the vdW potential retains its interpretive power even in the presence of a mesitylene guest.
3.3.Additional COF-1 Models.Our calculations for (COF1) 2 agree with previous theoretical studies that predict offsets of 1.4−2.8Å that have been computed for this and other COFs, 16−18,24−28 including COF-1. 26Those studies were performed either using supramolecular DFT or sometimes with force fields, in the interest of using relatively large model systems.(The COF-1 monomer unit in Figure 1d already contains 102 atoms.)We therefore wish to validate our results for the (COF1) 2 dimer model using larger model systems, and to that end we next consider dimers of the three monomer units in Figure 10.These were constructed from the basic COF-1 monomer in Figure 1d, without further optimization.The "half-COF1" structure in Figure 10a represents half of the Geometries considered include the eclipsed-cofacial structure, which we will call (0, 0) because it is the coordinate origin in Figure 8, along with a few displaced structures that we label in terms of lateral displacement coordinates (Δx, Δy).These include Δx = 1 Å = Δy, which is close to the global minimum structure for (COF1) 2 according to the XSAPT +MBD calculations, and also Δx = 7 Å = Δy as a model of ABstacking, where both lateral displacements are about half the pore size.Some intermediate displacements are considered as well.Interaction energies for these displaced structures were computed using supramolecular DFT, with functionals that include PBE+D3, 39,74 PBE0+D4, 75 PBE0+MBD, 42 and ωB97M-V, 76 with results presented in Table 1 alongside XSAPT+MBD results.Due to the possibility of numerical linear dependencies in these larger models, all of these DFT calculations were performed using τ ints = 10 −16 au and τ SCF = 10 −6 hartree.
DFT calculations in Table 1 were performed using the def2ma-SVP basis set 34 and include counterpoise correction.This procedure affords negligible basis-set convergence errors in supramolecular DFT calculations, 67 and tests on (half-COF1) 2 in Table S2 demonstrate that counterpoise-corrected DFT/ def2-ma-SVP interaction energies are converged to within ≲1 kcal/mol of DFT/def2-TZVPD results.Basis set effects are somewhat larger (∼3 kcal/mol) at the XSAPT+MBD level (Table S3), although the energetic ordering of various structures is conserved upon enlarging the basis set from def2-ma-SVP to def2-ma-TZVP.Since total interaction energies exceed 40 kcal/mol even for (half-COF1) 2 , we regard this as semiquantitative.The largest discrepancy between the two basis sets occurs between the (0,0) and (1 Å, 1 Å) structures, where the exchange-repulsion energy is largest.This is consistent with the observation that E exch exhibits the slowest convergence of the XSAPT+MBD energy components. 34e hex-COF1 monomer has twice as many atoms as the original COF-1 monomer unit, the 3COF1 monomer has almost three times as many atoms, and interaction energies for the dimers scale accordingly.The eclipsed-cofacial or (0,0) structure is consistently higher in energy as compared to the slip-stacked structure with Δx = 1 Å = Δy, which is lowest in energy at each level of theory that is examined.This consistency across various model systems suggests that the original (COF1) 2 model affords a reasonable description of short-range interactions in the real material.
For the three geometries having 7 Å displacements in one or both lateral directions, we do obtain somewhat different results using different model systems.For the original (COF1) 2 model, these three structures exhibit similar energetics (with relative energies within ≈3 kcal/mol of one another), regardless of which density functional is used, but for larger models, there is some variation in which of these AB-type geometries is most stable.That said, in the larger model systems these 7 Å-offset structures are significantly less stable as compared to the geometry with 1 Å offsets, by as much as 30− 40 kcal/mol.Even without relaxing these larger structures, these observations strongly suggest that something close to a 1 Å lateral offset represents the global minimum structure for larger COF-1 models, at various levels of DFT that are known to be reasonably reliable for noncovalent interactions.
The picture that emerges, across several levels of theory, is one of a slip-stacked AA motif with a small (<2 Å) lateral offset.On the other hand, it is unclear why XSAPT+MBD significantly destabilizes the AB-type structures, relative to DFT results.Some part of this is likely a basis-set effect, as results for (half-COF1) 2 in Table S3 indicate that the (7 Å, 7 Å) structure is stabilized by 4.5 kcal/mol, relative to the (0,0) structure, as the basis set is enlarged from def2-ma-SVP to def2-ma-TZVP.An investigation of these discrepancies is ongoing.

CONCLUSIONS
A survey of "porous" π-stacked systems including (porphyrin) 2 , (annulene) 2 and (COF1) 2 reveals that a slip-stacked geometry is preferred in all cases, to varying degrees.Among various energy components and combinations thereof, the vdW interaction potential (which encodes the competition between dispersion and Pauli repulsion) proves to be the best qualitative match for the total interaction potential.This provides a physical basis to understand the emergence of offset-or slip-stacking in these systems, just as it does in benzene dimer 9 and other π-stacked aromatic hydrocarbons. 10,11Other energy components are simply too small in magnitude or do not generate a potential surface that matches the topography of the total interaction potential surface.As such, we find no compelling rationale for the continued invocation of "Hunter−Sanders" in the context of π−π interactions Examination of (COF1) 2 and related model systems proves to be especially interesting.For (COF1) 2 (mesitylene), where a solvent molecule occupies the pore, steric repulsion severely limits the available lateral offsets, strongly favoring eclipsedcofacial or AA-type stacking.Upon removing the guest molecule, a ∼1.5 Å lateral offset becomes the minimumenergy geometry.This displacement is still best explained by the vdW interaction potential although the role of electrostatics becomes important for quantitative purposes.AB-type

Figure 1 .
Figure 1.Monomers examined in the present work, with an indication of the pore size for each: (a) benzene, (b) [18]annulene, (c) porphyrin, and (d) COF-1.Gray atoms are carbon, white are hydrogen, blue are nitrogen, red are oxygen, and pink are boron.

Figure 2 .
Figure2.XSAPT+MBD interaction potential and energy components for (benzene) 2 , along a cofacial sliding coordinate that is indicated, at a fixed intermolecular separation of 3.4 Å that corresponds to the parallel-displaced minimum-energy structure.(a) Total interaction energies (E int ) and vdW interaction energies, E vdW = E exch + E disp .(b) E int and E vdW plotted alongside additional energy components.

Figure 3 .
Figure3.XSAPT+MBD interaction potential and energy components for[18]annulene dimer, along a sliding coordinate that is indicated, at a fixed intermolecular separation of 3.2 Å that corresponds to the parallel-displaced minimum-energy structure.(a) Total interaction energies (E int ) and vdW interaction energies, E vdW = E exch + E disp .(b) E int and E vdW plotted alongside additional energy components.

Figure 4 .
Figure4.XSAPT+MBD interaction potential and energy components for (porphyrin) 2 , along a sliding coordinate that is indicated, at a fixed intermolecular separation of 3.7 Å that corresponds to the parallel-displaced minimum-energy structure.(a) Total interaction energies (E int ) and vdW interaction energies, E vdW = E exch + E disp .(b) E int and E vdW plotted alongside additional energy components.

Figure 5 .
Figure 5. XSAPT+MBD interaction potential and energy components for (COF1) 2 , along a sliding coordinate that is indicated, at a fixed intermolecular separation of 3.28 Å that corresponds to the parallel-displaced minimum-energy structure.(a) Total interaction energies (E int ) and vdW interaction energies, E vdW = E exch + E disp .(b) E int and E vdW plotted alongside additional energy components.

Figure 6 .
Figure 6.Two-dimensional XSAPT+MBD energy component surfaces for cofacial sliding of benzene dimer at a fixed face-to-face separation of 3.4 Å: (a) total interaction energy and (b) vdW energy, E vdW .The two panels use different energy scales.

Figure 7 .
Figure 7. Two-dimensional XSAPT+MBD energy component surfaces for cofacial sliding of [18]annulene dimer at a fixed face-to-face separation of 3.16 Å: (a) total interaction energy and (b) vdW energy, E vdW .The two panels use different energy scales.

Figure 8 .
Figure 8. Two-dimensional XSAPT+MBD energy component surfaces for cofacial sliding of (COF1) 2 at a fixed face-to-face separation of 3.28 Å: (a) total interaction energy (E int ) and (b) vdW energy (E vdW ).The two panels use different energy scales.See Figure S14 for a close-up view of the region up to 3 Å lateral displacement.

Figure
Figure Two-dimensional XSAPT+MBD energy component surfaces for cofacial sliding of (COF1) 2 (mesitylene) at a fixed face-to-face separation of 3.4 Å: (a) total interaction energy and (b) vdW energy, E vdW .The two panels use different energy scales.For these scans, the location of the mesitylene molecule is fixed while one COF1 monomer is displaced laterally.

Figure 10 .
Figure 10.Substructures and superstructures of the COF-1 monomer unit from Figure 1d: (a) "half-COF1", which is half of the monomer unit, (b) "hex-COF1", an extended unit containing twice as many atoms as the COF-1 monomer unit; and (c) "3COF1", a tricyclic version of the monomer unit.