Computational screening study towards redox-active metal-organic frameworks

The metal-organic framework (MOF) MFU-4 l containing Co(II) centers and Cl− ligands has recently shown promising redox activity. Aiming for further improved MOF catalysts for oxidation processes employing molecular oxygen we present a density-functional theory (DFT) based computational screening approach to identify promising metal center and ligand combinations within the MFU-4 l structural family. Using the O2 binding energy as a descriptor for the redox property, we show that relative energetic trends in this descriptor can reliably be obtained at the hybrid functional DFT level and using small cluster (scorpionate-type complex) models. Within this efficient computational protocol we screen a range of metal center/ligand combinations and identify several candidate systems that offer more exothermic O2 binding than the original Co/Cl-based MFU-4 l framework.


Introduction
Next to applications in gas storage or drug delivery, metal-organic frameworks (MOFs) receive increasing attention as redox active catalysts [1][2][3][4][5][6][7][8][9]. Combining bi-or multifunctional ligands (predominantly polycarboxylate ligands) and (transition) metal ions, moderately robust MOFs can be prepared. Advantages of MOFs as heterogeneous catalysts are enhanced catalyst stability due to the spatial separation of single catalytic sites in the framework [10], high porosities in the absence of any non-accessible bulk volume (dead volume) [6] and-most evidently-their pore size(s), and thus their substrate shape and size selectivity, that can be systematically tailored by employing different organic linkers [11]. Three conceptually different approaches have been used to gain catalytically active MOFs: introducing catalytic metal centres at the nodes of the constituting secondary building units (SBU), attaching catalytically active functional groups or coordination units at the MOF struts (i.e. the organic linkers), and embedding nanosized metal clusters within the pores of the MOF. Until now, the use of catalytically active metal centres positioned at the nodes of the MOF-constituting SBUs was predominantly explored for Lewis acid catalysis [12][13][14][15], whereas examples on redox catalysis are rare.
In terms of catalytic transformations especially MFU-4-type frameworks show several highly promising features [16]. Particularly large pore sizes were achieved by the use of bistriazolo-dibenzo-dioxin (BTDD) linkers, which then lead to the so-called MFU-4 l members within this class of MOFs (figure 1(a)) [17]. Their Kuratowski-type pentanuclear SBUs [M a(II) M b 4 × 4 (L 6 )] (figure 1(b)) offer four well accessible, tetrahedral metal centers M b each being coordinated by three N donor ligands. Most appealing for redox catalysis is that each such coordination site is stereochemically and topologically related to the typical coordination motif found in so-called scorpionates (metal complexes with hydrotris(pyrazolyl)borato ligands, see figure 1(c) and figure 1 in the supplementary data, available at stacks.iop.org/ NJP/15/115004/mmedia), a famous class of coordination compounds that was invented and developed by S Trofimenko from the early 1960s onwards [18]. Dioxygen activation has already been reported for cobalt scorpionates [19], as has the feasibility of oxygen insertion into aliphatic C-H bonds [20]. While cobalt complexes seem thus particularly promising for catalytic oxidation processes [21], the search for other catalytically active metal oxo complexes that can activate and insert oxygen into C-H bonds remains a flourishing area of active research [22]. Since at present there is no workable concept for incorporating scorpionate complexes into porous frameworks, MFU-4-type frameworks comprising this coordination motif could fill this gap, from which finally a novel class of technically and commercially attractive heterogeneous catalysts may evolve. Recently, a first step along this route was achieved by the demonstrated redox activity of a Co-based MFU-4 l framework [23]. Furthermore, this study revealed an important additional advantage of MFU-4 l frameworks in terms of the possibility for a postsynthetic metal center substitution, i.e. in principle there is a generic route to also explore other redox active metal centers. Notwithstanding, the iterative and time-consuming development of efficient syntheses leading to phase-pure MOFs, as well as their structural characterization imposes severe limitations on the number of compounds which can be examined in catalytic test reactions. This is particularly vexing, when considering that the non-radical catalytic oxidation or oxygenation activity will also critically depend on the interplay between the metal center and the electronic nature of the ligand X, cf figure 1, i.e. there is a quite large matrix of MOFs that needs to be explored.
In this situation, computational screening can step in as the most valuable tool to identify interesting ligand-metal center combinations on which ensuing synthesis activities would then be focused. An intuitive and well computable descriptor for the redox properties is in this context simply the O 2 binding energy to the metal center, reflecting the degree of bond activation. Such a screening study is only of use though, if it is much faster than experimental synthesis, yet still sufficiently predictive. The prior demand necessitates the use of numerically efficient approaches. While we will show below that highest efficiency can be reliably achieved with respect to system size by focusing the calculations on the small scorpionate-type moieties, the situation is less clear with respect to the level of theory necessary to reach predictive quality. The intricate electronic structure of metal-oxo complexes is generally known to severely challenge the present-day workhorse in terms of efficient electronic structure calculations, density-functional theory (DFT) with semi-local or hybrid functionals [24,25]. This is most clearly exemplified for the present case by the qualitative differences reported recently for the oxygen binding energy in a similar Co-based MOF within the MFU-1 structural family [26]. At the gradient-corrected level this binding was obtained as exothermic, while at the hybrid level it was endothermic with the absolute difference between the two computed binding energies exceeding 1 eV! In this work we address this problem by tracing it largely back to the self-interaction induced deficiency of the semi-local functional in properly describing the change of metal spin state concomitant to O 2 adsorption. As this is rather independent of the actual ligand, screening calculations at corresponding levels of DFT theory are still meaningful as they only require relative energetic trends to identify suitable MOF candidates, not the absolute binding energy itself. Demonstrating the additional possibility to perform these calculations on the small scorpionate-type units, this brings us into a position to efficiently screen a range of metal centers and ligands. The results clearly identify a small number of interesting combinations, among which particularly Mn and Fe centers with NH −1 2 ligands additionally prevent possible limitations due to spin-transition induced O 2 adsorption barriers.

