Electronic structure of porphyrin-based metal-organic frameworks and their suitability for solar fuel production photocatalysis

Metal-organic frameworks (MOFs) can be exceptionally good catalytic materials thanks to the presence of active metal centres and a porous structure that is advantageous for molecular adsorption and confinement. We present here a first-principles investigation of the electronic structure of a family of MOFs based on porphyrins connected through phenyl-carboxyl ligands and AlOH species, in order to assess their suitability for the photocatalysis of fuel production reactions using sunlight.


INTRODUCTION
The development of cheap, efficient techniques to carry out the photocatalytic splitting of water would permit the generation of a clean fuel (hydrogen) at low negative environmental impact. [1][2][3] Of similar interest is the photocatalytic reduction of CO 2 , which would allow the synthetic production of carbon-containing fuels (e.g. methanol) 4,5 and simultaneously contribute to recycle CO 2 from the environment. It is clear that the development of efficient technologies to carry out these energetically up-hill reactions using solar energy would be greatly beneficial, and therefore, many research efforts are being devoted to it, in particular to the search for adequate photocatalysts. [6][7][8] This research was initially focused on traditional inorganic semiconductors, such as TiO 2 and CdS 3,[9][10][11][12] , but it has now extended to a wider class of materials, including nanostructures such as fullerenes, nanotubes and graphene-like 2D solids. [13][14][15][16][17] On the other hand, water splitting and CO 2 reduction are at the core of natural photosynthetic reactions, 18 so the study of the related natural processes can help finding artificial routes for these reactions. 3,[19][20][21] Bioinspired molecular photocatalysts have been largely studied in recent years. [22][23][24] Besides Mn-complexes that are close to natural photosynthesis reactions, 23,25,26 porphyrins have also been identified as active molecular centres for artificial photosynthesis. 27,28 One drawback of molecular systems from a practical point of view is their recyclability, as the separation of the catalysts from liquid media is very difficult. To overcome this problem, one attractive route is their immobilization in solid hosts. 29,30 Metal-organic frameworks (MOFs) have appeared as promising hosts, [31][32][33][34] where catalytic centres can be encapsulated [35][36][37][38][39] or moreover be part of the constituents of the materials. [39][40][41][42] After pioneering work by Suslick and co-workers, 43 a number of porphyrin-based MOFs have been reported, 44-47 including some with photocatalytic properties. [44][45][46] The choice of MOFs as catalysts is also stimulated by their diverse porous architectures, 48,49 which confer them exceptional molecular adsorption properties. [49][50][51] In these systems the adsorbed reactant molecules can access the active sites embedded in a confinement field that favours the catalytic reactions 52 . Furthermore, it is possible to tune the adsorption properties of MOFs, by changing factors such as topology, metal composition or the nature of the ligand. 53,54 However, the impact of these modifications on the electronic structure of MOFs has not been widely studied. Gascon et al. 55 showed that by modifying the linker properties the overall bandgap of MOF-5 can be lowered. The dependence of the optical properties of a nanotubular MOF on the adsorbed guest molecules can be exploited in molecular sensing applications. 56 The magnetic and optical properties can also be tuned by the composition of the metal centres. 57,58 The catalytic behaviour of nanoporous solids strongly depends on both the structural and electronic features of the materials, so a rational design of these materials would make them very useful in catalytic applications. In contrast to semiconductor photocatalysts, in terms of their electronic behaviour, MOFs can be regarded as molecular-like catalysts. 59 Fateeva and co-workers 46 have showed that the optical and photocatalytic properties of a MOF containing porphyrins connected through phenyl-carboxyl ligands and AlOH species (called Al-PMOF here) are primarily determined by the porphyrin linker. For Al-PMOF these authors observed a strong absorption band at 415 nm (2.99 eV) and four lower energy bands, which are characteristic of the free porphyrin molecule in solution. The photocatalytic properties of porphyrin molecules can in principle be modified by the presence of metal (M) cations within the ring. 60 For example, metalation of Al-PMOF with Zn has been investigated to improve the efficiency of photocatalytic water splitting, 46 while CO 2 conversion is improved by Cu incorporation. 61 It is therefore relevant to study the electronic structure of metallated Al-PMOF materials in order to understand the effect of the porphyrin metal on photocatalytic reactions. In this paper we consider the incorporation of Fe, Co, Ni, Cu or Zn in the centre of the porphyrins in Al-PMOF, in comparison also to the di-protonated Al-PMOF. For this purpose, we use Density Functional Theory (DFT) to calculate the electronic bandgap, and the absolute positions of the band edges, as these are important factors in determining the suitability and efficiency of a material as photocatalysts.

