Effect of Unwanted Guest Molecules on the Stacking Configuration of Covalent Organic Frameworks: A Periodic Energy Decomposition Analysis

Effect on the Stacking Elucidating the precise stacking configuration of a covalent organic framework, COF, is critical to fully understand their various applications. Unfortunately, most COFs form powder crystals whose atomic characterisations are possible only through powder X-ray diffraction (PXRD) analysis. However, this analysis has to be coupled with computational simulations, wherein computed PXRD patterns for different stacking configurations are compared with experimental patterns to predict the precise stacking configuration. This task is often computationally challenging firstly because, computation of these systems mostly rely on the use of semi-empirical methods that need to be adequately parametrised for the system being studied and secondly because some of these compounds possess guest molecules, which are not often taken into account during computation. COF-1 is an extreme case in which the presence of the guest molecule plays a critical role in predicting the precise stacking configuration. Using this as a case study, we mapped out a full PES for the stacking configuration in the guest free and guest containing system using the GFN-xTB semi-empirical method followed by a periodic energy decomposition analysis using first principle DFT. Our results showed that the presence of the guest molecule leads to multiple low energy stacking configurations with significantly different lateral offsets. Also, the semi-empirical method does not precisely predict DFT low energy configurations, however, it accurately accounts for dispersion. Finally, our quantum-mechanical analysis demonstrates that electrostatic-dispersion model suggested Hunter and Sanders accurately describe the stacking in 2D COFs Abstract Elucidating the precise stacking conﬁguration of a covalent organic framework, COF, is critical to fully understand their various applications. Un-fortunately, most COFs form powder crystals whose atomic characterisations are possible only through powder X-ray diﬀraction (PXRD) analysis. However, this analysis has to be coupled with computational simulations, wherein computed PXRD patterns for diﬀerent stacking conﬁgurations are compared with experimental patterns to predict the precise stacking conﬁguration. This task is often computationally challenging ﬁrstly because, computation of these systems mostly rely on the use of semi-empirical methods that need to be adequately parametrised for the system being studied and secondly because some of these compounds possess guest molecules, which are not often taken into account during computation. COF-1 is an extreme case in which the presence of the guest molecule plays a critical role in predicting the precise stacking conﬁguration. Using this as a case study, we mapped out a full PES for the stacking conﬁguration in the guest free and guest containing system using the GFN-xTB semi-empirical method followed by a periodic energy decomposition analysis using ﬁrst principle DFT. Our results showed that the presence of the guest molecule leads to multiple low energy stacking conﬁgurations with signiﬁcantly diﬀerent lateral oﬀsets. Also, the semi-empirical method does not precisely predict DFT low energy conﬁgurations, however, it accurately accounts for dispersion. Finally,