Theory
All DFT calculations have been performed with the all-electron full-potential code FHIaims [27,28]. Its localized basis sets based on numeric atom-centered orbitals (NAO) allow one to perform the calculations both for finite molecular systems and infinite systems through the use of periodic boundary conditions (in the present case with -point sampling). Electronic exchange and correlation was treated on the level of the generalized gradient approximation (GGA) PBE functional [29] or on the level of the hybrid B3LYP functional [30,31]. Dispersive interactions lacking at these levels of theory were considered through the dispersion-correction scheme due to Tkatchenko and Scheffler [32], and we correspondingly denote this as PBE + TS and B3LYP + TS below. The central quantity used as descriptor in the computational screening approach is the O 2 binding energy defined as , and E tot (O 2 ) are the total energies computed for O 2 adsorbed at the MOF, for the clean MOF, and for the isolated O 2 gas-phase molecule, respectively. The calculation of this numerically inexpensive quantity is sufficient for the relative energetic trends required here, and could ultimately be augmented by vibrational calculations to allow comparison with experimentally measured reaction enthalpies. In the employed sign convention negative binding energies indicate exothermicity. Geometry optimization was generally performed at the PBE + TS level until residual forces fell below 10 −4 eV Å −1 . Hybrid B3LYP + TS energetics was subsequently obtained through single-point calculations on these optimized geometries. Systematic convergence tests show E b (O2) to be sufficiently converged for the intended relative energetic trend study (i.e. within ±50 meV) at the 'tier2' NAO basis set and when using tight integration grids. Spin and charge assignments in the text are based on calculated Hirshfeld projections.