COMPUTATIONAL METHODS
In our calculations the structures were represented by a primitive rhombohedral cell, which contains half the number of atoms as the conventional orthorhombic cell shown in Figure 1 (space group Cmmm) 46  We performed spin-polarized quantum-mechanical calculations using density functional theory (DFT) as implemented in the Vienna Ab-initio Simulation Package (VASP). [62][63][64][65] Geometry optimizations were performed using the generalized gradient approximation (GGA) with the Perdew-Burke-Ernzerhof functional (PBE) functional. 66,67 During relaxation, forces on atoms were minimized until they were all less than 0.01 eV/Å.  ) with one unpaired electron, whereas Zn 2+ is not spin-polarised. In each calculation, we constrained the difference in the number of up-and down-spin electrons to the corresponding integer number given above. Because there is only one metal atom per cell, all our spin-polarised calculations correspond to ferromagnetic configurations. Considering the relatively large distance between magnetic metal ions (ca. 6.7 Å in the direction perpendicular to the porphyrin plane) and the absence of common ligands, it can be expected that the strength of the magnetic coupling will be very small and that the actual structure is paramagnetically disordered at room temperature. However, the weak magnetic coupling implies that there must be very small energy differences associated with magnetic ordering in these structures, and therefore the simulation results should be independent on the magnetic order assumed. In order to obtain accurate electronic structures we carried out single-point calculations at the most favourable spin state for each composition, using the screened hybrid functional of Heyd, Scuseria and Ernzerhof (HSE06), 68,69 which generally provides bandgaps in closer agreement with experiment than those from GGA functionals. 70 Although the HSE06 calculations were single-point only, based on the PBE-optimised structures, test calculations showed that the final electronic structures were largely insensitive to the small geometric variations introduced by re-optimising the structures at more sophisticated levels of calculation, e.g. adding corrections to improve the description of d orbitals and of van der Waals interactions (see Supplementary Information). The projector augmented wave (PAW) method 71,72 was used to describe the frozen core electrons and their interaction with the valence electrons, i.e. those in level 4d for Fe, Co, Ni, Cu and Zn, and 2s2p for C, N and O. The kinetic energy cutoff for the plane-wave basis set expansion was set at 600 eV. A Γ-centred grid of k-points was used for integrations in the reciprocal space, where the smallest allowed spacing between k-points was set at 0.5 Å -1 , giving rise to 3 irreducible points in the Brillouin zone corresponding to the primitive cell. As in other periodic DFT codes, the band energies in VASP are given with respect to an internal energy reference (the average potential in the crystal). In order to align the band energies with the vacuum scale, it is necessary to evaluate the electrostatic potential in the vacuum region represented by an empty space within the simulation cell. In the present study we follow the methodology recently proposed by Butler et al. 73 to calculate the vacuum level in MOF structures, which consists of evaluating the average potential within a small sphere (radius of 2 Å) at different positions in the pore. By finding the point that is farthest apart from the framework atoms (the pore "centre"), where the potential is locally flat (no electric field), we can obtain a good approximation to the vacuum level. In Ref. 73, this procedure led to MOF ionization potentials in good agreement with experiment. A Python code provided by these authors was employed in our calculations to obtain the average potentials. 74

