Screening for high-spin metal organic frameworks (MOFs): density functional theory study on DUT-8(M1,M2) (with Mi = V,,Cu)

We present a first principles study of low-spin (LS)/high-spin (HS) screening for 3d metal centers in the metal organic framework (MOF) DUT-8(Ni). Various density functional theory (DFT) codes have been used to evaluate numerical and DFT related errors. We compare highly accurate all-electron implementations with the widely used plane wave approach. We present electronically and magnetically stable DUT-8(Ni) HS secondary building units (SBUs). In this work we show how to tune the magnetic and electronic properties of the original SBU only by changing the metal centers.


Introduction
New and efficient technologies are needed for the development of smart devices for data/signal-transfer, -manipulation and -storage. Nowadays many devices are based on magnetic materials. The general trend goes to smaller and smaller devices, but the downsizing and miniaturization in magnetic storage devices is physically limited. Below a critical domain size (superparamagnetic limit) thermal excitations can flip the spin orientation and data loss is the result. With this work we show how to realize a promising concept for solving this problem by integrating the so-called single-molecule magnets (SMMs) into the crystalline nature of so-called metal organic frameworks (MOFs). The research field of MOFs is currently under development in many directions, as an example for interesting synthetic work on magnetic MOFs we would like to mention the work by G. Christou et al. 1 SMMs are metal-organic compounds that show magnetic behaviour below a critical temperature (blocking temperature). Magnetism in such compounds is localized and no long-range order occurs. MOFs consist of metal-centers and organic-building blocks, which are repeated due to their crystalline nature. The secondary building units of MOFs may be structurally considered as SMMs. We want to show that a SMM building unit for a MOF combines the advantages of both the local magnetism and the periodic boundary conditions, which would allow to arrange local magnetic sites in a distinct three-dimensional order. DUT-8(Ni) is a so-called flexible MOF, 2-4 which electronic and magnetic structure has already been described. 5 Our goal is to find a high-spin (HS) solution for the magnetic coupling of the metal centers of this MOF. In other words we search for a MOF with stable local HS magnetic ordering which results in a stable ferromagnetic MOF. This may be interesting for a possible application in spintronics or a molecular switchable gas-sensor. In this work we present a study on how the magnetic properties of MOFs behave by changing the 3d metal centers. Density functional theory has limited accuracy due to several systematic shortcomings like e.g. self-interactionerror 6,7 and the problem of spin-contamination 8 (see ESI, † A3.1). However, it is a widely used and well established theory for the calculation of ground state energies. To determine exchange coupling constants, DFT is used for the calculation of the energy difference between low-spin and high-spin states, which should give qualitatively accurate results. Those energy differences are used as parameters for the calculation of the coupling constant J as given in a model Heisenberg Hamiltonian.

DFT calculations
All calculations presented in this work have been carried out in the framework of DFT. 9, 10 We used all-electron based program packages (codes) NRLMOL, [11][12][13][14][15][16][17][18] and ORCA 20 and compared the results with results from plane wave codes QUANTUM ESPRESSO (QE) 21 and GPAW. 22 Most of the calculations were performed with the generalized gradient approximation (GGA) 16 using the PBE 23 exchange-correlation functional. To verify that this functional is sufficient for our calculations, we repeated some calculations using the B3LYP 24,25 exchange-correlation functional as implemented in ORCA. The NRLMOL code uses an optimized Gaussian basis set, 18 numerically precise variational integration and an analytic solution of Poissons equation to accurately determine the self-consistent potentials, secular matrix, total energies and Hellmann-Feynman-Pulay forces. The FPLO code is a full-potential local-orbital minimum-basis code 19,26 to solve the Kohn-Sham equations on a regular lattice. 27 For the QE calculations, the projector augmentedwave (PAW) 28 method with the already mentioned GGA-PBE 23 exchange-correlation functional was used. The calculations were restricted to G point. The kinetic energy cutoff is 90 Ry for all systems. For further details we used the ORCA code with unrestricted Kohn-Sham DFT (UKS) and a triple z basis set with valence polarization function (def2-TZVP), the NRLMOL code with a DFT optimized Gaussian basis set, 18 the GPAW code with linear combination of atomic orbitals (LCAO) and a double z basis set with polarization function (DZP) and the QE code with the PAW method and a plane wave (PW) basis set ([element].pbe-(n)-kjpaw_psl.0. [1,2,3].[0,1,3].UPF, for details see theossrv1.epfl.ch/ Main/Pseudopotentials). All calculations were performed spinpolarized. An inclusion of the van-der-Waals interaction was considered. Test calculations show no effect on the magnetic properties. Thus it is neglected for all systems. In this work we only perform calculations on cluster models, thus finite systems. To be able to use codes implementing periodic boundary conditions (QE, GPAW) on such systems, we created unit cells with a sufficient amount of vacuum. The vacuum was chosen such that any cluster has a minimum distance of 10 Å to a cluster in any adjacent unit cell. This allows the application of QE and GPAW for our molecular systems. The usage of several codes ensures the reduction of methodological errors, which may occur in every implementation of DFT. Another reason is the comparison between all-electron and pseudopotential codes, as it becomes unfeasible to use all-electron codes for larger systems like MOFs. Thus it has to be proven that pseudopotential implementations reproduce all-electron results. This is important for further investigations, where the model systems can be used as SBUs for new MOFs.