Definition of screening space: ligand and metal centers
The starting point for the present study is the reported redox activity of the peripheral Co(II) centers contained within the coordination units of a MFU-4 l [M a Co b 4 Cl 4 (ta) 6 ] (M a = Zn or Co, tatriazolate ligand) framework [23]. Preserving the structural character of the MFU-4 l class of MOFs, obvious routes to be explored towards a further improved activity would be metal centers other than Co and ligands other than Cl − . In terms of ligand substitution we hereby concentrate on ligands that keep the metal center in the same formal charge state as for Cl − , namely OH − , These anions have been chosen because they are representative of ligands comprising different sigma and pi-metal bonding capabilities, i.e. they cover the whole range of different ligand types described by the so-called spectroscopic series of ligands found in classical monographs of ligand-field theory (vide infra). Among these OH − , Cl − , NO − 3 and CF 3 SO − 3 are considered as weak-field ligands providing different denticities (OH − , Cl − , monodentate; NO − 3 , bidentate; CF 3 SO − 3 , bi-or tridentate), whereas F − , NO − 2 , and CN − were chosen owing to their well-known propensity to interact via additional metal-ligand π -bonds and thus are regarded as strong-field ligands. NH − 2 and H − , finally, serve as strongly basic and reducing ligands. As to the metal centers significant reactivity with C-H bonds has particularly been observed for metal-oxo complexes that are significantly destabilized by occupied d-orbitals (i.e. possessing d n configurations with n = 4 or higher) [22]. Only then is the metal-oxo multiple bond sufficiently weakened to render alternatives, such as the formation of hydroxo groups, energetically competitive. Aiming to produce M = O units from M-O 2 precursors, we correspondingly focus within the 3d transition metal series particularly on Mn, Fe, Co and Ni centers.
For the metal-dioxygen interaction we can generally expect the existence of a weak physisorptive adsorption mode. In such a van der Waals bonded complex the dioxygen molecule will preserve its gas phase spin triplet state. In addition to this, direct coordination to a single metal center can occur in a side-on η 2 -mode or an end-on η 1 mode. We systematically tested both chemisorptive adsorption modes for the Co-based scorpionates described below. In full agreement with previous DFT work on the similar Co-based MFU-1 framework [26], the side-on η 2 -mode in which the peripheral Co(II) centers of MFU-4 l effectively switch to an octahedral configuration results as by far most favourable for all monodentate ligands (>0.4 eV, except for the NH −1 2 ligand where the difference is only 0.15 eV). For bidentate ligands (NO − 2 , NO − 3 , CF 3 SO − 3 ) this O 2 η 2 -binding mode implies a switching to monodentate ligand coordination, but we found this to be still more favourable or at least energetically en par to the alternative end-on η 1 O 2 coordination with bidentate ligand. All screening calculations are correspondingly restricted to the side-on η 2 chemisorptive mode.

