Exquisite control of electronic and spintronic properties on highly porous Covalent Organic Frameworks (COFs): transition metal intercalation in bilayers

Two-dimensional covalent organic frameworks (2D COFs) are crystalline organic porous materials stacked in a layered fashion. In general, these materials have excellent structural tunability, which can be achieved through the various tools of organic synthesis. Their layered and porous nature makes them attractive candidates for electronics, optoelectronics, and catalysis. However, their application is still limited due to relatively poor π-delocalization and practical applications require controlling and tuning their electronic structure. In this paper, using hybrid density functional theory, we computationally explore a novel 2D COF architecture, consisting of only two crystalline atomic layers made of benzene, boroxine, and triazine rings. We study the intercalation of first-row transition metals in the bilayer to enhance and fine-tune their electronic and magnetic behavior. This resulted in the development of one pristine bilayer, 63 intercalated bilayers, and one trilayer 2D COF. We found that the concentration and position of transition metals in the structure can drastically change the 2D COFs’ electronic, magnetic, and spintronic features. Based on their spin-polarized electronic properties, these transition metal-intercalated 2D COFs have potential applications as water splitting catalysts, direct semiconductors in the visible range, half metals, half semiconductors, and bipolar magnetic semiconductors.


Introduction
We report a first-principles computational analysis of transition metal (TM) intercalation in two-dimensional (2D) triazine-benzene-boroxine bilayer Covalent Organic Frameworks (COFs), providing a mechanism to tailor their electronic and magnetic properties. Due to the wide range of spin-polarized bandgaps and band alignments observed in these structures, we believe these materials have the potential to perform as two-dimensional porous optoelectronic, catalytic, and spintronic devices. Analysis of different magnetic states reveals that these structures can be stable in ferromagnetic (FM), ferrimagnetic(FiM), antiferromagnetic(AFM), and nonmagnetic(NM) arrangements, depending on the type and amount of TMs intercalated in the bilayer COF. In turn, their magnetic structure affects their potential electronic/spintronic properties, which we use to our advantage for fine-tuning their electronic structure.
COFs are polymeric porous materials that form extended network structures connected by organic linkers. These materials enable the integration of organic, light-element (e.g. C, H, N, B, O) building blocks with strong covalent bonds. Using dynamic covalent chemistry and slower reaction kinetics, one can obtain 2D or 3D crystalline periodic materials [1,2]; otherwise, we would get amorphous solids. In COFs, the faces and edges of their monomer units are exposed to pores, resulting in high surface area. Their structural, physical, and chemical features include tunable porosity, low density, high crystallinity, and strong covalent bonds. These properties establish these materials as an ideal platform for applications including gas storage [3,4], sensors [5], catalysis [6], energy storage [7,8], and optoelectronics [9].
In particular, 2D COFs are a family of novel crystalline materials comprised of covalently-linked layers stacked via non-covalent long-range interactions such as van der Waals (vdW) forces, Londondispersion forces, and aromatic interactions [10]. 2D COFs can be designed to a highly periodic atomic framework and give the uniaxial stacking of 2D macromolecular frameworks on an atomic level [11,12]. Furthermore, stacking interactions in these materials are essential to control their electronic properties [13], stability [14], crystallinity [15,16], and porosity [14,15]. The tunability of 2D COFs offers unique architectures that provide a wide array of electronic properties given their spatial stacking and interlayer interactions [17]. Moreover, their nanoporous and layered arrangement leads to an ultra-high surface area, offering a large number of binding sites for molecular sensing [18], gas storage [19], and catalytic reactions [20,21].
Traditionally, 2D COFs are formed by dynamic covalent linkages, such as imine and arylhydrazone bonds. These exhibit relatively poor π-delocalization and low stability, which arise from structural constraints (steric hindrance), functional groups that disrupt π conjugation, weak π-π interactions between adjacent aromatic rings, defects & imperfections, and large energy gaps due to the choice of building blocks; impeding their usage as conductors or semiconductors. However, their organic structure-and consequently, their electronic properties-can be modified or designed using the vast tools of organic synthesis. In recent years, researchers have been able to promote electron delocalization by generating conductive 2D COFs with electroactive organic building blocks and through chemical doping [22,23]. However, a technique that has not been widely explored for the modification of 2D COF electronic structures is intercalation.
Intercalation is the process of inserting guest species into the interlayer gap of 2D materials. It has been shown to provide methods to change layered materials' electronic band filling and Fermi level [24]; and can even induce exotic properties such as superconductivity [25,26] and enhance spin-orbit coupling gap opening in topological insulators [27]. Particularly, using transition metals as intercalants in 2D materials has resulted in interesting property-tuning mechanisms for electronic/optoelectronic, spintronic, magnetic, and topological insulator materials [27][28][29]. Although transition metal intercalation has been extensively explored for other 2D materials [30], in the particular field of 2D COFs, this mechanism has not been widely studied, with only a few publications on the matter. In the past, computational simulations have been carried out for intercalation in 'bulk' 2D COFs, i.e., structures with periodic boundary conditions, where layers repeat 'infinitely' in the z-direction [31]. This previous work predicted intercalation to be an excellent tool for modifying and enhancing COFs' electronic properties.
Here, we specifically deal with transition metal intercalation in an isolated, freestanding triazine-benzeneboroxine 2D COF bilayer, a material which-to the best of our knowledge-has not been conceptualized in the past. Our main hypothesis is that by controlling the intercalated d-block elements and concentrations, we will be able to modify their interlayer interactions (by promoting interactions between the COF π-orbitals and TM d-orbitals), which can lead to fine-tuning their electronic and magnetic behavior. Such changes include, for example, modifying the 2D COF's natural insulating behavior to that of a semiconductor by promoting the appearance of available states near the Fermi level via COF p-orbital and TM d-orbital interactions (illustrated in figure 1). Furthermore, electron confinement in two spatial dimensions is known to drastically change the electronic properties of materials when they are reduced from bulk to 2D [32][33][34][35][36]. Therefore, bilayer 2D COFs' electronic properties are expected to be drastically different from their previously-explored bulk counterparts [31,37].

