Band gap opening from displacive instabilities in layered covalent-organic frameworks $^\dag$

Covalent organic frameworks (COFs) offer a high degree of chemical and structural flexibility. There is a large family of COFs built from 2D sheets that are stacked to form extended crystals. While it has been common to represent the stacking as eclipsed with one repeating layer ("AA"), there is growing evidence that a more diverse range of stacking sequences is accessible. Herein, we report a computational study of layer stacking in two prototypical COFs, Tp-Azo and DAAQ-TFP, which have shown high performance as Li-ion battery electrodes. We find a striking preference for slipped structures with horizontal offsets between layers ranging from 1.7 \r{A} to 3.5 \r{A} in a potential energy minimum that forms a low energy ring. The associated symmetry breaking results in a pronounced change in the underlying electronic structure. A band gap opening of 0.8 - 1.4 eV is found due to modifications of the underlying valence and conduction band dispersion as explained from changes in the $\pi$ orbital overlap. The implications for the screening and selection of COF for energy applications are discussed.

Various forms of disorder exist in experimentally-synthesized COFs, such as bond breakage, pore collapse and stacking faults. Such imperfections can significantly affect the properties of 2D COFs, causing loss of crystallinity, porosity, and conductivity. [33][34][35][36][37] In particular, the interlayer stacking modes of 2D aromatic COFs play a critical role in determining their properties. The stacking behaviour of COFs is not thoroughly understood, however, due to difficulties in experimental characterisation of the dynamic, low-crystallinity materials. For instance, powder X-ray diffraction (XRD) can only detect the existence of crystalline domains, making the extraction of accurate results difficult in the presence of low long-range order and sizeable thermal dynamics. 7,38 As such, XRD measurements struggle to quantitatively distinguish crystalline structures from other similar aggregated structures, as a result of peak broadening in the diffraction pattern. 6,36,39,40 To achieve greater resolution of COF layer stacking, Kang et. al 39 recently used 13 C solid-state nuclear magnetic resonance (ssNMR) to distinguish different aggregated structures by studying the interactions between atoms and chemical groups from adjacent layers.
Five different stacking modes in 2D COFs have been reported: eclipsed, inclined, zigzag, staggered and random stacked. 36,41 The eclipsed stacking (AA) corresponds to zero horizontal (coplanar) offset between neighbouring layers in the ab plane, which has the highest symmetry and is the most often reported in experimental works. The inclined stacking (AA') corresponds to a constant, collinear offset between neighbouring layers. This stacking mode was observed using powder XRD and transmission electron microscopy (TEM) in SIOC-COF-8 and SIOC-COF-9. 42 Zigzag stacking (AB) corresponds to an alternating offset direction between layers but it still retains high porosity in the stacking sequence. Staggered stacking is a special type of AB stacking, whereby the offset between layers is sufficient to make one layer's skeleton centered directly above the other pore, e.g. a horizontal offset halfway along the ab unit cell diagonal. This large offset between layers would reduce the porosity completely in the structures. 36 These four stacking modes can be combined to form a random stacking sequence 41 , which is difficult to characterize experimentally or computationally due to limitations in equipment precision and computational demand.
Several studies have focused on stacking modes and their effect on properties for various 2D COFs. [33][34][35][36][37] It has been found that the AA stacking mode is the most energetically unfavorable as a result of strong repulsive interlayer orbital interactions. 33,34,38 Koo et al. studied the potential energy surface (PES) of 33 COFs using molecular mechanics (MM) and density functional theory (DFT) approaches, finding that COFs are preferentially stacked with 1-2 Å horizontal offsets between layers. 34 It has been reported that bulk COF structures have either inclined or zigzag stacking, which are more energetically favorable than eclipsed and staggered stacking. 6,43 The simulated XRD patterns of unidirectionally slipped (AA') and alternating slipped (AB) modes show a better agreement with the experimental XRD pattern than the eclipsed structures. 36,43 More precisely, in many studies 40,42,44 , the predicted diffraction patterns of inclined stacking are more consistent with experiment than other stacking modes.
COFs of Tp-Azo 45 and DAAQ-TFP 46 , whose experimentally reported crystal structures are shown in Fig. 1(a) and (c), have been reported with high energy capacity, good cycling performance and excellent stability as battery electrodes. 22,[46][47][48] It has been predicted that 30 Li + ions per unit cell can be inserted into and extracted from the porous Tp-Azo structure, using DFT simulations. 47 DAAQ-TFP COF linked by β -ketoenamines 49,50 was the first COF to exhibit reversible redox behavior in energy storage systems and has the highest surface area of all COFs linked by either imines or enamines. 13,46 However, the basic structural prop-erties, the stacking modes and the electronic structures in Tp-Azo and DAAQ-TFP COFs have not been reported. In this work, we present a theoretical study of the bulk properties and potential energy surface for stacking fault disorder of these two COFs. Furthermore, we investigate the effect of the stacking sequence on the electronic structure, rationalising the behaviour through consideration of the interlayer orbital interactions, and discuss the implications for COF material design for energy applications.

Methods
All electronic structure calculations were performed using Kohn-Sham DFT through the all-electron "Fritz Haber Institute ab initio molecular simulations" FHIaims package. [51][52][53][54] Both the semilocal functional of Perdew-Burke-Ernzerhof revised for solids (PBEsol) 55 and the hybrid Heyd-Scuseria-Ernzerhof (HSE06) 56 functional were used, and the Tkatchenko-Scheffler correction was implemented to account for van der Waals (vdWs) interactions between layers. The PBEsol functional was used for geometry optimisation, having been shown to predict atomic structures and energies of solid materials with good accuracy. 55 The HSE06 functional was used for calculations of electronic band structures, having been shown to accurately reproduce the electronic structure across a range of semiconductors. ? A k-point grid of 1 × 1 × 10 was used for the geometry optimisation with 6 × 6 × 6 sampling used for electronic structure analysis. An energy convergence criterion of 0.01 meV per unit cell was used with an atomic force tolerance of 0.01 eV/Å.
The initial crystal structure parameters of Tp-Azo and DAAQ-TFP were obtained from the CoRE-COF database 25 . These structures were firstly relaxed using the lighter Tier 1 numerical basis set, followed by a relaxation with the expanded Tier 2 basis set, before calculating the energetic and electronic properties. The well-converged conventional "intermediate" basis functions for each element species were used in the band structure calculations.
The relaxed structures were modified to study the stacking fault behaviour. Layer displacement was modelled by changing the angles of α and β of the unit cell, thereby shifting the individual pseudo-hexagonal layers along the ab plane to yield inclined stacking modes. 33 Fig. 2 shows the slip grid of one layer to the another layer along the a and b sides, with offsets of -6 Å to 6 Å and displacement steps of 0.5 Å. The distance between adjacent layers was kept fixed to that of the relaxed structures. The energies of the displaced structures were then calculated with fixed atomic positions and with relaxed atomic positions. The relaxed displaced structures were used to study the effect of the stacking faults on the physical properties.

Crystal structure optimisation
The crystal structure of Tp-Azo was assigned to a hexagonal P6/m space group, with eclipsed stacking of planar layers separated by a distance of 3.3 Å, on the basis of powder XRD measurements ( Fig. 1(a) and Supplementary Tab. S1(a)). 45 Similarly, the structure of DAAQ-TFP has been assigned to a P6/m space group from Pawley refinement of powder XRD patterns (Fig. 1(c) and Supplementary Tab. S1(c)). 46 Upon geometry optimisation, in both cases we find both a breaking of the planarity within layers through an undulating distortion, as well as relative coplanar displacements between layers as shown in Fig. 1(b) and (d). The space group symmetry lowers to P1. The interlayer distance of Tp-Azo it decreases from 3.30 Å to 3.23 Å, and in DAAQ-TFP it decreases from 3.60 Å to 3.33 Å during geometry relaxation from the reference structures. The layer shift of Tp-Azo along a is -2.63 Å and along b is 2.01 Å, the offset between neighbouring layers along the ab plane is 2.73 Å. The layer shift of DDAQ-TFP in the ab plane is 2.29 Å. The horizontal offsets of both Tp-Azo and DAAQ-TFP are higher than other COFs, which have been reported with offsets of 1-2 Å between neighbouring layers. 33,34

Binding between layers
Due to their non-covalent interlayer interactions, the structural properties of COFs can be modified through exfoliation or tuning of interlayer distances. 30,48 Single-or few-layer COFs are an emerging class of functional materials. 13,57,57 For example, nanosheets of DAAQ-TFP show promise in battery cathodes due to shorter ion/electron migration pathways and higher ionic/electronic diffusion rates. 48 Hence, knowledge of the binding energy between layers in COFs is important for tuning their performance in device applications.
The binding between layers of Tp-Azo and DAAQ-TFP was calculated from the total energy difference between the relaxed monolayer and the bulk COFs. Due to the requirement of periodic boundary conditions, the layer distance was increased to 30 Å to ensure negligible chemical interactions between repeating layers ( Supplementary Fig. S2). 58 The exfoliated COF layers were fully relaxed with this fixed interlayer distance. After relaxation, the undulating monolayer structure became planar again, indicating this distortion to be a result of interlayer interactions (Supplementary Fig. S2). The binding energy, γ, to form the monolayer can be calculated per unit area according to: where E monolayer is the total energy of the COF monolayer, N is the total number of atoms in the surface of the monolayer, and E bulk is the bulk energy per atom. 58 DAAQ-TFP COFs are 2.5 meV/Å 2 and 3.1 meV/Å 2 , respectively. Compared with the binding energies of other 2D layered materials such as graphite (13 meV/Å 2 ) and MoS 2 (20 meV/Å 2 ) 58 , Tp-Azo and DAAQ-TFP can be classified as easily exfoliable 2D materials (specifically, their binding energies are smaller than 30 meV/Å 2 ). 59,60 The high porosity in the COF structure greatly contributes to this low interlayer binding energy, with γ increasing to 8.5 meV/Å 2 and 9.2 meV/Å 2 , respectively, when the pores are omitted from the framework surface area A in Equation 1.

Potential energy surface for layer displacements
A series of displaced structures were generated with the layers offset to varying amounts along the a and b axes (Fig. 2). For each displaced structure, the internal geometry was relaxed and the energy minimum was set to 0 (Fig. 3). The PES exhibits a characteristic hexagonal shape for both COFs, resembling a "sombrero" potential. 34,61 A similar scan of rigid layers without geometry relaxation is shown in Supplementary Fig. S3; a steeper and more fragmented PES is produced. The interlayer π-π interactions give rise to a stable hexagonal ring (dark blue in Fig. 3(a) and (c)) where the relative layer displacements maximise the attractive electrostatic interactions.
Eclipsed stacking of layers is significantly less energetically favourable than the displaced arrangements and represents a local maximum on the PES. The center of the PES, corresponding to no displacement, is 0.20 eV/nm 2 higher than the minimum en-ergy for Tp-Azo, and 0.17 eV/nm 2 for DAAQ-TFP. This is shown most clearly from the 1D cross-sections in Fig. 3(b) and (d). The width of the low energy wells is approximately 1.8 Å along both the x and y axes.
The suggested behaviour is distinct from typical stacking faults associated with discrete local minimum configurations, e.g. mixtures of hexagonal (AB) and cubic (ABC) packing in close-packed crystals. Here, a continuous range of configurations are accessible. Random sampling of the low energy ring would produce an average structure that appears as eclipsed to macroscopic measurements, 35 yet in reality comprises locally offset COF layers. Moreover, the soft "sombrero" PES suggest a high sensitivity of the actual COF structures to the synthesis and processing conditions.

Electronic band gap opening
Next, we consider the impact of these displacive instabilities on the underlying electronic structure of the COFs. Remarkably, the band gap variation follows the inverse of the PES. The smallest band gap is exhibited by the eclipsed structure with no displacements along the a and b axes (at the center of the heatmaps). The HSE06 calculated band gap is 0.28 eV for eclipsed Tp-Azo and 1.29 eV for eclipsed DAAQ-TFP. A band gap opening of 1.37 eV (to 1.65 eV) and 0.75 eV (to 2.04 eV) is found for displaced Tp-Azo and DAAQ-TFP, respectively. A similar behaviour has previously been observed in COF-5. 7 We note that monolayers of Tp-Azo and DAAQ-TFP exhibit even larger band gaps of 2.06 eV and 2.36 eV as a result of quantum confinement in the 2D sheets ( Supplementary Fig. S4).
These band gaps suggest that Tp-Azo and DAAQ-TFP COFs are semiconducting materials. Fig. 4(b) and (d) show that the band gaps change sharply upon small displacement between layers within a deep well of 2 Å width. However, when the displacement is more than 2 Å but less than 6 Å, the band gap oscillates between 1.25 eV and 1.72 eV for Tp-Azo, and 0.71 eV and 0.97 eV for DAAQ-TFP. Our analysis suggests a relatively small variation within the low energy ring that should be populated at room-temperature in thermal equilibrium. The magnitude of the band gap plays an important role in battery applications. It is connected to the open-circuit voltage set by the electrochemical potentials of the anode and cathode, and also connects to charge transport (electrodes must conduct ions and electrons) as well as the stability windows. 62,63

Origins of strong electronic coupling to layer displacements
The electronic band structures of the eclipsed and slipped Tp-Azo and DAAQ-TFP COFs are compared in Fig. 5. Both COFs exhibit low band dispersion along the Γ-M-K-Γ path in reciprocal space, which corresponds to in-plane directions. The layer stacking direction, which is the shortest axis in real space, corresponds to the longer Γ-A line in the band structure.
For eclipsed stacking, the interlayer interactions produce dispersive bands with a band width 1.72 eV in the upper valence band (VB) and 1.56 eV in lower conduction band (CB) along the Γ-A path of Tp-Azo. The dispersion of DAAQ-TFP is slightly reduced, giving rise to a band width of to 1.35 eV in the VB and 0.47 eV in the CB in the interlayer direction. Eclipsed Tp-Azo and DAAQ-TFP both have strongly indirect band gaps arising from the interlayer interactions between between Γ and A.
Layer slippage results in a pronounced change in the band dispersion. The band structures remain weakly indirect between Γ and A, but the dispersion itself is inverted with the VB maximum changing location. Both the conduction and valence bands become much flatter upon layer displacement, particularly along the Γ-A path. In going from the eclipsed to displaced stacking mode, the widths of the topmost valence bands reduce from 1.72 to 0.53 eV (Tp-Azo) and from 1.35 to 0.28 eV (DAAQ-TFP), and from 1.65 to 0.15 eV (Tp-Azo) and 0.53 to 0.13 eV (DAAQ-TFP) for the bottom-most conduction bands.
The corresponding Γ point wavefunctions are shown in the insets of Fig. 5. They confirm that the band edges are formed from the C 2p z π subsystem. For the eclipsed structure, the interlayer interactions are strongly anti-bonding at the Γ point. This explains the strong downward dispersion towards A, where the phase of successive layers is reversed. In the displaced structure, stronger interlayer π bonding interactions are allowed at the Γ point and the band dispersion is suppressed along the Γ-A line. These changes result in a band gap that is weakly indirect and much larger in magnitude compared to the eclipsed structure.

Conclusions
It is convenient to represent and model covalent organic frameworks as an ordered sequence of eclipsed planar layers. However, by taking the examples of Tp-Azo and DAAQ-TFP, we have shown that the results can be misleading in line with recent observations for other layered COFs. A displaced stacking sequence of undulating layers both lowers the total energy of the frameworks and results in a large change in the electronic structure driven by interlayer π orbital overlap. Layer displacements produce a pronounced band gap opening in these frameworks. The unusual "sombrero" potential energy surface for layer displacements, which mirrors the variation in band gap, has important implications. Although macroscopically a given COF may appear to have an eclipsed structure, for example on the basis of diffraction measurements, locally a continuous range of stacking sequences are accessible. The strong coupling between layer orientation and electronic structure highlights the potential for COF twistronics where longer range modulations in the crystal potential are harnessed. Layer offsets may be controlled by various experimental approaches, such as chemical intercalation, synthetic modification of composition including aromatic ring size, and processing temperature. 36,37 These findings will be of particular importance when screening COFs for applications in energy storage and conversion where electrochemical and photochemical descriptors are significantly altered including accessible voltage ranges for batteries, stability windows for electrocatalysis, and visible light absorption for photoelectrochemical systems.

Conflicts of interest
The authors have no conflicts of interest to declare.