Reduction of system size
The primitive trigonal cell derived from the original face-centered cubic unit cell of MFU-4 l with a simple Cl − ligand comprises 162 atoms. While the localized basis sets within FHI-aims allow to efficiently treat such infinite systems within a periodic boundary condition setup, the computational costs are still noticeable, in particular at the hybrid functional level. Suspecting the O 2 binding to the tetrahedral metal centers to be predominantly determined locally, we therefore additionally consider two types of finite cluster models built from the MOF SBU. On the one hand this is the Kuratowski complex shown in figure 1(b), which is obtained by truncating the linkers of one entire SBU after the phenyl moiety and saturating the broken bonds with hydrogen (89 atoms). On the other hand this is the complex shown in figure 1(c) (and in a different perspective in figure 2), which only contains a single tetrahedral active center and three adjacent pyrazole units (27 atoms). This closely mimics the coordination environment of one  of the tetrahedral metal centers of the Kuratowski complex, but could similarly be viewed as a metal-modified scorpionate complex, in which the HB unit of a typical scorpionate is replaxed by a Zn(II) metal ion. For simplicity we will refer to this complex as scorpionate complex in the following, even though it is formally rather a 'metal-modified' scorpionate. Table 1 summarizes the calculated PBE + TS O 2 binding energy at the two cluster models for a Co(II) center and a series of ligands. For the O −1 2 side-on binding mode already the difference between the two cluster models is only of the order of ∼0.1 eV, and the local bonding geometries around the metal center are essentially identical. Equivalent calculations for physisorbed O 2 , i.e. the van der Waals bonded complex, show even smaller energetic differences of the order of 40 meV. For Cl − and CN − ligands both adsorption modes have also been computed in the fully periodic MFU-4 l framework, revealing insignificant energetic and geometric differences to the results obtained for the Kuratowski complex (<20 meV). Finally, O −1 2 side-on binding energies at scorpionate and Kuratowski complexes were also calculated at the B3LYP + TS level, and exhibited differences of the same order of magnitude as those found at the PBE + TS level. We correspondingly conclude that the oxygen binding at the smallest scorpionate complex mimics the one in extended MFU-4 l quite faithfully, with the indicated uncertainty of ∼0.1 eV in the binding energies not relevant for the intended screening approach. All calculations reported from now on are therefore conducted most efficiently for this small cluster model.

Energetic trends at gradient-corrected and hybrid DFT level
As mentioned in the introduction the conflicting results for gradient-corrected and hybrid DFT reported by Tonigold et al for the highly similar MFU-1 system [26] question the reliability of   either level of theory for the description of the O 2 -MOF interaction. Before embarking on the actual screening study we therefore further analyse this difference for the present system. Within the perspective of a screening approach that primarily intends to identify energetic trends over a series of possible ligands or metal centers, we thereby concentrate on Co based scorpionates and compute the O 2 binding energy for a range of ligands. The results summarized in table 2 reveal a similar disagreement between gradient-corrected and hybrid DFT as observed for the MFU-1 system [26], with PBE + TS yielding a significantly stronger binding. Intriguingly, however, this difference of about 1.1 eV is rather independent of the various ligands, which suggests that its origin lies primarily in the treatment of the Co metal center itself. Thinking along this line we recall that O 2 binding in the side-on mode effectively changes the Co coordination from tetrahedral to octahedral. Figure 2 exemplifies the corresponding optimized geometries for the Cl − ligand, with equivalent results obtained for all other ligands tested. In agreement with the expectations from ligand-field theory this coordination change is accompanied with a change in the Co spin state, cf figure 3: in the tetrahedral coordination of the plain MOF both functionals predict a quartet high-spin state for the Co(II) center, whereas a singlet low-spin state is predicted for the octahedrally coordinated Co(III) after dioxygen binding. Suspecting the concomitant d-orbital rehybridization to largely contribute to the discrepancy between the two levels of DFT theory, we perform calculations in which we already force the plain MOF structure to the low-spin state, namely a doublet state for Co(II). As apparent from figure 2 this immediately breaks the tetrahedral symmetry, with the Cl − ligand moving towards the position it also exhibits after dioxygen adsorption. Table 3 compiles the energetic difference between this low-spin state and the true highspin state as obtained by PBE + TS and B3LYP + TS. Interestingly, we again find a rather large Our analysis thus traces a large part of the discomforting discrepancy between the two DFT functionals back to the different ability of the gradient-corrected and hybrid functional in describing the high-spin state of the tetrahedrally coordinated Co(II) center in the plain MFU-4 l structure. Specifically, we suspect the generally acknowledged electron delocalization problem of the gradient-corrected functional [33] to lead to a wrong account of the fully-occupied d e orbitals. This interpretation gets support from the observation that in the calculated PBE + TS density of states (DOS) these orbitals are partly found at higher energies than the d t2 orbitals, in contrast to the predictions from ligand-field theory, cf figure 3. In contrast, in the B3LYP + TS DOS this is no longer the case, i.e. here all d e orbitals are consistently located below the d t2 orbitals. In fact, we can further corroborate the central relevance of this d orbital ordering by artificially enforcing the correct ordering at the PBE + TS level through the application of an appropriate U operator [34]. In corresponding PBE + U + TS calculations, in which the U value of 1.2 eV has been set to fix the d e orbitals at the position they have in the B3LYP + TS calculations, we then indeed obtain O 2 binding energies that are very close to those obtained with the latter functional (within 0.1 eV).
We therefore conclude that the discrepancy between the gradient-corrected and hybrid functional calculations are predominantly due to the wrong d orbital positioning predicted by the prior electron-delocalization riddled functional. As this is mostly independent of the actual ligand employed, already calculations at the computationally most efficient PBE + TS level should yield reliable energetic trends within the space of screened ligands. On the other hand, we cannot expect the same to hold for different metal centers, i.e. electron delocalization induced errors in the d orbital positions will generally be different for different metals. This view is confirmed by evaluating the difference in computed O 2 binding energies at PBE + TS and B3LYP + TS level for scorpionates with Mn, Fe and Ni centers. Consistent with the results presented for Co based scorpionates in table 2 we each time find the energetic difference between the two levels of theory to be rather independent of the actual ligands. Notwithstanding, Table 4. DFT B3LYP + TS O 2 binding energies (eV) at scorpionates with different metal centers and ligands. Interesting candidates with exothermic binding energies below −0.3 eV are marked in bold. this difference varies largely for the different metal centers. Whereas it was ∼1.1 eV for the above discussed case of Co, it is ∼0.55 eV, ∼0.25 eV and ∼0.9 eV for Mn, Fe and Ni, respectively. Such variations preclude a comparison of the O 2 binding energy descriptor for different metals at the PBE + TS level. We correspondingly perform the actual screening calculations over the full range of metals and ligands at the B3LYP + TS level. At this level the electron delocalization error is largely removed through the admixture of exact exchange [35] and we expect reliable energetic trends also across the different metal centers tested.