Results and discussion
2.1. Pristine triazine-benzene-boroxine bilayer 2D COF Previous publications demonstrate the stability of 2D COFs comprised of triazine, benzene, and boroxine rings. Prakhira S et al theoretically predicted the stability of this 2D COF architecture in its bulk configuration [31], and Gontarczyk K et al experimentally synthesized this structure (labeled BTA-COF4 in their publication) [38]. In these publications, the 2D COF was calculated/synthesized as a bulk structure comprised of many stacked layers. However, in our work, we focus on a freestanding 2D COF bilayer, which has not been theoretically predicted or experimentally synthesized in the past. This material is expected to present different electronic properties compared to its bulk counterpart, due to electron confinement in two dimensions. Figure 2 shows the structural features and electronic properties of the pristine bilayer triazine-benzene-boroxine 2D COF. Geometry optimization using unrestricted hybrid Density Functional Theory with Grimme -D3 dispersion corrections (HSE06-DFT-D3) reveals the following structural parameters:  Figure 2 shows a top view of the optimized structure, where the eclipsed 2D COF stacking can be observed. This figure also shows the spin-polarized band structure and atom-projected density of states (DOS) (middle and right-hand side of the figure), showing that the distribution of spin-up and spin-down electrons is symmetrical in this material, which is expected from its non-magnetic configuration. The material presents a bandgap of 3.40eV, indicating insulating electronic behavior. Comparatively, its bulk counterpart was predicted to have a smaller bandgap of about 2.6 eV [31]. Furthermore, the DOS reveals that carbon atoms are the main contributors to the electronic states around the Fermi level, i.e. near the conduction band minimum (CBM) and valence band maximum (VBM). The dynamical stability of the pristine bilayer 2D COF was confirmed through vibrational calculations at the Γ, M, and K reciprocal lattice points, which revealed no imaginary frequency modes, indicating a stable structure.
The pristine optimized structure was taken as a starting point to generate the intercalated 2D COFs. Transition metals were positioned at the centroid of the bilayer COF rings in different concentrations. Intercalants and concentrations will be referred to as TM-01 through TM-05, where TM = [Co, Cr, Cu, Fe, Mn,  Ni, Sc, Ti, and V]. There are two variations of TM-01 and TM-04, which we distinguish by adding 'a' or 'b' to the name. Therefore, we generated seven different atomic configurations each with nine TM intercalants, resulting in a total of 63 intercalated bilayer 2D COFs. Configuration TM-01a has one transition metal per unit cell, positioned at the centroid of boroxine rings. TM-01b also has one transition metal per unit cell, positioned at the centroid of triazine rings. TM-02 has two transition metals per unit cell, placed at the centroid of boroxine and triazine rings. TM-03 has three transition metals per unit cell, located at the centroid of all benzene rings. TM-04a has four transition metals per unit cell, intercalated at the centroid of all benzene and boroxine rings. TM-04b also has four transition metals per unit cell, intercalated at the centroid of all benzene and triazine rings. Finally, TM-05 has transition metals located at the centroid of all rings in the structure (triazine, benzene, and boroxine).