Crystal structures
We first examine the crystal structures resulting from the substitutions of different cations in the porphyrin. Upon replacement of the two H atoms by the transition metal atoms, only small variations of cell parameters are observed (actions between 0.1% and 0.6%.  The metal-nitrogen interaction seems to be the main factor controlling the small changes of the cell volume. The calculations systematically overestimate all the cell parameters in comparison with experiment in Ref. 46 by an average of ~1%. However, the calculations are able to reproduce well the small variations observed experimentally in the cell parameters from the protonated to the Zn-substituted structure: a small contraction of the a parameter, a small expansion of the c parameter, while the b parameter remains roughly the same. Overall, a small contraction of the cell volume by is observed in both the experiment (0.31%) and the calculations (0.13%) upon Zn substitution.

Thermodynamics of metal substitutions from aqueous solution
In order to discuss the stability of the material with the different metal cation (M 2+ ) substitutions, we have calculated the enthalpies and free energies for the process of cation exchange with aqueous solution, according to the following reaction: using a mixed theoretical-empirical approach for the treatment of the aqueous cations. For each cation, the enthalpy in aqueous solution is approximated using the following quantities: a) the DFT energy of the neutral atom in the gas phase, calculated with VASP at the experimental spin groundstate, using a large supercell and the same functional and precision parameters as in the MOF calculations; b) the sum of the experimental first and second ionization energies, 75,76 which is the energy needed to ionize the neutral atom to a M 2+ cation (of course, only the first ionization energy is used in the case of hydrogen); and c) the experimental hydration enthalpy, 77 which is the enthalpy change in the process of moving the M 2+ cation from the gas phase to aqueous solution. The addition of these three contributions gives the enthalpy of the aqueous cation, which will be used in the calculation of the enthalpy change of reaction (1), ∆ . On the other hand, the reaction free energy ∆G can be estimated as: where ∆ 0 is calculated using the same procedure as for ∆ , but employing the experimental hydration free energies of the cation (instead of the hydration enthalpy), i.e., considering the entropy contribution to hydration. 78 The second term, where k B is Boltzmann's constant, takes into account the effect of the relative concentrations of the cations in aqueous solution. The obtained values for ∆ and ∆ 0 are reported in Table 2. As shown in Figure 2, the variation with cation concentration and pH is relatively weak. Increasing the solution pH (i.e. decreasing the proton concentration) or increasing the metal concentration in solution, makes the exchange reaction slightly more favourable. The calculated free energies for reaction (1)    agreement with the observed experimentally easy introduction of Cu 2+ in Al-PMOF. 61 The preference for particular spin states, which can be seen from Table 2, will be discussed below in terms of the electronic structure.

Electronic structure
The total electronic density of states (DOS) and their projections on the 3d orbitals of the metals are shown in Figure 3. All structures are semiconductors, with bandgaps in the range between 2.0 and 2.6 eV. A bandgap of around 2 eV is generally considered to be ideal for single-semiconductor water splitting photocatalysis. 79 Anatase TiO 2 , one of the most widely investigated photocatalysts for water splitting, has a In order to understand the position in the DOS of the metal 3d contributions we need to refer to their splitting due to the porphyrin ligand field. Although the position of porphyrin as a ligand in the spectrochemical series is considered to be ambiguous, it is known that in metallo-porphyrins with square planar coordination of divalent cations, porphyrin tends to act as a strong field ligand. 83 That means that the d x 2 -y 2 level is very high up in energy compared with the others, as shown schematically in Figure 4. But in this MOF, the degeneracy of the two lower d levels is broken, because the crystal does not have the 4-fold rotation axis that would be present in the isolated gas-phase porphyrin. The distortion could be then referred to as rectangular distortion, to indicate the lowering of the symmetry to a 2-fold rotation axis.  ). There is only one possible spin state for this cation (μ=1), since the four low energy levels are filled and one electron occupies the d x 2 -y 2 level. e) Zn 2+ (d 10 ). In this case, the five d orbitals are filled, so μ=0.
Regarding the bandgap, we can see in Figure 3 that all materials (except Fe-Al-PMOF) have similar values of bandgap, between 2.3 eV and 2.6 eV. The projection of the DOS on individual atoms and orbitals (not shown in the figure) reveals that the HOMO is associated with the p orbitals of the N and C atoms of the porphyrin moieties, while the LUMO is associated with the p orbitals of the C and O atoms of the carboxyl ligands. In the case of Fe-Al-PMOF, the Fe atoms introduce the d xy levels into this otherwise unoccupied region (close to the valence band), leading to a decrease in bandgap to 2.02 eV. Note that the DOS plots in Figure 3 are aligned with respect to the Fermi level, which means that the introduction of a peak in the band gap region shifts all the peaks to the left in the case of Fe. The bandgaps of the materials clearly make them good candidates to carry out photocatalytic reactions, since they would be able to absorb most of the solar radiation. However, in order to assess the ability of a material to perform photocatalytic water splitting or carbon dioxide reduction, we also need to investigate if the semiconductors exhibit the correct alignment of the bands with respect to the half-reaction potentials, which we do in the following section.

