Electrochemical Potential of the Metal Organic Framework MIL-101(Fe) as Cathode Material in Li-Ion Batteries

We discuss the characteristic factors that determine the electrochemical potentials in a metal-organic framework used as cathode for Li-ion batteries via density functional theory-based simulations. Our focus is on MIL-101(Fe) cathode material. Our study gives insight into the role of local atomic environment and structural deformations in generating electrochemical potential.


Introduction
Metal-organic frameworks (MOFs), first synthesized in 1999 [1], are composed of metal-based nodes with organic linkers. They provide porous crystalline structures available in many compositions and geometries. These unique structural properties make MOFs suitable for hydrogen storage, catalytic reactions, gas separation, clean air, and other applications [2,3]. MOFs are also promising for manipulating positronium atoms [4,5], and are of interest in connection with superconductivity and topological properties [6]. Moreover, MOFs have been used as battery materials, particularly for Li-ion batteries [7], as they can reduce the cost for grid-scale battery production [8,9]. MOFs are also attractive because they can be synthesized at lower temperatures from relatively inexpensive precursors [10].
Considering the advantages of MOFs and the low cost of iron, Combelles and Doublet [11] studied the iron-based Fe[OH(BDC)] MOF as a cathode material for Li-ion batteries and reported it as the first example of MOFs with good cycling life and reversible Li insertion. Shirley Shen's group proposed a MOF based on Fe called Matériaux de l' Institut Lavoisier 101 [MIL-101(Fe)] as a cathode material for Li-ion batteries, in 2015 [12]. Shin et al. [12] found that the Fe 3+ /Fe 2+ redox pair is active in MIL-101(Fe), as is the case in the well-known olivine LiFePO 4 material [13][14][15][16]. However, reversibility of this MOF degraded rapidly with cycling, so that appropriate functionalization of MIL-101(Fe) will be needed to address this problem [12]. In 2017, Yamada et al. [17] re-investigated the cathodic properties of MIL-101(Fe) for Li-ion batteries, and reported that the performance can be optimized by changing the mixing ratio of carbon or adhesive materials to incorporate a larger number of Li ions in the material and achieve good cycle ability. More recently, many promising MOFs have been proposed both for anodes and cathodes [18]. Moreover, these typically insulating materials have been made electrically conductive [19].
In this paper, we discuss the behavior of Fe 3+ /Fe 2+ redox orbitals in MIL-101(Fe), by focusing our density functional theory (DFT) calculations on its secondary building unit (SBU). The electron distribution of the SBU, and its electrochemical and magnetic properties have been reported in other studies [20][21][22]. Since the secondary building units (SBUs) are quite well separated and known to be the active sites in MIL-101 [23], it is reasonable to replace the extended periodic structure of the solid with a non-periodic cluster model that can capture the binding properties near the active site. Such an approach has been used in previous studies of MIL-101 materials [23][24][25][26][27] and other Fe-containing organometallics [20][21][22]. Beyond the van der Waals materials, cluster models have been successfully applied for the optimization of reduction activity of covalently bonded Fe-N4@graphene for fuel-cell fabrication [28][29][30]. In our work, this technique successfully unravels the mechanism of redox reactions in this MOF and opens a pathway for identifying irreversible Li intercalation reactions.