Magnetism
After generating the sixty-three intercalated structures, we performed a preliminary structural optimization. This procedure determined if our compounds possess spin-polarized states, and gave us a starting geometry to work with. During these early optimizations, the TM atoms' electronic states were populated in a ferromagnetic alignment, and their spin was locked for the first 3-5 iterations of the first self-consistent field (SCF) cycle. After this initial spin locking, the system was left unconstrained and allowed to drift toward the magnetic state that minimized its energy.
We further investigated the structures' magnetic behavior starting from their pre-optimized geometries (obtained from the discussion above). In order to find the most stable magnetic structure at each transition metal concentration, we tested all magnetic states that were compatible with the P6 m2 symmetry of the system. These magnetic states include NM, FM, and AFM/FiM configurations for some structures. Table S1 in the supplementary information shows the different magnetic configurations tested at each TM concentration. For every magnetic configuration, we ran one single-point energy (self consistent field) calculation on the preoptimized geometries, where the TM magnetic states were spin-locked for the first 10 iterations of the cycle. After the single-point calculations, if an AFM, FiM, or NM structure was more energetically favorable than its initial FM configuration, we computed another geometry optimization for that particular magnetic structure. Even though the atomic structure remains the same, changing a structure's magnetic configuration can result in vastly different electronic properties. Therefore, it was crucial to find the most stable magnetic state in our materials. Table 1 shows the preferred magnetic configuration for all intercalant concentrations calculated in this work. To analyze the magnetic configurations, we use α − β Mulliken Population Analysis (MPA). This analysis revealed that the magnetic state of intercalated TMs is dependent on the concentration: Co and V intercalations are FM for configurations TM-01 through TM-04, but FiM for the TM-05 configuration. Cu intercalation is non-magnetic throughout all configurations. Finally, Cr, Fe, Mn, Sc, and Ti intercalations can be FM, AFM, or FiM, depending on TM concentration. Furthermore, in these compounds, the magnetic behavior of intercalated atoms can be highly dependent not only on their concentration, but also on their position. In particular, we highlight that some TM atoms have reduced unpaired electron density depending on their position. For example, for Cr intercalation, the α − β (unpaired) electron population is larger for Cr atoms intercalated at the centroid of boroxine or triazine rings, compared to those intercalated at the centroid of benzene rings, with the exception of Cr-02 where magnetic moments are reduced in triazine rings. In the Ni-04 and Ni-05 configurations, MPA indicates reduced spin density for Ni atoms intercalated in boroxine and triazine rings, compared to benzene rings. Sc-02 and Sc-04 have a reduced number of unpaired electrons for atoms intercalated at the centroid of triazine rings, and Sc-03 has a very small amount of unpaired electron density for all Sc atoms. Finally, for Ti intercalation, the only atoms that show open-shell configurations are those intercalated at the centroid of boroxine rings.
It is important to note that, even though MPA gives a good indication of the type of magnetic configuration (FM, FiM, AFM, NM) for the intercalated TM atoms, it is known to yield inaccurate results for the number of unpaired electrons due to an overestimation of bonds' covalent character. Therefore, we use this tool to qualitatively find if atoms have closed-shell/open-shell character and if they present smaller or larger unpaired electron density relative to other atoms in the structure. Magnetic moment calculations would be necessary to numerically determine the atoms' spin.

Structural properties and binding energy
In this subsection, we focus on the structural properties of the optimized intercalated 2D COFs. Structural features, including average bond distances and average angles for the sixty-three intercalated bilayers, are displayed in table S2 of the supplementary information. To summarize the information in table S2, we find that C−C, B−O, and C− N average distances range from 1.40 to 1.42 Å for d C−C , 1.38 to 1.42 Å for d B−O , and 1.33 to Table 1. 2D COFs' magnetic configuration, interlayer distance (ILD), α and β electron bandgaps (E g a and E g b , respecively), and bandgap type: conductor (C), indirect (I), or direct (D).