Band edge positions with respect to electrode potentials
The alignment of the bands with the vacuum level allows us to explore the thermodynamic feasibility of the photocatalytic processes. A single-semiconductor photocatalyst requires certain important characteristics in its electronic structure. For example, for the water splitting reaction, the positions of the conduction and valence band edges should straddle the redox potentials for water photolysis, 79,85,86 i.e. the valence band edge should be below the energy of the oxygen evolution reaction (OER): and the conduction band edge should be above the energy corresponding to the hydrogen evolution reaction (HER): The energy scale has the opposite sign of the potential scale, so lower energy means higher potential, and vice versa. The bandgap must therefore be wider than 1.23 eV (difference between the HER and the OER levels). After loss mechanisms are accounted for, a bandgap of 2 eV or more is generally considered as necessary, 79 but the bandgap should not be too wide, in order to allow the adsorption of photons from the visible part of solar radiation. It is known that, in the vacuum scale and at pH=0, the HER level is located at -4.44 eV, and the OER level is located at -5.67 eV. 87 At temperature T and pH > 0, these energy levels are shifted up by pH×(k B T×ln10). By referencing the electronic levels in our semiconductor solids with respect to the vacuum level (taken here as the electron potential at the centre of the largest pore), we can assess whether the band edges of the semiconductor are in a favourable position to catalyse the solar splitting of water under a given set of conditions. In the case of carbon dioxide conversion to fuels (e.g. methane, CH 4 , methanol, CH 3 OH, or formic acid, HCO 2 H), the position of the conduction band of the semiconductor photocatalyst has to be above the redox potential for the CO 2 reduction half-reaction, which depends on the specific fuel produced. Since the CO 2 /CH 4 , CO 2 /CH 3 OH and CO 2 /HCO 2 H levels are above the HER (H + /H 2 ) level, the photocatalyst for the CO 2 reduction reactions requires a minimum bandgap that is wider than for water splitting. For all materials, the conduction band edge is at roughly the same position (ca. -3 eV in the vacuum scale), which is above the energies corresponding to the H + /H 2 and CO 2 /CH x Oy levels. That means that the MOF would be thermodynamically able to donate an excited electron from the conduction band for the reduction half-reactions to proceed. The conduction bands are slightly higher than desired for water splitting, but they are in a nearly ideal position for the carbon dioxide reduction reactions. As for the oxygen evolution half-reaction, in all materials, except Fe-Al-PMOF, the valence band edge lies at about -5.5 eV in the vacuum scale, below the energy of the O 2 /H 2 O level. In the case of Fe-Al-PMOF, it lies slightly (~0.2 eV) above. Of course, due to the approximations made in our calculations, such small energy differences are not reliable. But it is clear that the valence band edge for all compositions is in the correct energy range in the absolute scale. Small deviations from the ideal band edge positions in a semiconductor with the adequate bandgap can always be corrected via the application of a weak bias voltage in a photoelectrochemical device, shifting both band edges with respect to the redox levels (applying a voltage increases device complexity and also consumes energy, so a large bias voltage should be avoided). For example, in a recent first-principles screening of materials for water splitting photocatalysts, Wu et al. set the threshold of the allowed bias voltage to 0.7 V. 89 Our present calculations show that the Al-PMOF, in all the different compositions explored here, would be able to operate as single-semiconductor photocatalyst with little or no bias voltage applied at neutral pH.