Introduction
Over the years, a considerable amount of effort has been geared towards exploring the potential application of nanoporous materials in gas storage and sieving, filtration, extraction, separation, and catalysis. [1][2][3][4][5][6][7][8] The main breakthrough in this research dates from 1994, when Yaghi and co-workers introduced a new synthetic approach, commonly referred to as reticular chemistry, which uses carefully designed building blocks to create extended crystalline structures. [9][10][11] The peculiarity of this approach is the ability to assemble predefined crystalline structures from judiciously designed molecular building blocks whose structural integrity and rigidity are preserved throughout the assembly process. The successful implementation of this synthetic approach has resulted in an ever-increasing scientific and industrial interest in the synthesis and application of a plethora of nanoporous materials, which can be classed as Metal-Organic Frameworks (MOFs), [12] Metal-Organic Polyhedra (MOPs), [13] Zeolite Imidazolate Frameworks (ZIFs), [14] and Covalent Organic Frameworks (COFs). [15] COFs are an interesting class of nano-porous molecules, wherein organic linker molecules are stitched together by connector molecules, which are also organic, to form a framework that is held together by strong covalent bonds. In comparison to other porous materials, COFs present the advantage that they are made mainly from lightweight elements such as carbon, oxygen, hydrogen, nitrogen and boron.
Consequently, they have relatively low mass densities with a comparably large surface area, tunable pore sizes and structures, easily functionalisable and versatile covalent combination of building units. [16,17] Despite the advantages that COFs present over other nano-porous materials, there is a significant challenge in elucidating their packing sequence, which is of paramount importance for their optoelectronic, catalytic and sorption properties. For instance, in two-dimension, the organic linkers are topologically stitched together to form 2D monolayers. These monolayers in turn stack together in a third dimension via non-covalent interactions to form layered structures. [18] However, these weak non-covalent interactions between layers often lead to the formation of layered structures with ill-defined polytypes, which are either polycrystalline or amorphous solids. Unfortunately, a detailed atomic characterisation of these solids is not possible since they are insensitive to most analytical techniques such as solid-state nuclear magnetic resonance (NMR) or single X-ray diffraction analysis. [19] Hence, the characterisation of COFs predominantly rely on powder X-ray diffraction (PXRD) analysis. Although essential dataset for structural analysis and simulations can be obtained from PXRD analysis, the absence of single crystals renders the identification of a suitable stacking structure elusive. [20] Consequently, PXRD analyses are mostly coupled with other analysis to get a glimpse of the geometry and stacking arrangement of any newly synthesised COF. For example, transmission electron microscopy, TEM, analysis is often used to determine the composition of the COF. [21] Meanwhile, computational calculations are used to predict their geometric parameters and packing arrangement. This is generally done by computing educatively-guessed periodically bound systems and comparing computationally determined PXRD patterns to those obtained experimentally.
Early attempts in the elucidation of the stacking arrangement of 2D COFs led to the identification of two stacking arrangement, the eclipse (AA) and the staggered (AB) arrangement. In the eclipse arrangement, layers lie directly on the top of the other consequently aligning the hexagonal pores of each layer to generate a 1D pore system. [22][23][24] However, the perfectly eclipsed structures are not energetically favourable due to repulsive interactions between like-atoms from adjacent layers. Therefore, off-centred or laterally displaced energy minima, also referred to as slip-stacked structures, are commonly assumed as the eclipse stacking arrangement. [25][26][27] However, the exact slip-stacked structures for any COFs are so far elusive since layer offset can occur in any direction. Consequently, higher symmetry structures based on eclipsed layer stacking are usually assumed as an average structure model. [18] On the other hand, the staggered stacking arrangement involves a layering in which, each layer is offset by half unit cell resulting in three connected vertices laying directly over the centre of the pore of the adjacent layer. [22] This stacking arrangement is similar to layering in graphite in which three connected vertices, carbon-atoms, lie directly over the centre of a sixmembered ring of the adjacent layer. So far, such stacking arrangements have been observed only in COF-1. [15] To this day, the unprecedented staggered stacking arrangement observed in COF-1 is highly contended amongst COF chemists. The crystal structure of COF-1 is known to possess a guest mesitylene at the centre of the pore in each periodically bound layer, which could be the reason for this AB stacking arrangement. This can be justified by the fact that the COF was reported to loss its crystallinity when heated to 200 • C to remove the guest molecule. [15] From a series of computational study performed by Heine and co-workers, [26] the authors showed that in the absence of the guest molecule, the inclined and serrated eclipse structures were the energy minima. In 2013, Ravikovitch and coworkers performed an indepth study on the irreversible phase transition resulting from the removal of the guest mesitylene molecule. [22] More recently, Jiang and co-workers investigated the host-guest interaction of the mesitylene guest molecule in the confinement of pore channels in COF-1 as well as the effect of the guest molecule on the variation in PXRD patterns. [20] 5 The ambiguity surrounding the stacking pattern of COF-1 depending on the presence (AB) or absence (AA) of the guest molecule -is an extreme case where the effect of an unwanted guest molecule is critical in correctly predicting stacking configurations. For this reason, COF-1 represents an excellent candidate for exploring the effect of unwanted guest molecules (which are present in most experimentally synthesised COFs) in simulated stacking structures. We believe that this can apply to all COFs and should be taken into account when computational simulations have to be used in predicting the stacking configuration of newly synthesised COFs.