Magnetism and exchange coupling constant
Magnetism can be described by the coupling of local spins at distinguishable magnetic centers (e.g. i and j). A possible description is given by the so-called Heisenberg-Dirac-Van Vleck (HDVV) Hamiltonian 29-31 where J ij is the coupling constant between neighboring spins and -S i / -S j are spin operators. A high-spin state (parallel aligned spins, ferromagnetic coupling) is indicated by a positive sign of the coupling constant J ij while a negative sign refers to a low-spin state (anti-parallel spins, anti-ferromagnetic coupling). For dimers (i = 1 and j = 2) the coupling can be expressed with the total energies of these two different magnetic orderings 32,33 For dimers which include the same metal centers or metal centers with the same number of unpaired electrons, eqn (2) reduces to eqn (4), because h -S 2 i LS becomes 0. Thus, eqn (4) is only valid for those cases and not usable for other kinds of mixed dimers. Further approximation can be done by assuming that S HS 2 c S HS , which transforms eqn (4) into eqn (3). This is not suitable for our study, as this assumption is not valid in all systems. With that, all calculated coupling constants were derived from eqn (2). A more detailed discussion about the used calculations scheme, i.e. restricted open shell Kohn-Sham (ROKS for NRLMOL, FPLO and QE) and unrestricted Kohn-Sham (UKS for ORCA) as well as a detailed discussion about the broken symmetry method is given in the ESI † (see A3.1) or see the excellent review paper of Neese. 8

Variation of transition metal centers
Recently we described the electronic and magnetic properties of the flexible metal organic framework DUT-8(Ni) (see Trepte et al. 5 ). Flexible means that this MOF has an open and a closed phase. In this previous work, 5 a detailed discussion of the influence on the coupling constant of different structural parameters was carried out using several model systems. The model system M1 (see Fig. 1 and ESI, † A1) is a good approximation to study the magnetic behaviour of the open crystalline system. In this M1 model system the Ni atoms are bipyramidally coordinated with four O and one N each. The interatomic Ni-Ni

