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-Pack k-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 G = 2k 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 d-subshells were projected onto the total DOS plots showed in the Supplementary Information. Geometry visualization and analyses were performed using VESTA. Absolute band-edge alignments were carried out following the procedure described on the "CRYSTAL tutorials" website. 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).