Modelling the stacking geometries
The variation of the non-covalent intermolecular interactions in the stacking of COF-1 was evaluated by plotting out full potential energy surfaces for both the guest free and guest containing structures of COF-1, which are represented in For an optimal computational cost, full potential energy surfaces (PES) to effectively evaluate interlayer heights between monolayers were first computed for both structures in Fig. 1 using the Geometry, Frequency, Non-covalent, eXtended Tight Binding (GFN-xTB) method. [28] The GFN-xTB method is a computationally and numerically robust semi-empirical tight-binding method that has been designed for computing highly accurate molecular properties for elements across the entire periodic table, including some lanthanide and actinide. In comparison to other semi-empirical methods that employ minimal atomic orbital basis sets, which limit the accuracy, GFN-xTB uses atom centred minimal basis set, which are Slater functions that are augmented with a second s-function and a d-polarisation function to enable a proper description of hydrogen bonding and hyper-valent bonding arrangements respectively. Moreover, it is properly parametrised to inherently account for dispersion correction using the D3 method, [29] which is effective in accounting for non-covalent interactions between 2D COF layers.
To generate the PES, bilayers resulting from the optimised monolayers were used.
Atoms from one layer of the stacked system were rigidly held in place meanwhile the other layer was translated over the surface of the rigid layer in increments of 0.1 Å along the a and b axes. At every a and b offset, the interlayer distance, c-axis, was span from 2.8 Å to 4.00 Å in increments of 0.01 Å and the intermolecular potential energies were computed. The configurations with the lowest energy along the c-axis were projected onto the ab plane to enable visualisation of the 3-dimensional data in two dimensions. Also, each of these configurations alongside their optimal interlayer heights was saved for use in computing the periodic energy decomposition analysis (PEDA). The range of layer translation and spacing were probed wide enough to include all minimum energy structures. However, layer rotations were not considered while generating the PES for the simple reason that it could result in a loss in order, which would be computationally intractable to model.
A PEDA was then computed for each of the previously computed configurations that constitute the GFN-xTB PES. In setting up this computation, each monolayer was used as a molecular fragment from which their interlayer interactions energies at each configuration were computed using the optimal interlayer heights. [30] A summary of this method as applied to this study is found in the Electronic Supporting Information (ESI) while interested readers are recommended to consult the following references for more details. [30][31][32][33][34][35] The PEDA PES was computed using first-principle density functional theory, DFT. To be precise, we used Grimme3 dispersion corrected [29] generalised gradient corrected Perdew, Kieron and Burke exchange-correlation functional PBE-D3, [36] alongside the doubly po-8 larised triple-zeta, TZ2P, basis set. [37] The PBE-D3/TZ2P, functional has been shown to provide an accurate representation of intermolecular interactions as well as lattice parameters, unit cells and bonding patterns of COFs, [30,38,39] thus making it an ideal functional for this studies. All computations were perform using the Amsterdam Modelling Suits (AMS) package version ADF2019.305. [40] 3 Results and Discussion

GFN-xTB generated potential energy surfaces
The PES computed using the GFN-xTB method for the guest free and guest containing structures of COF-1 are presented in Fig. 2a   Meanwhile, the staggered configuration, with a lateral offset corresponding to (7.5, 4.3), is observed to be destabilised from the minima by approximately 0.20 eV.
On the other hand, the PES of the guest containing systems of COF-1 is highly corrugated with multiple energy minima at different translational offsets. These energy minima, grouped in regions A, B, C and D in Fig. 2b, correspond to 3 discrete sets of energy minima with different unique stacking configurations. These numerous stacking configurations could potentially impede experimental reproducibility as a result of synthesised molecules self-assembling into any of these low energy configurations. For instance, only the stack layers in region C has PXRD patterns that matched those from Yaghi's original synthesis. [15] Moreover, numerous synthesis have been reported for the synthesis of COF-1 in literature, wherein authors reported lattice parameters that matched those of the originally synthesised COF to within 80-90 % accuracy. [41][42][43] Notwithstanding, most of the PXRD pattern from these syntheses do not fully match that of the originally synthesised COF, consequently suggesting that these newly synthesised COFs might have adopted one of the low energy stacking configurations.