Intercalant
Magnetic configuration Similar trends as the ones described above can also be seen in the average angles between atoms. For benzene rings, we find that C-C-C average angles (∠ C−C−C ) are equal or slightly larger for intercalated structures compared to the pristine COF, with Ni-01b, Ni-02 and Cu-05 being the only outliers. Similar to the C-C bonding distances, these angles are noticeably larger in configurations TM-03, TM-04, and TM-05 since these configurations have intercalants close enough to the benzene rings to cause distortions. On the other hand, the angles in boroxine rings suffer a strong rearrangement in all intercalations, except TM-01b, TM-03, and TM-04b, following the same general trend as B−O distances (TMs are far enough from boroxine rings to avoid distorting the lattice significantly). Finally, the angles in triazine rings are generally slightly distorted for all configurations, but distort the most when intercalating atoms at the centroid of triazine rings (TM-01b, TM-02, TM-04b, TM-05), which is a similar behavior to that observed for C−N distances.
Next, we pay attention to two important features that describe TM interactions with the 2D COF: (1) binding energy (BE) of transition metal intercalants in the 2D COF structure and (2) average interlayer distance (ILD) between the two COF layers. The third column of table 1 shows the average interlayer distance (ILD) between COF layers for all intercalant concentrations. As a reference, the pure 2D COF has an ILD of 3.34 Å. Furthermore, in figure 3 we plot the relative binding energy per intercalant for all 2D COF TM intercalations calculated in this work.
Binding energy was estimated using the following equation: is the BE per transition metal atom intercalated in the structure, E TM2DCOF is the total energy of the optimized TM-intercalated 2D COF, E layers is the energy of the 2D COF layers without intercalants (keeping the same ILD as in E TM2DCOF ), and E intercalants is the energy of the transition metals. Due to the negative sign on the right side of equation (1), more positive BEs indicate stronger interactions. From figure 3, it can be observed that TM-03 configurations (intercalation at the centroid of benzene rings only) in general have larger BE than other configurations, with the only exception being Fe, which has equal BE for TM-03 and TM-04b. Furthermore, from table 1, we observe that-for most configurations-TM-03 has a very small change in ILD with respect to TM-02, and in some cases is even smaller than TM-02, even though it has more intercalated atoms in its structure. This indicates high stability and stronger interactions for intercalants at the centroid of benzene rings. Another important trend observed throughout all configurations is that TM-01a and TM-04a are less energetically favorable than their TM-01b and TM-04b counterparts, which can be observed from their binding energies. This trend is also evident in the ILD (see table 1), which is smaller when intercalating atoms at the centroid of triazine rings ('b' configurations) compared to boroxine rings ('a' configurations), indicating stronger interactions for TM-01b and TM-04b. From these results, we predict that-in general-during intercalation the centroid of benzene rings will be populated first, followed by triazine rings, and boroxine rings will be populated last.
From table 1 column 3, it can be observed that the ILD, in general, was modified for all TM-intercalated 2D COF materials with respect to the pure 2D COF (3.34 Å), due to the interaction between 2D COF layers and the TM atoms. Initially, we expected the ILD to increase proportional to the number of intercalated atoms; however, this was not the case. For structures with Cr, Ti, and V intercalants, we observe a trend where the TM-03 configuration has an equal or smaller ILD than TM-01 and TM-02 configurations. Even though TM-03 COFs have more atoms intercalated in the bilayer, their position at the centroid of benzene rings (instead of boroxine or triazine rings) reduces their ILD. This reduced distance indicates highly favorable binding between the TMs (Cr, Ti, and V) and benzene rings. On the other hand, structures intercalated with Cu atoms follow the expected proportional increase with the number of intercalated atoms, except for Cu-02 and Cu-03, which have the same ILD of 3.44 Å. A similar phenomenon occurs in Mn-intercalated structures, where Mn-03 (3.45 Å) presents very similar ILD as Mn-02 (3.44 Å). Furthermore, Fe-intercalated structures present an increasing ILD for Fe-01, Fe-02, and Fe-03 configurations; however, there is a decrease from Fe-03(3.46 Å) to Fe-04a(3.39 Å). The only difference between Fe-03 and Fe-04a configurations is an additional Fe atom intercalated at the centroid of boroxine rings; therefore, the decrease in ILD leads us to propose that Fe has a similar affinity to benzene, boroxine, and triazine rings. This is reflected in the binding energies, which are very similar for configurations Fe-03, Fe-04, and Fe-05. This behavior also occurs in Sc-03 and Sc-04a, where the ILD decreases by adding a TM at the centroid of boroxine rings. Interestingly, the Fe and Sc intercalation also have similar binding energy trends for the different concentrations intercalation patterns. On the other hand, Mn-and Ni-intercalations are the only groups that follow the initally-expected trend, where the ILD is proportional to an increasing number of intercalated atoms. Another interesting result is that the ILDs of configurations Co-01b (3.33 Å), Co-04a , and V-04b (3.29 Å) are smaller than the ILD of the pristine structure, suggesting weaker steric repulsions and potentially indicating high affinity between TMs and COF layers for these configurations. Finally-in terms of intercalants moving in the x − y plane-we observe that for most structures, the TMs tend to stay at or near the centroid of the rings. Cu, Ni, and Sc are the only outliers (see figures S3, S6, and S7 in the supplementary information). The Cu atoms at the centroid of benzene rings tend to drift toward triazine rings, and the Ni & Sc atoms drift toward boroxine rings.