Conclusions
We have carried out a theoretical study of the electronic properties of a porphyrin-based metal organic framework, including both protonated and metallated (Fe 2+ , Co 2+ , Ni 2+ , Cu 2+ , Zn 2+ ) porphyrins. We have found that the protons would be spontaneously exchanged with these cations, when the materials are submerged in aqueous solutions of these cations. The analysis of the electronic bands reveals that the bandgaps of all materials are in the favourable range for efficient adsorption of solar light (2.0 to 2.6 eV). Furthermore, the alignment of the bands is also favourable in all cases for the photocatalysis of water splitting and carbon dioxide reduction, which means that a device using these materials would require little or no bias voltage to function. Our calculations only show small variations in the electronic band edges of the metal-substituted structures. The only exception is the Fe-substituted one, where an occupied d xy state is introduced in the gap region above the valence band. In all the other cases studied here, the band edges are determined by the porphyrin electronic structure. This means that the choice of metal at the porphyrin centre could be used to optimise other properties, e.g. molecular adsorption, affecting the photocatalytic process. Our calculations demonstrate that these porphyrin-based MOFs are very promising candidates for efficient photocatalysis of fuel production reactions. Further studies about the adsorption and diffusion of molecules in these structures would be helpful to select the best composition for performing the photocatalytic reactions studied.

Effects of Hubbard and dispersion corrections on the results
In our simulations, we used the GGA-PBE functional to optimise the geometries, while the electronic structures were calculated using single-point runs with the HSE06 functional. We have followed this approach because we are dealing with large simulation cells, where it becomes very computationally expensive to carry out the optimisations at the HSE06 level.
Here we provide evidence to show that the errors in the PBE optimisation are not expected to carry over to the electronic structure calculation. We have selected one of the structures (the one with Ni, which in this case is non spin-polarised, and therefore easier to calculate) for the test calculations.
In Table S1, we compare the geometries obtained with the following methods: a) PBE, b) PBE-D2, where the van der Waals interactions are taken into account via the DFT-D2 method of Grimme, c) PBE+U, where the U eff parameter that describe the on-site Coulomb interaction of d electrons is set at 3 eV, and d) PBE-D2+U (i.e., both corrections applied). Table S1. Calculated lattice parameters, cell volume and the two perpendicular N-N distances inside the porphyrin for the Ni-Al-PMOF structure, using different DFT functionals for the optimisation. As in Table 1  The cell parameters exhibit little deviation from one calculation to the other (maximum discrepancy is ~1% along the b axis). The small effect of the U correction on the cell parameters is due to the fact that it only affects the metal at the centre of the rigid porphyrin ligand. On the other hand, the effect of the D2 correction is also small (although larger than for the U correction) because the whole structure is covalently linked.
The question remains on whether these small geometric changes can affect the final electronic structure. Figure S1 shows the positions of the band edges, as calculated with an HSE06 single point on the geometries optimised with different functionals. Clearly the variations in the positions of the band edges are very small (less than 0.04 eV). Figure S1. Bandgaps and band edge positions of the Ni-Al-PMOF system with respect to the vacuum level, as calculated with the HSE06 functional on geometries optimised at different levels of calculations (PBE, PBE-D2, PBE+U with U=3 eV, and PBE-D2+U). Energy levels corresponding to redox potentials of water splitting and carbon dioxide reduction reactions producing methane, methanol, and formic acid at pH = 7 are also shown with dotted lines.