Periodic energy decomposition analysis
A periodic energy decomposition analysis was performed to fully understand the interactions occurring between the COF layers. In the PEDA scheme, the total interaction energy is decomposed into electrostatic, Pauli repulsion, orbital interaction and dispersion term as fully described in the ESI. The PES for the total interaction energies between COF layers for both guest free and guest containing systems are presented in Fig. 3. For both COF systems the total interaction energy is observed to be slightly varian from the semi-empirical GFN-xTB counterparts.
In the guest free system, Fig. 3a, the minima are observed to be significantly misaligned from the fully eclipsed structure with the global minima having displacement coordinates of (1.80, 3.30, 3.10) in the a, b and c directions as represented in Fig. 4a.  Overall, the presence of the guest molecule provide a plethora of stacking configurations with significantly varied orientations, which are all within 0.2 eV. This implies that these stacking patterns that ranges from the AB stacking to the various slip-stacked patterns are all susceptible patterns that can be adopted from actual experiment. Although, this system represents an extreme case wherein the effect of the guest molecule is critical, we still believe that the effect of solvent molecules and other potential guest molecules present in COF crystals should be given more considerations especially for predicting their stacking patterns. Also, the difference in the sampled minima energy configurations between the total interaction energy computed using the first principle DFT method and the semi-empirical GFN-xTB (and the commonly used DFTB method reported by Heine and co-workers, [26] ) further highlights the need for more robust computational methods for simulating COF structures. The underlying short coming of the semi-empirical method alongside the effect of the guest molecule can further be investigated by exploring the various contributions to the total interaction energy.
The electrostatic energy is defined as the interaction between static charge densities in the COF layers. In general, it includes attractive Coulomb interactions between the nuclei of one layer and the electrons in the adjacent layer as well as the nuclei-nuclei and electron-electron repulsive Coulomb interactions between adjacent layers. The PES for the electrostatic energy contribution to the interaction energy is presented in Fig. 5. In stacked COF layers, the electrostatic energy is highly destabilised in the perfectly eclipsed configuration as observed in Fig. 5. This is because nuclei-nuclei and electron-electron repulsion are maximum in this configuration since like-atoms 14 from adjacent layers lie directly adjacent to each other. For this reason, the perfectly eclipse configuration is seldom observed in stacked structures with the same repeating layers. In the guest free system, the AB stacked configuration with coordinate (7.50, 4.30, 2.75) is the global minima. In this structure, each boraxine ring of the top layer lie at the centre of the pore of the adjacent layer, consequently preventing the highly repulsive interactions from electropositive boron-boron atoms and electronegative oxygen-oxygen atoms of adjacent layers as shown in Fig. 6a. In the guest containing system, the global minima on the electrostatic PES is a slip-stacked structure with displacement coordinates of (2.50, 2.60, 3.01), which is closely related to the global minima of the total interaction energy. In this structure, the layers are stacked such that a boron atom in the 6-membered boraxine ring on one layer lay at the centre of a 6-membered phenyl ring of the adjacent layer, similar to the stacking arrangement in graphite. This configuration maximises the attractive electrostatic interactions through the interactions of electropositive boron atoms of one layer and electronegative oxygen atoms of the adjacent layer.
The orbital energy term also known as the polarisation energy accounts for chargetransfer between the occupied and empty orbitals of the of adjacent layers as well as the reorganisation of charge in these layers. The PES for the orbital energy contribution to the total interaction energy is presented in Fig. 7

16
The sum of the electrostatic and orbital energy represents the quantum-mechanical Coulomb interaction between polarised charge distributions for the two monomers, without any multipole approximation, and fully accounts for interpenetration of the monomer charge densities. [44] The PES for this polarised Coulomb potential, which is represented in Fig S4 in the ESI, fully agrees with the Hunter and Sanders model that suggests that π-π interactions between stacked layers should always favour the slip-stacked configuration over the perfectly eclipse configuration. [45] The Pauli term contains the exchange and repulsion energies, which are stabilising and destabilising, respectively. The exchange interaction arises due to the antisymmetric nature of the wave function that allows the exchange of electrons between monomers, meanwhile the repulsion interaction originates from the interactions between electrons of the same spin on adjacent layers. The PES for the Pauli energy contribution to the interaction energy is presented in Fig. 8. For both systems, the Pauli interaction energy is nearly flat along all phenyl co-facial systems including the perfectly eclipse system. The displacement coordinates for the global minima are (7.50, 0.10, 3.41) and (0.00, 0.00, 3.49) for the guest free and guest containing systems, whose structures are presented in Fig. 9. Carter and Herbert [44]. In this study, which was performed on benzene dimer, the authors suggested that the Pauli repulsion, rather than electrostatics, provides the driving force towards the slip-stacked arrangement. However, we acknowledge that this discrepancy could be a consequence resulting from the underlying differences in EDA schemes [46] or the inability of this results to be applied in π-π interaction of extended aromatic systems such as 2D COFs. [47]  The dispersion energy term is defined here as the difference in correlation energy between the super-molecule and the individual monomer fragments. It takes into account all the long-range interactions and its accuracy depends on the density functional used. In this study, dispersion correction from the PBE functional is improved using the Grimme-D3 semi-empirical scheme. [29] The PES for the dispersion energy contribution to the interaction energy is presented in Fig. 10. In the guest free system, all energy minima below 0.4 eV are located within 4.0 Å meanwhile, the dispersion energy is significantly destabilised beyond this region.
On the other hand, the presence of the guest molecule significantly stabilises the dispersion energy beyond this region consequently stabilising the AB configuration.
Finally combining the dispersion energy with polarised Coulomb (Electrostatic + Orbital) energy provide a more sensible representation of the total interaction energy in comparison to the Pauli + dispersion energy as shown in Fig. S5 and Fig. S6 in the ESI. This agrees with the classical Hunter and Sander models, which 20 posit that π-π stacking in layered systems is governed by a competition between electrostatic and dispersion. [45]

