Introduction of model systemsWe here consider COF model systems that are experimentally accessible and that are constructed from well-known precursors, which are summarized in Fig. 1. We study a diverse set of structures, including triphenylene (TP), phthalocyanine (Pc) and fluorinated Pc derivatives with imide linkages (for TPQ− and PcQ−) or diamine linking groups (for TPQ+, PcQ+, and PcF8Q+). These molecules are frequently used in COF synthesis45,46,47. For our purpose it is highly relevant that they exhibit a large variation in their quadrupole moments (QPMs), represented by their in-plane QPM component (of the traceless QPM tensor) as indicated in Fig. 1. The sign of these QPMs is also used as superscript (either Q+ or Q−) for compound names in this work for clarity. Some of these building blocks, namely Pc and fluorinated Pcs, have been employed successfully in crystalline thin films and blends to engineer electrostatic potentials48,49. Their characteristic potential distribution in gas phase is overlaid with the structure in Fig. 1. For better illustration, the potential is evaluated at a distance of 2.5 Å above the molecular plane to avoid the strong atomic core potentials at smaller distances that would otherwise dominate the images.Fig. 1: Monomer building blocks.Molecular structures of TPQ−(a), TPQ+(b), PcQ−(c), PcQ+(d), and PcF8Q+ (e) with their DFT electrostatic potentials plotted 2.5 Å above the molecular plane for better visibility (see details in the “Methods” section). Positive values of the potential (red) refer to electron attraction. All values are given with respect to the reference potential, which is the vacuum potential whose value is set to zero. In-plane quadrupole moments are given below each figure, representing both in-plane directions.For the sake of clarity, we first discuss the results of the three TP-based COFs and generalize the concept of pore potential engineering by QPMs to other structures later. In these TP-based COFs, the TP units act as vertices (nodes), which are connected by benzene rings. Figure 2 shows the structures of the COF monolayers as relaxed by density functional theory (DFT) [cf. “Methods” section]. All three materials are constrained to a hexagonal lattice and the same symmetry group P6/mmm. However, they were chosen to have different linkages between the TP units and the benzene linker, which can be realized by the different monomeric units used. The structures in Fig. 2 are overlaid with the electrostatic potential from DFT. We choose the convention that positive potential values indicate attraction of electrons while negative values represent repulsion, as compared to the vacum level (zero potential, see “Methods” section). We observe that the potential within the pore is strongly negative if TP units are linked by an imide linkage group (Fig. 2a), but strongly positive if they are linked by two amine groups (Fig. 2b). We therefore denote these structures COF-TPQ− and COF-TPQ+, respectively. The potential of COF-TPQ− closely resembles the potential of TPQ− (Fig. 1a) and, analogously, the potentials of COF-TPQ+ and TPQ+ are also similar. Therefore, the sign-flipped molecular quadrupole moment between both molecules is successfully inherited to the potential of the corresponding COFs and may therefore qualify as tuning knob in COF engineering. The third COF, denoted COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\), is chosen to be a variant of COF-TPQ+ that is derived from it by formally abstracting two hydrogen atoms per linkage and thereby introducing formal conjugation on the entire structure. As a result, the potential in the pore turns negative, seemingly at odds with the potential of the TPQ+ monomer, indicating a somewhat more involved relationship between monomer QPM and COF potential.Fig. 2: Structure of TP COF monolayers and DFT electrostatic potential.The potential is calculated at a distance of 2.5 Å above the structure. Positive values of the potential (red) refer to electron attraction. All values are given with respect to the reference potential, which is the vacuum potential whose value is set to zero.In order to predict the COF’s pore potential from the monomer properties, we develop a simple ab initio-based model. Indeed, a model that avoids expensive DFT calculations of bulk COFs but uses simple DFT parameters of the constituents for predicting the potential in the bulk, would be highly desirable and could be used in the COF design process or to formulate guiding principles. Towards this goal, we rationalize the potential at the center of the pore by considering the charge distribution over the framework. Given that the electrons are located at the framework structure and that the electron density is exponentially suppressed toward the center of a pore, a multipole expansion for the potential at the TP’s center-of-masses is permitted. Without actually performing the expansion, its form can easily be anticipated and the leading terms will be discussed. Since the COFs are charge neutral, there are no monopole contributions in this expansion. Also, the dipole contribution vanishes due to the inversion symmetry of the moieties (see further discussion of other scenarios at the end of this paper). Thus, the QPM is the first non-vanishing contribution in the expansion. Therefore, the minimal model for the pore potential we strive for, consists only of quadrupole moments arranged periodically in the lattice. Possible higher multipole moments in the expansion can be neglected for the model. The resulting electrostatic potential Vmodel(r) is then obtained as a sum of quadrupole fields,$${V}_{{\rm{model}}}({\boldsymbol{r}})=\frac{1}{8\pi {\epsilon }_{0}{\epsilon }_{r}}\sum _{k}\frac{{({\boldsymbol{r}}-{{\boldsymbol{r}}}^{(k)})}^{T}\cdot {\hat{{\boldsymbol{Q}}}}_{{\rm{model}}}^{(k)}\cdot ({\boldsymbol{r}}-{{\boldsymbol{r}}}^{(k)})}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}^{(k)}{| }^{5}},$$
(1)
where \({\hat{{\boldsymbol{Q}}}}_{{\rm{model}}}^{(k)}\) is the traceless quadrupole tensor of the k-th monomer at position r(k). The high symmetry of the COFs allows for additional simplifications. In the present case, the in-plane components \({Q}_{{\rm{xx}}}^{(k)}={Q}_{{\rm{yy}}}^{(k)}\) are equal due to symmetry and \({Q}_{{\rm{zz}}}^{(k)}=-2{Q}_{{\rm{xx}}}^{(k)}\) to ensure that the entire tensor is traceless. In addition, the off-diagonal elements are zero. This strongly simplified tensor structure occurs for all sites and essentially reduces the complexity of our model to a single parameter. For concreteness, we introduce Qmodel for the xx-component of \({\hat{{\boldsymbol{Q}}}}_{{\rm{model}}}^{(k)}\). This means that for a given lattice, this single parameter characterizes the electrostatic potential, while its value varies with the chemical structure of the COF and its vertex.We first verify the model potential at the center of the pore by comparing Vmodel(r) with the ab initio potential from DFT calculations. For a simple and practical model, we find that it suffices to include only the quadrupole moments that are closest to the pore’s center (illustrated by blue dots in Fig. 3a) in the sum of Eq. (1), while an extension to more distant QPM centers would be straightforward. The resulting potential is evaluated at the center of the pore for various out of plane distances d. A comparison of the quadrupole model with the ab initio electrostatic potential is shown in Fig. 3b and reveals an excellent agreement for all three TP COFs over a wide range of distances. We therefore conclude that the potential inside a pore can be well described as a superposition of quadrupole fields.Fig. 3: Modeling of the potential.a The potential is evaluated at the center of the pore (blue cross) and for certain distances d to the framework (along the blue path starting from the center at d = 0). Vmodel is calculated according to Eq. (1) with the included QPM centers represented by the blue dots. b The DFT potential is compared to Vmodel for COF-TPQ−, COF-TPQ+, and COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\). All values are given with respect to the vacuum potential.The fitted values of the QPMs for the model potential are summarized in Table 1. The sign of the fitted QPM reflects the sign of the DFT potential, i.e., a positive in-plane QPM of COF-TPQ+ leads to a positive potential, whereas COF-TPQ− and COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\) show negative values for Qmodel. However, when comparing the fitted Qmodel with Qmono, i.e. the xx-component of the tensor \({\hat{{\boldsymbol{Q}}}}_{{\rm{mono}}}\) obtained from DFT for the monomers (cf. Table 1), they coincide only for COF-TPQ+ (within 5%), indicating that the current model based on the TP building blocks is not complete. In particular the differences between COF-TPQ+ and COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\) sharing the same TPQ+ unit cannot yet be captured. The choice of both these COFs for our study showcases the importance of the linkages in engineering the pore potential.Table 1 Quadrupole moments and model parameters for investigated COFs and their constituentsThe impact of the COF linkages on the potential is traced back to the different charge densities in this region (away from the TP units). Figure 4a shows the difference between charge densities \(\Delta \rho ={\rho }_{{{\rm{TP}}}^{Q+}}-{\rho }_{{{\rm{TP}}}_{{\rm{conj}}.}^{Q+}}\) at the same geometry apart from differing H atoms. The yellow isosurface represents regions with more electrons in the non-conjugated COF, whereas the blue features indicate a higher density for the conjugated COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\). This analysis reveals that most of the change in the electronic structure occurs around the linkage, while Δρ is small at the vertices. This indicates that taking into account this effect of Δρ, i.e., considering also the linkages within the COF, should improve the model.Fig. 4: Modeling of the potential with COF linkage.a Electron density difference of COF-TPQ+ and COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\). Yellow isosurfaces indicate lower electron density in the conjugated structure, the blue regions exhibit higher electron density in the conjugated structure. b Illustration of the two parameters \({\hat{{\boldsymbol{Q}}}}_{{\rm{mono}}}\) (quadrupole of the monomeric TP units) and the linker quadrupole contribution \(\delta \hat{{\boldsymbol{Q}}}\) defining the potential (see main text). Both are located by definition at the vertex (solid circle). The linker contribution arises from point charges δq located at the linkages with distance r to the vertex (indicated by arrows) and reflect the change in charge density upon linkage of the monomeric building units.These differences in the charge density between the two TPQ+ COFs can be simply included when considering them as resulting from auxiliary charges δq near the linkage. We introduce these effective charges at a distance r from the TP vertices as defined in Fig. 4b (open circles). These charges replace the more complex DFT-derived charges of Fig. 4a and depend only on the linker, while r is simply defined by the closest distance to a vertex and may flexibly account for the specific COF geometry. By the same symmetry arguments, these charges δq give rise to another quadrupole tensor \(\delta {\hat{{\boldsymbol{Q}}}}^{(k)}\) for each vertex. Grouping three of them (indicated by arrows in Fig. 4b) allows choosing their center at the vertex (filled circle), independently of whether they are caused by the linker. As a result, the total quadrupole tensor is the sum of both contributions$${\hat{{\boldsymbol{Q}}}}_{{\rm{model}}}^{(k)}={\hat{{\boldsymbol{Q}}}}_{{\rm{mono}}}^{(k)}+\delta {\hat{{\boldsymbol{Q}}}}^{(k)},$$
(2)
where \({\hat{{\boldsymbol{Q}}}}_{{\rm{mono}}}^{(k)}\) are the tensors of the monomers and \(\delta {\hat{{\boldsymbol{Q}}}}^{(k)}\) are determined by auxiliary charges δq near the linkages. The latter two can be simply related by using the symmetry of TP-COFs according to$$\delta {\hat{{\boldsymbol{Q}}}}^{(k)}=\frac{3}{2}{r}^{2}\delta q\left(\begin{array}{ccc}1&0&0\\ 0&1&0\\ 0&0&-2\end{array}\right).$$
(3)
Since the involved distance r is given by the geometry of the COF, δq is the only free parameter that needs to be determined from DFT for each type of linkage. These auxiliary charges δq mostly depend on the chemical structure and the nature of the bonds in the immediate vicinity of the linkage, e.g., functional groups near the linkage. In the vast majority of COFs, this does not affect the overall charge of COF building blocks. Exceptional cases, however, are conceivable in which a stronger electron transfer between building blocks leads to their charging like in charge-transfer salts. In such cases, the model requires extensions to include these scenarios with long-ranged charge transfer. In absence of this, δq is determined for a certain linkage and could be applied to other COFs with that same linkage without further DFT calculations.We finally use Eq. (2) (resolved for \(\delta \hat{{\boldsymbol{Q}}}\)) and fix the auxiliary charge δq. They are placed directly at the site of the linkage between TP vertex and benzene linker, i.e. in between the two nitrogens for COF-TPQ− (cf. Fig. 4b) and on the nitrogen atom for COF-TPQ+. The resulting values are compiled in Table 1. This finalizes the model for the electrostatic potential of a COF’s pore that only depends on two parameters Qmono and δq, where Qmono describes the contribution from the monomers and δq the contributions from the linkage. This means that the potential is traced back to its basic building blocks. Once the parameters for those building blocks are known, one can easily estimate the electrostatic potential for new COFs without further ab initio calculations, which is a huge simplification.Generalization for other COFsWe first validate our model by extending our study to a family of Pc-based COFs with different symmetry, which are shown in Fig. 5 together with their electrostatic potential. The Pc monomeric building units (cf. Fig. 1) are connected via benzene and the resulting COFs exhibit four-fold symmetry47,50. As a validation test, we determine, as previously, \(\delta \hat{{\boldsymbol{Q}}}\) from these COFs and the quadrupole moments of the building blocks, thereby deducing the δq anew. All parameters are summarized in Table 2, confirming the robustness of these values with minor modifications in the Pc COFs. A graphical overview over the values of these bond parameters for the entire set of COFs is provided in Fig. 6a by plotting δq of each COF against the averaged value for the specific type of linkage. It shows a very good agreement within each set for both the imide and diamine-linked COFs, i.e., δq is characteristic for the bond linkage because these specific values vary much less between COFs with the same linkage than between the different linkages. For example, the values of δq for the imide linkage of COF-TPQ− and COF-PcQ− (−0.027 e and −0.025 e) are in good agreement. They are clearly distinguished from the diamine-linked COFs, whose average δq value for COF-TPQ+, COF-PcQ+, and COF-PcF8Q+ is much smaller (δq = 0.0035 e ± 0.01 e). This clustering of δq for the same type of linkage, is also found for the conjugated COFs COF-\({{\rm{TP}}}_{{\rm{conj}}.}^{Q+}\), COF-\({{\rm{Pc}}}_{{\rm{conj}}.}^{Q+}\) and COF-\({\rm{PcF}}{8}_{{\rm{conj}}.}^{Q+}\) with (δq = −0.13 e ± 0.03 e). The larger scatter here may be attributed to the strongly electro-negative fluorine atoms close to the site of the linkage. The relative scatter of 23% indicates a second-order effect. That is δq is, besides being dominated by the linkage, also slightly influenced by the neighboring PcF8 unit. This is also visible in the difference of the charge density (Supplementary Fig. 1). If not considered, as we do for simplicity, this scatter can be taken as an error estimate for the parameter. The minor influence on the potential, however, is of greater importance and will be shown below.Fig. 5: Structure of Pc-based COF monolayers and DFT electrostatic potential.The potential is calculated at a distance of 2.5 Å to the plane of phthalocyanine (Pc) and fluorinated Pc COFs with and without formal conjugation at the linker unit. All values are given with respect to the reference potential, which is the vacuum potential whose value is set to zero.Table 2 Summary of quadrupole moments and model parameters.Fig. 6: Linkage parameters and model verification.a δq obtained for every COF is compared to the average value for respective type of linkage. b Electrostatic potentials calculated by DFT are plotted against Vmodel at distances 0, 5, and 10 Å to the layer. The color represents the monomer building unit, the symbol denotes the type of linkage, and the size of the symbol indicates its distance from the layer.In order to verify the model’s validity for the potential of COF monolayers, Fig. 6b compares the DFT potential to Vmodel for the investigated systems at various distances d to the layer. The modeled potential is in excellent agreement with the ab initio results throughout all COFs and distances. In any case the influence of the bond linkage (represented by δq and therefore \(\delta \hat{{\boldsymbol{Q}}}\)) is strong and their consideration in the model leads to Vmodel ≅ VDFT and, particularly, to the corrected sign of all potentials.Multilayer COFs and bulk systemsBesides single monolayers, our interest is in the tuning effect of multilayers and bulk systems as these systems dominate the COF literature. Most of the experimental COFs consist of stacked 2D layers in which the pores form large tubes51,52, and the potential inside such tubes and at its entrance is of high relevance. Owing to the additive property of electrostatic potentials, the potential inside a COF tube should be the superposition of the potentials of each layer and hence describable by our model. This is because the π-stacking only marginally changes the charge density of a single layer. To verify this, we compare the pore’s potential of the three Pc-based systems to a quadrupole field as described above and add two layers to represent a three-layer stack. The ab initio potential of the COFs is calculated with eclipsed (AA) stacking and agrees very well with the model potential as demonstrated in Supplementary Fig. 2. The obtained quadrupole moments are similar to the monolayer ones and only the dielectric screening due to adjacent layers reduces the strength of the potential, which is not the case for monolayers. Based on the superposition of QPMs when stacking multiple layers and taking the obtained model parameters including the dielectric screening, we extend our study to larger stacks and predict the potential for bulk materials.Figure 7 shows the potential of Pc-based systems when approaching the bulk limit, for which we include up to 55 layers (20 nm) and 150 × 150 QPM centers within each x–y plane (~3.3 μm in each direction) to construct Vmodel. Each individual layer is marked by either a white or a gray stripe to provide a clearer understanding of the proportions. Comparing Fig. 7a and b for AA stacking, we observe that the potential within COF-PcQ− and COF-PcQ+ can be successfully tuned to extreme opposite values. Herein, the sign and magnitude of the potential reflects their different QPMs. It should be emphasized that the energy difference between the absolute minimum and maximum in Fig. 7 can exceed 2 eV for thicker materials. Interestingly, in both systems the potential approaches a plateau with a constant value over more than 15 layers inside the COF. In addition, a steeply increasing/decreasing potential (within a few layers) is present in both cases at the entrance/exit of the tubes, independently of the magnitude of the respective QPMs. Further analysis shows that the potential’s shape only depends on the pore size and stacking distance of the COF, which are very similar in both cases and only enter the model via the geometry factors r(k) (cf. Eq. (1)).Fig. 7: Prediction of bulk potentials.a, b Potential for bulk systems with increasing number of layers, made from a COF-PcQ− and b COF-PcQ+. Number of layers are indicated in the legend. The potentials are described by the quadruple model. All potentials approach zero for infinite distance to the samples (not shown for illustration purpose). Qxx was fitted here for a three-layer reference bulk to be −23.9 DÅ (COF-PcQ−) and 45.6 DÅ (COF-PcQ+). c, d Analogous simulations for COF-PcQ− in different stacking geometries as indicated in the insets with c AA’ serrated and d inclined stacking fashion. Gray and white stripes both indicate position of individual layers.In order to investigate the influence of stacking, we compare exemplarily for COF-PcQ− the eclipsed stacking in Fig. 7a with serrated AA’ and inclined stacking fashions in Fig. 7c and d, respectively. We use a lateral shift of 1.5 Å in all cases. For AA’ stacking, the path for plotting the potential is chosen to go through the averaged pore centers as indicated in the inset. Neither qualitative nor quantitative changes are found for the potential when compared to the AA stacking. In case of inclined stacking, the path passes through the center of the pores of each layer (cf. inset of Fig. 7d). The tilt of the tube axis against the principle axis of the QPM tensor leads to small decrease of the effect by ~10%, since in-plane and out-of-plane components partly cancel along this direction. Also small variations are visible close to the surface. In general only minor changes are observed. This can be rationalized by the large unit cell of COFs as compared to the lateral displacement that is typically observed between the layers, which is a result of the van der Waals bonding. We therefore expect that a random lateral displacement will also marginally change the results.It is finally very instructive to examine how the potential wells observed in few-layer samples gradually transition into plateaus, with their potential values reaching saturation at a specific thickness. At first glance, this appears counter-intuitive, as one might expect that adding more layers (including more QPM centers and their potentials) would result in a logarithmic increase of the overall potential that would not be limited. An increase is indeed observed for systems with less than 15 layers. However, with increasing bulk size, in-plane and out-of-plane components of the quadrupole tensors counteract each other due to their different signs. To illustrate this, we reformulate Eq. (1) in terms of in-plane and out-of-plane components, expressing the spatial difference r − r(k) through their Cartesian components Δx(k), Δy(k) and Δz(k). Using further that \({\hat{{\boldsymbol{Q}}}}_{{\rm{model}}}^{(k)}\) only depends on a single parameter Qxx, which is equal for every k (see above), we express the potential as$${V}_{{\rm{model}}}({\boldsymbol{r}})=\frac{{Q}_{xx}}{8\pi {\epsilon }_{0}{\epsilon }_{r}}\sum _{k}\frac{{\left(\Delta {x}^{(k)}\right)}^{2}+{\left(\Delta {y}^{(k)}\right)}^{2}-2{\left(\Delta {z}^{(k)}\right)}^{2}}{| {\boldsymbol{r}}-{{\boldsymbol{r}}}^{(k)}{| }^{5}}.$$
(4)
This form highlights the action of the in-plane components \({\left(\Delta {x}^{(k)}\right)}^{2}+{\left(\Delta {y}^{(k)}\right)}^{2}=:{\left(\Delta {\rho }^{(k)}\right)}^{2}\) and the out-of-plane components \(-2{\left(\Delta {z}^{(k)}\right)}^{2}\) of every site k. They always have different signs because the underlying QPMs are traceless. Hence, the numerator describes a double cone with 109.5° apex angle as illustrated in Fig. 8. QPMs inside the cone decrease the potential because \({\left(\Delta {\rho }^{(k)}\right)}^{2}-2{\left(\Delta {z}^{(k)}\right)}^{2} < 0\), i.e. the negative out-of-plane component dominates. Conversely, QPMs located outside the cone increase the potential, i.e. in-plane components dominate. Figure 8 further shows schematically COF systems of different bulk sizes (layer numbers) with the QPM centers indicated by blue and orange dots for a dominating in-plane and out-of-plane QPM, respectively.Fig. 8: Interplay of quadrupole components.The COF system is schematically illustrated by QPMs distributed in layers. The 109.5° apex angle of the double cone determines whether their out-of-plane (orange) or in-plane component (blue) dominates the potential at the cone apex. Compared are situations with few layers (a), a bulk system (b) and a COF surface (c).This double cone structure can be used to rationalize the different regimes. In monolayer and few-layer cases, almost all QPMs lie outside the cone (indicated by blue dots) and therefore contribute via their in-plane components. These components add together with the same sign, increasing the magnitude of the potential. This no longer applies to many layers as they approach bulk systems, as illustrated in Fig. 8b. When more layers are added, the layers near the middle primarily contribute with their in-plane component, while only a few QPM centers lie within the cone. With increasing vertical distance to the cone apex, more QPM centers lie within the cone (orange dots) and therefore the potential in the center of the COF’s pore decreases by adding those layers. As a result of an increasing number of layers, the electrostatic potential decreases at the center and eventually assumes a plateau because the relevant environment of QPMs remains unchanged for shifts of the cone apex.Interestingly, the potential changes its sign with varying d and shows a barrier at the exit of the tubes before it approaches zero for infinite distances to the sample. The barrier increases with the number of layers and converges to a finite value for the case of infinite many layers. We can illustrate this transition by shifting the double cone apex in Fig. 8c. At a certain distance to the surface, most QPMs are located within the cone giving rise to an inverted potential compared to the bulk. The magnitude of the potential values inside the pore and the barrier heights depend on the geometry, i.e. pore size, symmetry, and stacking distance. As the pore size of COFs is typically much larger than the layer distance, it is expected that the potential shape found in Fig. 7 will be prototypical for COF structures.
Engineering the electrostatic potential in a COF’s pore by selecting quadrupolar building blocks and linkages