Screening the ligand/metal center matrix
Within the established computational screening approach (scorpionate cluster model, DFT B3LYP + TS) we proceed to compute the O 2 binding energy in both the van der Waals physisorptive mode and the side-on chemisorptive mode for a range of ligands and different metal centers. As expected, physisorption is possible at all tested metal centers. The bond strength varies little with the actual ligand and metal center, and lies between −0.1 and −0.2 eV throughout. Suspecting that such a weak binding does not allow for sufficient bond activation, we thus concentrate on the chemisorptive mode and aim to identify metal center-ligand combinations that yield more exothermic binding energies than those offered by physisorption. yields generally a very weak binding, and concomitantly low bond activation. From the calculations in the hypothetical low-spin state for the plain scorpionate we attribute this to a rather high energy cost for the spin change and rehybridization to octahedral coordination in the presence of this ligand. This rehybridization cost is generally also the reason for the quite large number of in fact endohedral binding energies computed: The actual O 2 -MOF interaction is certainly stronger than in the physisorptive mode and we validated that in all cases the chemisorptive binding mode is metastable, i.e. corresponds to a minimum on the potential energy surface. Nevertheless this energetic gain is counterbalanced by the rehybridization cost to octahedral coordination, which in many cases leads to a net endothermic binding energy. Encouragingly, other ligands than the original Cl − seem to allow for more flexibility, cf table 3, and correspondingly show more exothermic O 2 binding. Particularly interesting in this respect are the H − , OH − , and NH − 2 ligands. There is, however, a non-trivial correlation between ligand and metal center. A ligand that shows high exothermicity for one metal center does not necessarily show this for another metal center, too. In other words, the energetic trend across the tested ligands is different for the different metal centers. This dictates the explicit computation of the full ligand-metal center matrix, rather than a selective optimization among the ligands for one metal center and subsequent optimization of the metal center for this optimized ligand. Distributed over the matrix we find several metal center-ligand combinations that yield a binding energy more exothermic than −0.3 eV, and thus promise a bond activation significantly stronger than in the purely physisorptive adsorption mode. Before concluding on these combinations as interesting candidates for experimental synthesis, one has to recognize though that computational screening is only as useful as the actual descriptor employed. The computed binding energy is an intuitive measure of bond activation and ensuing redox activity. Yet, it does not tell about unwanted side reactions that e.g. compromise the structural stability of the entire MOF. In this respect, the use of H − as ligand could be critical, since we can expect a further reaction of bound O 2 molecules leading to hydroperoxide species, which cannot be transformed back to hydride and thus no catalytic activation of O 2 can be achieved.
Furthermore, a computed thermodynamic metastability of the chemisorptive binding mode does not tell about possible kinetic limitations in reaching this state. Here, we particularly note the largely differing spin properties of the Co centers compared to the other metal centers tested. Consistent with the expectations from ligand-field theory we find the change from tetrahedral to octahedral coordination concomitant with O −1 2 side-on binding to be accompanied with a high-spin to low-spin transition in the case of cobalt. Specifically and independent of the actual ligand, the Co(II) center in the plain framework changes from its quartet spin state to a singlet state as Co(III) in the octahedral configuration after O −1 2 binding. Together with the spin transition of the O 2 molecule from its triplet gas-phase state to the doublet state after adsorption, this amounts to a total change in spin multiplicity of S = 4. In contrast, for Mn and Ni we obtain smaller spin multiplicity changes of S = 2, and for Fe a S = 0, since Mn and Fe metal centers stay in their high-spin state after adsorption. In ligand field theory such a high multiplicity change in the case of Co necessarily imposes a pairing energy penalty on the reaction, since two electrons formerly in singly occupied orbitals need to be paired into one orbital, cf figure 3. If this happens before the energy gain due to the O 2 chemisorption sets in a barrier results along the reaction path even in an electronically adiabatic picture, where additional spin selection rules are not even considered [36][37][38]. As illustrated in figure 4 for a NH −1 2 ligand this is indeed what we find in our adiabatic DFT calculations, i.e. we obtain a sizable activation barrier of the order of 0.6 eV to reach the chemisorptive state. Equivalent calculations for other ligands confirm that such a high activation barrier is a specific property of the Co centers, whereas in particular for the promising combination Fe or Mn center with NH −1 2 ligand we obtain vanishingly small activation barriers. In light of these findings these two systems appear finally as most promising candidates for experimental synthesis among those included in our screening study: both yield exothermic binding energies in the side-on mode that are significantly stronger than pure physisorption, cf table 4, and this adsorption mode can be reached without kinetic limitations.

Conclusions
We have performed a computational screening study to identify promising alternative metal center-ligand combinations for the MFU-4 l metal-organic framework that could improve on the redox activity achieved recently with Co metal centers and Cl − ligands [23]. Using the O 2 binding energy as descriptor to assess this redox activity, we first established that calculations at the DFT hybrid functional level and using small Scorpionate-type cluster models represent a numerically efficient and reliable computational protocol for the required relative energetic trends. Screening Mn, Fe, Co, and Ni centers together with a range of ligands we identify several candidates that afford a chemisorptive O − 2 side-on binding mode that is notably more stable than mere van der Waals physisorption. Allowing to reach this chemisorptive state without significant activation barriers, particularly Fe or Mn centers with NH −1 2 ligands appear as most promising for an ensuing experimental synthesis. A major limiting factor towards even larger O 2 activation is the high energy cost connected with the rehybridization from tetrahedral to octahedral coordination concomitant to side-on O −1 2 adsorption. This suggests to extend the screening approach either to metal(II) centers with higher reactivity towards oxygen (e.g. V II , Cr II ) or to ligands that additionally stabilize the coordinated dipolar dioxygen e.g. through intramolecular hydrogen addition.