Computational Method
Electronic structure calculations in MOFs have been recently reviewed by Mancuso et al. [31]. Coulomb correlation effects are expected to be important, especially at the transition metal atom sites [32,33].
The electronic structure calculations were performed on a simplified cluster model of the large MIL-101(Fe) unit cell adapted from the family of trimetallic [M 3 O(OOCR) 6 L 3 ] + n paddlewheel complexes (M=Fe, R=H, n = 1; no L ligands), representing the SBU of the MOF shown in Figure 1. This model with/without the axial anions/molecules has been successfully applied in many earlier MIL-101 studies [23][24][25][26][27]. Notably, in order to understand various magnetic orders present in the MOF, the system could be treated also with a periodic approach using a very large supercell containing several cluster models, as shown by Trepte et al. [34] in the case of DUT-8(Ni). Our calculations were performed through two main approaches: conventional quantum chemical modeling (geometry optimization and frequency calculation) with Gaussian 16 revision A.03 [35], and single-point energy calculations with RESPECT (relativistic spectroscopy) code [36]. In the conventional calculations, the model structure was optimized using the unrestricted wave function formalisms of the B3LYP, 1 M06-2X, 2 M06-L, 3 PBEPBE, 4 PW91PW91, 5 and ωB97XD 6 methods, which are the methods commonly used in modeling of various MOFs, particularly MIL-101. 78910 11 1213 14151617181920 21222324 252627 (each separated set of references refers to one of the methods; PBEPBE and PBE0; just the first wB97xD ref is for MIL-101). These methods were combined with the def2-TZVP 2829 basis set to increase the accuracy of the calculations while keeping the computational cost low. 30 Because of the importance of van der Waals interactions in MOF systems, the Grimme's D3 (GD3) semi-empirical dispersion correction 31 was used for the B3LYP, M06-2X, M06-L and PBEPBE calculations. For PW91PW91, the required parameters were missing. In the case of ωB97xD, no dispersion correction was needed because the functional inherently includes long-range exchange-correlation correction. 323334 After geometry optimization, the harmonic vibrational frequencies were checked, and the lack of any imaginary frequency approved the stability of the obtained structures. Accordingly, the M06-L and M06-2X calculations were retarded after observing some imaginary frequencies that could not be eliminated using conventional techniques. Furthermore, geometry and spin symmetry breaking were respectively assessed by checking the geometrical parameters and the total spin values before (S 2 ) and after (S 2 A) annihilation of the highest spin contaminant. Finally, to calculate the electrochemical potential as a Li-battery cathode, the interaction between one to three Li + ions with the model clusters (which results in the reduction of the Fe 3+ metal centers to Fe 2+ ) 35 was simulated by inserting the Li + cations in adjacency of the metal centers and optimizing the resultant structures.
For the ReSpect calculations, which provide high accuracy and efficiently for studying materials with In our Gaussian calculations, the model structure was optimized using the unrestricted wave-function formalisms of Becke, 3-parameter, Lee-Yang-Parr (B3LYP) [37], Minnesota M06-2X and M06-L [38], Perdew-Burke-Ernzerhof (PBE) [39], Perdew-Wang (PW91) [40], and Head-Gordon (ωB97XD) [41], which are the commonly used exchange-correlation functionals used in modeling various MOFs, particularly MIL-101 [24,26,[42][43][44][45]. These DFT approximations were combined with the def2-TZVP basis set [46,47] to increase the accuracy of the calculations, while keeping the computational cost low [48]. Because of the importance of van der Waals interactions in MOF systems, the Grimme's D3 (GD3) semi-empirical dispersion correction [49] was used for the B3LYP, M06-2X, M06-L and PBE calculations. For PW91, the required parameters describing the van der Waals interactions are not available. In the case of ωB97XD, no dispersion correction was needed, because the functional inherently includes long-range exchange-correlation correction. After geometry optimization, the harmonic vibrational frequencies were checked. The absence of imaginary frequencies was used to assess the stability of the calculated structures. Accordingly, the M06-L and M06-2X structures were discarded since they produced some imaginary frequencies that could not be eliminated using conventional techniques. Furthermore, geometry and spin symmetry breaking were assessed by checking the geometrical parame-ters and the total spin values, respectively. Finally, in order to calculate the energy gained by using the model cluster as a Li-battery cathode material, the interaction between one to three Li + ions with the model clusters, which results in the reduction of the Fe 3+ metal centers to Fe 2+ [12], was simulated by inserting the Li atoms near the metal centers and optimizing the structures.
For the RESPECT calculations, which provide high accuracy and efficiency for studying materials with heavier atoms, such as transition metals, based on the four-component fully relativistic Dirac-Kohn-Sham approach [36], the single-point energy calculations were performed using the very reliable B3LYP functional only. Because the def2-TZVP basis set is not available in RESPECT, the uncontracted cc-pVDZ basis set was used instead. As we will see in the next section, the basis set difference does not affect the results significantly. We should mention that Gaussian calculations use the VWN3 parametrization of the local density approximation part of the B3LYP functional, while ReSpect results are based on VWN5 parametrization [50].
A detail to keep in mind is that the energy results of RESPECT exclude zero-point energy (ZPE) correction, while the Gaussian results include it.

Spin State of the Pristine Phase
Before undertaking the redox calculations, we checked the ground state of our simplified MIL-101(Fe) model. For this purpose, we examined several spin states using the robust B3LYP functional. Since the lowest and highest spin states are respectively doublet and 16tet, we considered the doublet, quartet, hextet, octet, dectet, 12tet, 14tet and 16tet spin states. However, the geometry optimization of the quartet and hextet structures failed to converge. This is not surprising because these spin configurations involve excited states, where they are not expected to be described well within our DFT-based computations. Table 1 summarizes the results. The optimized structures are shown in Figure S1 of the Supplemental Material (SM). The doublet, octet, 12tet and 14tet states were found to break geometric symmetry of the structure, whereas a realistic model should conserve the C3h symmetry, at least to some extent. The doublet structure is also associated with an imaginary frequency and spin contamination. We assessed the extent of spin contamination by considering the deviation of the eigenvalue of the total spin-squared operator before (<S 2 >) and after (<S 2 A>) annihilation from the expected s(s + 1) eigenvalue, where total spin s = 1 2 n and n is the number of unpaired electrons. Usually, up to 10 % spin contamination is acceptable [51,52]. Table 1. The electronic energy (E), Gibbs free energy (G) at 298.15 K and 1 atm, imaginary frequencies and spin eigenvalues (before (S 2 ) and after (S 2 A) annihilation) of the model at the B3LYP/def2-TZVP level. The energy values in parentheses are the single-point electronic energies calculated by ReSpect at the B3LYP/cc-pVDZ level on the 16tet model (i.e., B3LYP/cc-pVDZ//B3LYP/def2-TZVP). The geometries are shown in Figure S1 of the SM. Abbreviations are as follows. Not available: N/A. Energy convergence failure in ReSpect: (ECF). Geometry optimization convergence failure: GOCF. Geometry symmetry breaking: GSB. Spin contamination: SC. Imaginary frequency: IF. All our calculations predict that the 16tet and 14tet states are the ground and first excited states, respectively, in agreement with the results of Shin et al. [12]. In Gaussian, the electronic energy difference between these two states is 0.43 eV at the B3LYP/def2-TZVP level. The corresponding value for RESPECT at various levels of approximation including relativistic effects is given in Table 2.  Table 3 shows how the energy of the MIL-101(Fe) model varies with the addition of electrons to the system in the ideal case where the lattice relaxation is neglected. This model helps to separate the contributions of structural relaxation and electronic degrees of freedom [14]. If we add one electron to the system, we obtain a gain in energy ∆E ideal = −6.89 eV, as shown in Table 4. When a second electron is added, more energy (2.33 eV) is gained. However, no further energy is gained by the addition of a third electron. In fact, ground state of the model with charge −1 does not contain any empty orbital with negative energy. The calculations were performed at the fixed 16tet geometry with the RESPECT code using the B3LYP exchange-correlation functional and uncontracted cc-pVDZ basis set.

Iron Reduction upon Lithium Intercalation
The results obtained in the ideal case can be extended to model more realistic cathodic conditions by adding one Li + ion per additional electron to the system, and then relaxing both the atomic positions and the electronic structure. The gained energy was calculated using the equation Li,atom (1) in which m is the number of intercalated Li atoms, and E MIL-101(Fe)/mLi , E MIL-101(Fe) and E Li,atom are ZPE-corrected electronic energies of the Li-intercalated MOF, free MOF and free Li atom, respectively. The energy given by Equation (1) can be converted to the voltage with respect to the Li reference anode via the Faraday constant (F) and replacing Li atom energy with the energy of Li metal as shown by Zhou et al. [53] by using the cohesive energy of Li metal E coh = E Li,atom -E Li,metal = 1.65 V [54]: The intercalation geometries and energies are reported in Figure 2 and Figure S2 of the SM. Figure 2 shows that at the B3LYP/def2-TZVP level, the addition of the first electron (i.e., intercalation of the first Li atom) gives a voltage of 2.79 V while the value obtained in the ideal case (without geometrical deformation) is 5.24 V. The reduction peaks in the experimental cyclic voltammogram reported by Shin et al. [12] for MIL-101(Fe) with axial Cl anion (MOF/carbon black weight ratio of 3:7) are at 2.99, 2.59, 2.42, 2.27 and 2.13 V, with respect to the Li/Li + reference, while the corresponding peaks reported by Yamada et al. [17] for water-bound MIL-101(Fe) (mixed with Ketjen black and PTFE adhesive) are at 2.94, 3.08 and 3.50 V. The value obtained from our more realistic model (2.79 V) is in reasonable agreement with the experimental values but lower than the ideal value of 5.24 V. The difference between the ideal voltage of 5.24 V and the model value of 2.79 V indicates that part of the energy gained in the reduction reaction is lost in deforming the MIL-101(Fe) SBU in the presence of Li atom. A similar deformation effect has been found by Hafiz et al. [14] in their modeling of the FeO 6 octahedron in the lithium iron phosphate cathode. The negligible difference between the experimental and theoretical values can be attributed to the sensitivity of the results to axially bound anions and molecules [12] and the matrix composition [17], in addition to the uncertainties inherent in computations (see Figure S2). For the intercalation of the second Li atom, both the ideal and non-ideal cases give the total gained energy of 9.22 eV (i.e., ∆E ideal = ∆E non-ideal = −9.22 eV), showing no net loss of energy, due to the irreversible restructuring of the MOF model. In this case, the energy used for restructuring helps reach a more stable configuration with the release of energy. Intercalation of the third Li atom, however, releases 13.86 eV energy (∆E ideal = −6.93, but ∆E non-ideal = −13.86 eV). Clearly, the energy gained is mostly the result of the structural disintegration of the fragment rather than an electron transfer process. In fact, intercalation of three Li atoms can form a Li-O-C ring-like structure on top of the SBUs, resulting in an irreversible destruction and decomposition of the cathode material when the Fe 3+ atoms reduce beyond a certain point (e.g., 0.6 Li atoms per unit MIL-53(Fe) [55]). This is consistent with the X-ray absorption spectroscopy (XAS) and cyclic voltammetry results of Shin et al. [12], showing a decay in the Li insertion capacity caused by irreversible accumulation of lithium. Significant structural distortion around Fe upon Li + ion-intercalation/Fe-reduction has also been reported in the periodic unit cell modeling study of Combelles et al. on the MIL-53(Fe) MOF [56]. Similarly, other exchange-correlation schemes shown in Figure S2 of the SM yield an energy gain lower than the ideal case upon the addition of the first Li atom, and disintegration of the model structure with the addition of the third Li. However, according to the experimental study of Yamada et al. [17], at some MIL-101(Fe)/carbon mixing ratios, the disintegration does not occur. A possible healing mechanism similar to the one proposed by Barbiellini and Platzman [57] could exist in this case. Li intercalation is accompanied with some changes in the partial charge of the atoms (especially Fe) and variations in the energy of the highest occupied molecular orbital (E HOMO ) and the energy of the lowest unoccupied molecular orbital (E LUMO ), the Fermi level (E F = E HOMO + (E LUMO − E HOMO )/2 and the band gap (E g ) energies. The changes obtained from natural bond orbital (NBO) analysis are summarized in Table S1 in the SM. Table S1 shows that the intercalation reaction induces a redistribution of charge in the system: the decrease of charge on Fe atoms reflects reduction of Fe +3 to Fe +2 . Notably, absolute partial charges are not equivalent to formal oxidation states [58]. However, as the density of states (DOS) plots in Figure S3 show, the changes in orbital energies and atomic charges are consistent with the changes seen in the DOS plots. Furthermore, the differential charge densities (Fukui functions [59]) shown in Figure S4 help to visualize the redox orbitals [14] in the present case. In addition, the NBO results related to the intercalation of one Li atom at the B3LYP/def2-TZVP level show that the intercalation of the Li atom mainly relies on three charge transfer processes from the lone pair (LP, Lewis type) orbital of the bridging O atoms to the anti-bonding LP (LP * , non-Lewis) orbital of Li

Conclusions
We have presented a comprehensive theoretical study that provides insight into the mechanism of lithiation/delithiation and the potential variation in a MOF cathode. DFT simulations including structural distortions offer a quantitative understanding of the changes in total energy induced by the modification of the structure due to electron orbitals of transition metals and ligands. Our first principles approach enables us to gain a molecular-level understanding of the relationship between structural distortion, potential changes and redox orbitals in battery MOFs. Our method can thus elucidate structuralactivity relationships of working MOF-based battery materials to improve their performance.
Supplementary Materials: The following are available online at www.mdpi.com/xxx/s1, Figure S1