Results
In preliminary calculations we calculated the coupling constant of the original model system M1 containing Ni using various DFT codes with different implementations and levels of precision (see Table 2). We investigated the dependence of the calculated coupling constants on the choice of the exchange-correlation functional by comparing results obtained with PBE and B3LYP. The value of the coupling constant changes, but the sign and the magnitude are retained (see ORCA PBE/B3LYP results Table 2). Thus we decided to use the PBE functional for the HS screening. It should be considered that only in the UKS formalism spincontamination is explicitly calculated (see ESI, † A3.1). For this reason we used eqn (3) for the evaluation of the results from ORCA to make those accurate ORCA values comparable with all other results. This corresponds to the ideal spin operator expectation values in the LS and HS state. All calculations of the coupling constant with the PBE functional give the same sign and a very similar absolute value, besides GPAW where the basis set might be too small. Additionally, the QE results are comparable with the all-electron calculations (see Table 2). Due to reasons of reproducibility we use NRLMOL and FPLO to have independent all-electron codes for the calculation of J considering the variation of metal centers. For comparison of two different implementations of DFT (all-electron and plane wave) and because of the demonstration of the preliminary results, we additionally perform the screening with QE (see Table 3).
For a better visualization of the exchange coupling constants we introduce a coupling map (see Fig. 2). On the x-and y-axes the corresponding metal centers A/B are drawn, where the color  V and Fe in MOFs usually form chains and no paddle wheels, especially not the kind which is found in DUT-8(Ni). Ni, Co and Cu on the other hand are able to form such paddle wheel structures. This might lead to a stable crystalline structure by doping the DUT-8(Ni) with e.g. Fe. Another possibility is to use the proposed model systems as SBUs for a new kind of MOF, in which the magnetic characteristics of the model system are retained. For a deeper insight in the electronic structure of the HS SBUs we visualized the highest occupied molecular orbital (HOMO) and the lowest unoccupied molecular orbital (LUMO) for all Fe containing M1 HS systems (see Fig. 3). It is clearly visible that the HOMO/LUMO levels for all models are dominated by the d-orbitals of the metal centers. These results are confirmed by density of states (DOS) calculations (see Fig. 4). The DOS show clearly that the levels around the Fermi energy come solely from d-orbitals of the metal center. This result can be obtained from pseudopotential calculations (QE) as well (see ESI, † A4). Conversely, the model systems M1-[Fe,Cu] shows levels directly at the Fermi energy, where the other HS Fe models do not. In order to investigate this finding further, we carried out some further model calculations using NRLMOL. The model M1 derived from experimental data shows only a two-fold symmetry in the plane orthogonal to the metal center connecting axis. However, the deviation from a four-fold symmetry is rather small. We carried out calculations on a model system with higher   four-fold symmetry and found a doubly degenerate state at the Fermi level which was not completely filled. This indicates that the four-fold symmetry is indeed unstable and one would expect lower symmetry to cancel this degeneracy. The states at the Fermi level might be a hint to a metallic SBU and with that to the possibility to construct a metallic MOF (see M1-[Fe,Cu] states at E F in Fig. 4). However, further investigation in this respect are needed to prove whether such a metallic behaviour occurs in the crystalline structure or not. To investigate the magnetic stability of the HS states we carried out fixed magnetization calculations (see Fig. 5) using FPLO and QE. For each magnetization between the ones for the HS and LS state, the initial spin-orientations per atom were chosen to be either parallel or anti-parallel to each other, thus referring to the spin-polarization in the HS and the LS state, respectively. This was done to take into account different spin-orientations which lead to the same total magnetization. In When one of the Fe HS SBUs is used for a new MOF, the SBUs would be separated with mechanically rigid and stable bridges from each other and their geometry would be fixed (M1 -organic bridge -M1). With that it might be possible to retain the magnetic as well as the potential metallic behaviour (see M1-[Fe,Cu]) of the SBU inside this MOF.

Conclusions and outlook
A parameter free screening for the determination of HS solutions for the DUT-8(Ni) model system M1 has been performed using DFT. A comparison between different exchange-correlation functionals showed that PBE is sufficient for the calculation of the exchange coupling constant J. To summarize, we have shown that it is possible to find HS SBUs for a given MOF and a stable M1 Fe HS family has been derived for DUT-8(Ni). Furthermore these HS solutions are stable with respect to all magnetizations. We showed that all used DFT codes agree qualitatively as well as quantitatively. Thus it is possible to use the plane wave method for further investigations. Based on the results of this work it should be considered to insert the obtained HS SBUs into the crystalline structure of DUT-8(Ni), either as a complete replacement of the original SBUs to gain a fully ferromagnetic MOF or as a ''magnetic doping'' by replacing individual SBUs with some HS SBUs to introduce local HS sites. Additional calculations should consider the mechanical and thermal stability of the resulting structures (molecular dynamics and vibrational analysis) and the effect of structural relaxation on the exchange coupling constant. For first insights, we optimized the M1-[Fe,Fe] system and recalculated the coupling constant. We obtain J = 165.1 cm À1 for NRLMOL and J = 167.7 cm À1 for QE, respectively. Thus the magnitude of J only changes slightly due to geometry optimization and the high-spin character of the coupling is kept. Additionally we relaxed the M1-[Fe,Cu] model and the resulting HOMO-LUMO gap does not change significantly. With that the metallic behaviour of this model is stable against relaxation. Further investigations might include the search for a new kind of MOF, where the HS character of the SBUs is retained. Fig. 5 Relative energies of the Fe HS SBUs with fixed magnetizations calculated with FPLO and QE. E min is the energy of the most favourable state for each system. HS and LS mark the high-spin and low-spin states.