Electronic properties
In this section, we highlight the electronic properties of our intercalated 2D COF bilayers. Our main initial hypothesis was that, by changing the type of transition metal and its concentration, it is possible to fine-tune the bilayer COFs' electronic properties. Table 1 columns 4-7, show the spin-polarized bandgap values and their type (direct, indirect, conductor) for all 2D COFs calculated in this work. As an initial observation, some of our materials (Co-01a/b, Co-03, Co-05, Cr-01a/b, Cr-04a, Cu-01a/b, Cu-03, Cu-05, Fe-03, Fe-04b, Fe-05, Ni-03, Ni-04a/b, Ni-05, Sc-01a/b, Sc-04b, Ti-01a, Ti-02, and Ti-05) have conducting electronic behavior in at least one of their spin channels. Therefore, our studies confirm that intercalation can be used to enhance the conductivity of these materials. This is a drastic change from the pure 2D COF, which presented insulating behavior (3.40 eV bandgap) and significant energy gaps below the Fermi level. The data in table 1 also indicates that most of our 2D COFs are semiconducting; moreover, our spin-polarized calculations reveal that most of these materials have different bandgaps for α(spin-up) and β(spin-down) states. To further explore the intercalated 2D COFs' electronic properties, we plot the band structures and DOS of all intercalated structures aligned with respect to vacuum. These plots can be found in the Supplementary Information figures S1-S9). Through these plots, we observe a general trend where the addition of TMs in the structure reduces not only the bandgap, but also the electron density gaps within the structures (below the valence band minimum). We will further expand on how this affects the potential applications of these materials in the following subsections.
To obtain an overview of the electronic properties of all semiconducting and insulating COFs, we plot the spin-polarized absolute band-edge alignment for materials with both spin-up and spin-down bandgaps greater than zero (see figure 4). This plot allows us to easily see spin-polarized conduction band minima and valence band maxima; furthermore, helping to detect which of these compounds are potential candidates for catalytic, optoelectronic, and spintronic applications. Materials with a VBM above SHE and a CBM below SOE can generate electrons/holes with favorable energies to catalyze the water-splitting reaction. The CBM for all of our semiconducting materials falls within a range of −1 to −2 eV versus SHE (refer to energy axis on the right side of the plot), which is more negative than that of H + /H 2 (0 eV), resulting in the ability to reduce hydrogen through photoexcited electrons in these materials. On the other hand, only one of our materials (Ni-01a) has its VBM for both spin channels below the oxygen evolution level (1.23 eV versus SHE), making this material a potential candidate for oxygen evolution catalysis. Furthermore, for this intercalation pattern, we observe that the Ni d-band center for Ni-01a is relatively close to the Fermi level compared to other Ni-intercalations, which is known to be a descriptor of strong binding interactions with reactants and enhanced catalytic activity [39,40].
The wide range of electronic properties achieved through intercalation allows for versatility in our 2D COFs' applications. Next, we move to optoelectronic applications, where we pay specific attention to direct bandgap semiconductors. In direct semiconductors, electronic excitation from valence to conduction band does not require electrons' momentum to change, e.g., via phonon interactions. Therefore, this type of semiconductor is preferred for applications that require efficient photon absorption and emission, for example, LEDs, solar cells, and photodetectors. From the information in table 1, we see that twenty-eight out of the sixty-three intercalated 2D COF structures have direct bandgaps in at least one of their spin channels (Cr-02, Cr-03, Cr-04a, Cr-04b, Cr-05, Cu-04a, Fe-01, Fe-04a, Fe-04b, Mn-01a, Mn-03, Mn-04a, Mn-04b, Mn-05, Ni-04a, Sc-02, Sc-03, Sc-04a, Sc-04b, Sc-05, Ti-03, Ti-04a, Ti-04b, V-01a, V-03, V-04a, V-04b, V-05). It is important to note that in all of these cases, the direct bandgaps correspond to the smallest spin-polarized bandgap value. Our direct semiconducting materials have a wide range of bandgap values, from 0.20 to 2.02 eV. Three of these materials have the ability to absorb and emit photons in the visible range: Cr-03, Ti-03, and V-03, since their smallest spin-polarized bandgap value falls between 1.65 and 3.26 eV (visible light spectrum). Interestingly, these materials have a similar structure, as they present a TM-03 intercalation pattern (with the intercalants positioned at the centroid of all benzene rings). The spin-polarized band structures and projected DOS for these materials are plotted in figure 5. In this figure, we observe the important role of the transition metals, which provide electron density near the Fermi level from their 3d-orbitals. Transition metals and carbon atoms provide the largest electron density at the CBM/VBM. Thus, it is evident that the carbon 2p-orbitals in the COF layers hybridize with the TM d-orbitals, providing the unique direct-semiconducting electronic structures observed. These results confirm our hypothesis that promoting p − d orbital interactions between the COFs and TMs leads to a modification in the electronic structure and properties of these systems. Furthermore, Ti-04a, V-04a, and V-01a-which are also direct semiconductors-have their larger spin-polarized bandgaps in the visible range and, although less likely, these visible-range direct transitions can also occur for these materials.
From table 1, we also observe that eighteen of our compounds are semiconducting/insulating in one of their spin channels and conducting in the other, i.e. they present half-metallic behavior (Co-03, Cr-01a, Cr-01b, Cr-04a, Cu-01a, Cu-01b, Fe-04b, Fe-05, Ni-03, Ni-04a, Ni-04b, Ni-05, Sc-01a, Sc-01b, Sc-04b, Ti-01a, Ti-02, Ti-05). Thus, intercalation in 2D-COFs can induce spintronic half-metal behavior. Due to their conductivity in only one spin orientation, half metals have applications in spin-based electronics, for example, as spin-injectors in spin-based transistors and diodes [41,42]. Moreover, three of these compounds (Cr-04a, Fe-04b, and Ni-04a) have a direct bandgap in their semiconducting/insulating channel, while the rest are indirect. The half-metal bandgaps range from 0.52eV to 3.58 eV, with several intermediate values. Thus, it is possible to achieve a wide range of bandgaps, depending on the concentration and position of intercalants, opening their application in next-generation spin-based electronics. As an example, in figure 6 we plot the structure and spin-polarized electronic properties of 2D-COF-Ni-05. This material has conducting spin up states, while the spin down channel shows a semiconducting 1.32 eV bandgap. We observe that the spin up conducting states mainly come from carbon contributions. The spin down semiconducting channel, however, mostly presents carbon contributions at the CBM, and the VBM is primarily comprised of nickel and nitrogen orbitals. The interactions between COF d-orbitals and TM p-orbitals again become essential to provide states that result in enhanced electronic/magnetic properties, in this case, half metal behavior.
Spin-up and spin-down band edges are plotted separately in figure 4, allowing us to easily observe which materials are potential candidates for other spintronic applications. Within the semiconducting 2D COFs that were analyzed, we find some materials with off-shifted spin-up and spin-down band edges, i.e., spin-up and spin-down VBM and CBM have different relative positions. Within our intercalated structures, we find five materials that clearly present half-semiconducting character: Co-04a, Co-04b, Cr-02, Mn-04a, and V-01a. These materials present a wide bandgap in one of their spin channels and narrow bandgap in the other spin channel. Half semiconductors have their VBM and CBM in the same spin channel, allowing them to produce 100% polarized electrons through photon or thermal excitation. To exemplify this half-semiconducting character, in figure 7 we plot the band structure and DOS for Mn-04a. This material presents a direct narrow bandgap of 0.57 eV in its spindown channel and an indirect wide spin-up bandgap of 2.1 eV in its spin-up channel, making spin-down electron transitions much more likely to occur not only due to their bandgap values but also their bandgap type (small direct  bandgap vs wide indirect bandgap). Therefore, this material is an ideal candidate for the generation of 100% spindown polarized electrons. Similar to our optoelectronic and half-metal materials, this spintronic behavior is heavily mediated by Mn 3d-orbitals, which are hybridized with carbon, boron, and oxygen 2p-orbitals and are important contributors to the DOS near the band edges. Figure 4 also reveals that the band alignments for eleven of our materials have potential as bipolar magnetic semiconductors (BMS): Cu-02, Fe-01a, Fe-01b, Mn-01a, Mn-01b, Sc-02, Ti-04a, Ti-04b, V-02, V-04b and V-05. These spintronic semiconducting materials have their VBM and CBM fully polarized in the opposite spin direction, with their spin-up and spin-down channels off-shifted. In BMSs, it is possible to generate fully spinpolarized electrons by applying a gate voltage. Taking V-05 as an example (see figure 8(a), under zero gate voltage it presents semiconducting behavior. Further, at an appropriate negative gate voltage, the Fermi level is shifted down (in the vacuum scale), resulting in a conductive spin-up channel and semiconducting spin down-channel. Additionally, at an appropriate positive gate voltage, the Fermi level is shifted up, resulting in a conductive spindown channel and a semiconducting spin-up channel.
Even though 2D-COF-V-05 presents appropriate band edges for a BMS, its density of states has large gaps within the band structure, making it difficult for electrons to travel throughout the material. To remedy this issue, it is necessary to increase the available electronic states in the system. Therefore, as proof of concept, we  added another fully-intercalated V-05 2D COF layer. We present these results in figure 8(b). In this plot, we observe a rearrangement of electronic states when adding a third intercalated COF layer, where now a positive gate voltage would generate a conductive spin-up channel and a negative voltage would generate a conductive spin-down channel. Most importantly, the large gaps below the VBM and above the CBM are substantially reduced, which results in enhanced available electronic states within the material. This is an important result, as we see that spintronic behavior in 2D COFs can be fine-tuned by the number of layers, which is an effect that has not been researched in the past for this family of compounds.
We conclude this section for electronic and spintronic properties by mentioning that 2D-COF-Mn-05 presents a peculiar band structure, with spin-down bands that resemble a double Dirac cone near the Fermi level (see figure S5). A similar phenomenon has been theoretically predicted in S-Graphene [43]. With the Mn-05 intercalation, we potentially find a non-trivial topological material. However, analyzing this phenomenon accurately requires the introduction of relativistic effects (namely, spin-orbit coupling), which we have not considered here. Thus, we leave this for future work.

Conclusions
In this work, we performed a comprehensive hybrid-level computational study of transition metal intercalation in triazine-benzene-boroxine 2D COFs for catalytic, magnetic, and electronic/spintronic property enhancement. We first describe the structural and electronic properties of the pristine bilayer 2D COF. This structure presents insulating electronic character (bandgap of 3.40 eV). In this material, carbon atoms are the main contributors to the electronic states surrounding the Fermi level. Next, we generated sixty-three intercalated bilayer structures with nine first-row transition metals. Our analyses reveal that the type, concentration, and position of intercalants changes their structural, magnetic, and electronic features. Structurally, intercalation results in changes of the interlayer distance, which can be directly related to interaction strength between the transition metals and the bilayer. The structural analysis was complemented with a binding energy analysis, where we observe that transition metal intercalation is energetically preferred at the centroid of benzene rings, followed by triazine and boroxine. This indicates that-generally-during intercalation, the centroid of benzene rings will be populated first, followed by triazine, and finally boroxine. Magnetically, we observe a wide range of stable spin configurations (FM, FiM, AFM, NM), which depend on the type, position, and concentration of intercalants. This wide range of magnetic features allows to fine-tune and enhance 2D COFs' features for photocatalytic, electronic, and spintronic applications. In general, we find that transition metal d-orbitals have important contributions to the density of states around the Fermi level and can even promote the contribution of C, N, B, and O atoms around the Fermi level.
Within our sixty-three intercalated bilayer 2D COFs, we find that twenty-four of our materials are conductors in at least one of their spin channels, while the rest are semiconductors or insulators with a bandgap ranging from 0.2 to 3.24 eV. In terms of photocatalysis, our absolute spin-polarized band edge alignment shows that our semiconducting materials have a band alignment suitable for hydrogen reduction through photoexcited electrons, and one of our materials (Ni-01a) has a band alignment suitable for oxygen evolution photocatalysis. Within the semiconductor COFs, we find three (Cr-03, Ti-03, and V-03) that have direct bandgaps that fall within the visible range; therefore, they are ideal for optoelectronics. Finally, we explore the spintronic properties of our materials, where we observe that eighteen of them (Co-03, Cr-01a, Cr-01b, Cr-04a, Cu-01a, Cu-01b, Fe-04b, Fe-05, Ni-03, Ni-04a, Ni-04b, Ni-05, Sc-01a, Sc-01b, Sc-04b, Ti-01a, Ti-02, and Ti-05) present half metal behavior, five of them (Co-04a, Co-04b, Cr-02, Mn-04a, and V-01a) have appropriate band alignments for their use as half semiconductors, and eleven of them (Cu-02, Fe-01a, Fe-01b, Mn-01a, Mn-01b, Sc-02, Ti-04a, Ti-04b, V-02, V-04b and V-05) have potential application as bipolar magnetic semiconductors. However, bilayer BMSs present large gaps within their band structure below the VBM, which inhibits the movement of electrons within the material. To remedy this issue and as proof of concept, we added another fully intercalated COF layer to the V-05 structure, which provided more intermediate electronic density levels while preserving its bipolar magnetic semiconducting character. In general, this work provides a new architecture of isolated bilayer and trilayer TM-intercalated 2D COFs, which have the potential to be used in a wide array of fields, including catalysis, optoelectonics, and spintronics.

Materials design
The pristine 2D bilayer COF was prepared with boroxine (B 3 O 3 ), triazine (C 3 N 3 ), and benzene (C 6 H 4 ) rings. It was built using two layers with P6 m2 symmetry (space group 187). In the unit cell, each COF layer is comprised of one triazine, one boroxine, and three benzene rings, and the two layers are stacked directly on top of each other (eclipsed stacking). The reason for this combination of building units is to maximize the environment of different TM atoms and organic rings, while maintaining synthesis feasibility which can explore the effects of TM atoms in the different rings. Seven different types of TM-intercalated 2D bilayer COFs were computationally prepared through the addition of TMs at the centroid of the bilayer COF rings: • 2D-COF-TM-01a. One TM per unit cell, intercalated at the centroid of the boroxine rings.
• 2D-COF-TM-01b. One TM per unit cell, intercalated at the centroid of the triazine rings.
• 2D-COF-TM-03. Three TMs per unit cell, intercalated at the centroid of all benzene rings.
• 2D-COF-TM-04a. Four TMs per unit cell, intercalated at the centroid of all benzene and boroxine rings.
• 2D-COF-TM-04b. Four TMs per unit cell, intercalated at the centroid of all benzene and triazine rings.
In these 2D COF architectures, TM = Co, Cr, Cu, Fe, Mn, Ni, Sc, Ti, and V. Therefore, a total of one pristine structure and sixty-three intercalated bilayer structures were generated. Additionally, a trilayer V-05/V-05 2D COF structure was generated by adding another intercalated COF layer to the previouslyoptimized V-05 bilayer. All magnetic configurations compatible with the P6 m2 symmetry were tested, and electronic properties were plotted for the most stable magnetic configuration.
The optimized structures for all sixty-five 2D COFs are provided in CIF format in the supplementary information. The CIF files were symmetrized using the FINDSYM software suite [44,45].

Periodic density functional theory (DFT)
The equilibrium geometries and electronic properties of all 2D layered structures were investigated through dispersion-corrected unrestricted hybrid Density Functional Theory. The exchange-correlation functional was defined through the Heyd-Scuseria-Ernzerhof 2006 (HSE06) hybrid functional [46]. Weak, long-range interactions (vdW and London dispersion forces) were taken into consideration through Grimme's third-order dispersion correction (-D3) [47]. All 2D COFs have P6 m2 hexagonal symmetry, which was constrained throughout the optimization of atomic positions and lattice parameters. All structures were fully optimized in atomic positions and lattice parameters. Calculations were performed using the CRYSTAL17 code [48], which uses atom-centered Gaussian basis sets to describe atomic orbitals, allowing an efficient implementation of Post-Hartree-Fock methods (such as hybrid DFT). To avoid self-interactions due to periodic boundary conditions in the slabs, we have considered a vacuum space of 500 Å, resulting in essentially no periodicity in the z-direction. This is possible due to the use of localized Gaussian Type Orbitals (GTOs). GTOs present an advantage over plane wave basis sets, as the latter become significantly more resource-intensive with the addition of vacuum space. Triple-zeta valence with polarization (TZVP) quality basis sets [49] were used for all atoms in the 2D COF calculations. The threshold used to evaluate the self consistent field (SCF) energy convergence, root mean squared (RMS) forces, maximum forces, RMS atomic displacements, maximum atomic displacments, and between-geometry convergence energies were set to 2.72 × 10 −6 eV, 1.54 × 10 −2 eV Å −1 , 2.31 × 10 −2 eV Å −1 , 6.35 × 10 −4 Å, 9.53 × 10 −4 Å, and 2.72 × 10 −6 eV. Vibrational modes are highly sensitive to structural parameters so we set stricter convergence criteria for frequency calculations. The stricter thresholds for SCF energy convergence, RMS forces, maximum forces, RMS atomic displacements, maximum atomic displacments, and between-geometry convergence energies were set to 2.72 × 10 −10 eV, 1.54 × 10 −3 eV Å −1 , 2.31 × 10 −3 eV Å −1 , 6.35 × 10 −5 Å, 9.53 × 10 −5 Å, and 2.72 × 10 −10 eV. Unrestricted wave functions were considered due to the partially filled 3d orbital electrons of the TM atoms intercalated in the structures. Electronic spin states were allowed to relax in the unrestricted wave function calculations after a set amount of SCF steps, described in section 2.2.
The electronic properties of the 2D bilayer COFs were calculated at the same level of theory (i.e. DFT-HSE06-D3). Integrations inside the first Brillouin zone were sampled on 8 × 8 × 1 Monkhorst-Packk-mesh grids for all 2D COFs during geometry optimization and single-point energy calculations, this is, at a resolution of approximately 2π × 1/60 Å −1 (ak a = 40-60, bk b = 40-60, k c = 1, where a, b, c are direct lattice vectors and k a , k b , k c are integer shrinking factors). Moreover, to compute the density matrix and Fermi energy, a denser 'Gilat net' sampling of k-points was carried out with a resolution of approximately 2π × 1/120 Å −1 (k k 2 G max = , where k G is the Gilat net shrinking factor and k max is the largest shrinking factor among the previously defined k a , k b , k c ). The band pathway followed the Γ − M − K − Γ high symmetry reciprocal lattice points. The electronic density of states was calculated by including the corresponding atomic orbitals for all elements (C, N, O, B, H, Co, Cr, Cu, Fe, Mn, Ni, Sc, Ti, and V) in the 2D COFs. The contributing components of the transition metal dsubshells were projected onto the total DOS plots showed in the Supplementary Information. Geometry visualization and analyses were performed using VESTA [50].
Absolute band-edge alignments were carried out following the procedure described on the 'CRYSTAL tutorials' website [51]. The alignment was done by running single-point calculations using the optimized structures, with ghost atoms added in the near vicinity of the slabs' outermost layers to more accurately describe the electrostatic potential in vacuum. The valence band maximum and conduction band minimum were shifted relative to the vacuum level to obtain a better description of the materials' absolute band alignment. The vacuum level was defined as the asymptotic value of the plane-averaged electrostatic potential sufficiently far away from the slab (≈50 Å for our calculations).