Conclusion
In this study, we explored the different metastable stacking configuration at all possible lateral offsets for the guest free and guest containing systems of COF-1.
We showed that the presence of the guest molecule leads to multiple low energy stacking configurations at significantly different lateral offsets, which are all experimentally feasible. Although COF-1 represents an extreme case where the effect of the guest molecule is critical in predicting the stacking configuration, we contend that these results can easily be applied to any COF or crystalline system.
For this reason, we recommend that newly synthesised COFs should be accurately purified and their physical properties should be explicitly discussed with computational collaborators (since most of these crystals could easily absorb gases from the atmosphere, which might affect their experimental PXRD patterns) in order to ease simulations and prediction of stacking patterns. In addition, we showed that care should be taken when using semi-empirical methods in simulating stacking configurations. Although this semi-empirical methods quite accurately account for dispersion, they still miss out on precisely capturing the first principle lateral offsets. Finally, we showed that π-π stacking in the 2D layered system is correctly governed by the simplistic electrostatic vs dispersion model posited by Hunter and Sander.

Conflicts of interest
The authors declare no conflict of interest. proposed by Marokuma , [2] and later on developed by Gernot/ Mortitz Von for organic systems [3] and Ziegler/Rauk for transition states. [4] In the (P)EDA scheme, the intrinsic bond energy for the interaction of two molecular fragments A and B that form a molecule AB is obtained by separating the bond formation process into several sub-steps. The bond dissociation energy ∆E bond (negative of the dissociation energy without zero-point vibrational corrections De) is composed of a preparatory step and an interaction step with dispersion correction included as given equation 1.
Where ∆E prep is the preparation energy, ∆E int is the interaction energy and ∆E disp is the dispersion energy, whose inclusion is the main difference between the PEDA and EDA. This dispersion energy is included via the semi-empirical correction scheme DFT-D3 put forward by Grimme and co-workers through the difference of fragment and super-molecule dispersion energies. [5] The preparation energy , ∆E prep , represents the energy involved in distorting the configuration of the interacting fragments from their isolated ground state configuration A GS , B GS to the reference configuration A, B that they adopt in the molecule AB. This energy is composed of both geometric distortion and electronic excitation from the relaxed to the distorted fragments as presented in equation 2.
The promoted fragments can be represented by Slater determinants Ψ A and Ψ B built from fragment orbitals at A and B, respectively, which are interacting to form the molecule AB. Hence the interaction energy, ∆E int , results from the energy difference between the energy of the supermolecule, E AB , and those from the prepared fragments (monomers), E A and E B . Consequently, the basic idea in any EDA scheme is to partition the interaction energy, ∆E int , into well-defined terms as given in equation 3.
Where ∆E elstat correspond to the electrostatic interaction energy, which represent the static charge densities of each fragments within the supermolecule. This term includes the attractive Coulomb interactions between nuclei of one fragment (e.g

Setting-up PEDA Job
In order to set-up the PEDA job on AMS, each monomer was used as a separate fragment from which the stacked bilayer was considered to be the super-     The displacement coordinates for 100 energy minima configuration for the total interaction energy are presented in Table S1. Coordinates and energies for the entire PES can freely be downloaded from http://doi.org/10.5281/zenodo.4058818 Continues on next page Continues on next page