# Microscopic modeling of a spin crossover transition

###### Abstract

In spin crossover materials, an abrupt phase transition between a low spin state and a high spin state can be driven by temperature, pressure or by light irradiation. Of a special relevance are Fe(II) based coordination polymers where, in contrast to molecular systems, the phase transition between a spin and a spin state shows a pronounced hysteresis which is desirable for technical applications. A satisfactory microscopic explanation of this large cooperative phenomenon has been sought for a long time. The lack of X-ray data has been one of the reasons for the absence of microscopic studies. In this work, we present an efficient route to prepare reliable model structures and within an ab initio density functional theory analysis and effective model considerations we show that in polymeric spin crossover compounds magnetic exchange between high spin Fe(II) centers is as important as elastic couplings for understanding the phase transition. We discuss the relevance of these interactions for the cooperative behavior in these materials.

###### pacs:

75.30.Wx, 71.15.Mb, 71.20.Rv, 75.30.Et^{†}

^{†}: New J. Phys.

## 1 Introduction

An intensively debated class of materials with potential applications as optical switches, sensors or memory devices [1, 2, 3, 4] are spin-crossover polymer (SCP) systems involving transition metal ions linked with organic ligands [5]. These systems show a sharp transition triggered by variation of temperature, pressure or by light irradiation between a low-spin (LS) ground state and a high-spin (HS) excited state with a thermal hysteresis loop [6]. Especially important in these materials is the large cooperativity shown at the HS-LS transition in contrast to molecular spin crossover systems, for example. The origin of this transition and its cooperativity has been mainly discussed in the frame of elastic models [6, 7, 8, 9, 10], and only recently a possible role of magnetic exchange was suggested [11, 12]. Still, a conclusive ab initio microscopic study where all important interactions are considered is missing and the origin of the large cooperativity has not been completely settled. It is our purpose to investigate this issue in what follows. Ab initio theoretical studies for SCP systems are faced with major difficulties due to the nonexistence of well-determined crystal structures. To our knowledge, electronic structure calculations have only been performed at the level of semiempirical extended Hückel approximation for an idealized triazole-bridged Fe(II) chain [13]. In the present work, we overcome the unavailability of structural data by predicting a crystal structure for a Fe(II) spin-crossover polymeric crystal using known experimental constraints and a combination of classical force field and quantum mechanical molecular dynamical methods. We analyze with Density Functional Theory (DFT) calculations the LS-HS phase transition and show that there exists an interplay between magnetic exchange and elastic properties that is responsible for the large cooperativity in these systems. We also corroborate the quality of our designed structure by comparing with magnetic experiments done on a real sample. Our methodology and results provide a new perspective on the parameters underlying the traditional theoretical approaches.

There have been a number of attempts to theoretically account for the features of the HS-LS transition in spin crossover materials. Most of the theoretical work is based on elastic considerations. Two types of elastic models that focus on the faithful reproduction of macroscopic quantities like the HS fraction, are the following: In the first approach, the cooperativity in the HS-LS transition is defined in terms of local distortions which interact with one another elastically causing a long range effective interaction between spin states. This leads to an Ising-type Hamiltonian [14, 15, 16] which describes the elastic interaction between spin states (LS and HS) in terms of fictitious spin operators ( = -1 (1) for LS (HS)) coupled via a nearest neighbor interaction . The coupling constants are parameters of the theory and have not yet been determined from a microscopic model. In an alternative approach, the free energy of spin crossover systems is calculated based on an anisotropic sphere model that describes volume and shape changes of the lattice at the transition [7, 8, 9]. One of the few attempts to include magnetic interactions is the recent proposal done in Ref. [12] where the author considers a one-dimensional model for HS-LS systems which contains elastic and magnetic Ising exchange interactions. The ground-state phase diagram is then obtained by the transfer matrix technique for different relative elastic to magnetic coupling strengths.

In the present work we concentrate on the electronic and magnetic degrees of freedom in a spin crossover polymer, and investigate their influence on the microscopic origin of the cooperativity in the HS-LS transition. While it has been assumed in the past that the Fe(II) center nearest neighbor interaction is mostly of phononic origin, our study indicates that a significant part of this interaction arises from magnetic exchange.

## 2 Crystal Structure

For this investigation it is indispensable to obtain a reliable crystal structure suitable for DFT analysis. We aim at describing the complex Fe[(hyetrz)](4-chlorophenylsulfonate)HO [17] ( stands for ) (see compound 1 in Fig. 1), which was synthesized from 2-hydroxyethyltriazole and iron(II)-p-chlorobenzenesulfonate as described in Ref. [18].

This compound precipitates as a fine, polymeric powder and single crystals cannot be grown because the polymer is insoluble in water and organic solvents [19]. Melting or sublimation attempts result in decomposition. Hence, X-ray structure analysis is not possible and even the X-ray powder diagram consists of a few broad peaks only which prevent the structure to be determined from X-ray powder data. Other polymeric Fe-triazole compounds have similar properties and no single crystal structure for these systems is known. In the literature only single crystal structures for trimeric Fe compounds (e.g. Ref. [20]), and polymeric Cu-triazoles [21, 22] were found. In these compounds, the metal ion is coordinated by six nitrogen atoms and neighboring metal ions are connected by three pyrazole bridges. For the polymeric Fe triazoles a similar structure is assumed, see Fig. 1. This structure is also supported by spectroscopic methods [22] including solid state NMR [23].

In view of the above and the requirement of having reliable crystal structures for microscopic studies, we design on the computer a model system of polymeric Fe triazole, as a basis for calculating the electronic and magnetic properties. We consider all experimental information available and construct a crystal structure as close as possible to the actual structure. We employ a method that has been developed and tested on the coordination polymer Cu(II)-2,5-bis(pyrazole-1-yl)-1,4-dihydroxybenzene which has a simpler structure and less atoms per unit cell than polymeric Fe-triazole [24, 25].

Since we aim at understanding the HS-LS transition with accurate all-electron DFT calculations, which are computer intensive, we keep the essential features of the material and simplify those elements that are secondary to the transition, like the nature of the substituent R and the counter ion X (spin transitions are observed for a wide range of different substituents R and counter ions X). We consider compound 2 in Fig. 1 with R = CH and X = F. This model structure has the short range environment of the Fe(II) centers exactly as in structure 1 in Fig. 1 while the longer range environment ( Å) of the Fe centers is significantly simplified. The molecular geometry of compound 2 was constructed according to data from single crystal analysis of dimeric and trimeric Fe-triazole compounds [26]. A hypothetical, but sensible crystal structure of compound 2 was built up with the minimum number of atoms per unit cell (72 atoms). The iron-triazole chain itself has symmetry. A crystal structure with hexagonal symmetry is in principle possible, e.g. in space group , but we chose a structure with symmetry with two formula units per unit cell. A similar arrangement of chains is also found in a corresponding Cu polymer, [Cu(hyetrz)](CFSO)HO [21]. (This structure is triclinic, with space group , , but the deviations from monoclinic symmetry are probably caused only by the symmetry-breaking hydroxyethyl ligands and CFSO counterions. Otherwise the symmetry would be .) In our structure, the Fe ions are located on crystallographic inversion centers, whereas all triazole units contain a crystallographic mirror plane between the neighboring nitrogen atoms. All Fe ions are crystallographically equivalent and we enforce, for simplicity, a perfectly octahedral environment of the Fe ions.

In the design and analysis of the model structures we employ four distinct classical and ab initio methods. For preoptimization of the model structures we consider a classical modified Dreiding force field [27] with atomic charges calculated by the Gasteiger method [28] as implemented in the Cerius package [29]. We then perform quantum mechanical first principles calculations within DFT with three different basis sets, each for a distinct purpose: The plane wave ab initio molecular dynamics (AIMD) method [30] is used for determination of precise equilibrium structures. The linearized augmented plane wave (LAPW) method [31] is used to determine accurate electronic and magnetic properties and the linearized / Nth order muffin tin orbital (LMTO/NMTO) methods [32] are used to calculate effective Fe Wannier functions and to understand the low energy excitations of the system.

For the AIMD calculation [30] we considered a plane wave cutoff of 30 Ryd for the plane wave part of the wave function and we used the following sets of (s,p,d) projector functions per angular momentum: Fe(2,2,2), F(2,2,1), N(2,2,1), C(2,2,1) and H(2,0,0). We employed a mesh and the symmetry was preserved during the relaxation with the help of 131 constraints.

For the LAPW calculations [31], we employed a mesh and a value of that is sufficiently large due to the small radii of the hydrogen atoms.

The NMTO-downfolding [32] calculations, which rely on the self-consistent potentials derived out of LMTO [33] calculations, were carried out with 40 different empty spheres in addition to atomic spheres to space fill. The convergence of LMTO calculations in each case was cross-checked with full potential LAPW calculations.

All DFT calculations were performed within the Generalized Gradient Approximation [34].

We first proceed with the relaxation of the model structure with the classical force field within the space group (non-standard setting of ), keeping fixed the Fe-N distances of the FeN octahedra to the values Å, Å, Å, Å, Å, Å, Å, Å. This constraint determines the value of crystal field splitting from the onset and constitutes our control parameter as shown below. As the N-N bond of the triazole (see Fig. 1) has experimentally a well defined bond length of Å, the choice of a Fe-N distance leads automatically to lattice parameters linear in . The relative change of the lattice parameter (the chain direction) is larger than that for and . This means that the volume of the unit cell also increases linearly from Å at Å to Å at Å. It should be noted that the force field we use is not suitable for relaxing the F anion into an equilibrium position due to the lack of appropriate parameters for F anions. Therefore, the F counter-ion is put into a likely position and the optimization is left to the next step where we perform a precise relaxation of the structure with the help of ab initio molecular dynamics. This AIMD step is essential as we find that the force field relaxed structures still show no or very bad convergence in LAPW, indicating inappropriate positions of at least some atoms. We relax all unconstrained atoms within AIMD until the forces that initially are of the order of mRyd/Bohr until they are far below mRyd/Bohr. More significantly, we converge bond lengths to within Å and bond angles to within 0.1. We cross check the final AIMD relaxed structures by calculating the LAPW forces and making sure that they are of the order of mRyd/Bohr or less (with the only exception of the difficult to place counter ions F for which we allow mRyd/Bohr as they do not play a crucial role for the Fe(II) center interactions). The stipulation that forces are very low in two fundamentally different DFT codes is a high standard and gives us confidence that our conclusions about the electronic and magnetic structure are drawn on a solid basis. The precision to which we converge bonds and angles means that we can then proceed to predict interactions between Fe(II) centers with the necessary precision of a few Kelvin. In Fig. 2 a representative of the resulting structures is shown. The top panel pictures the chain of FeN octahedra with alternating orientations, and in the bottom panel we demonstrate the arrangement of the Fe(II) chains in the crystal.

## 3 Energy scales

At the atomic level, two sets of energy scales are responsible for the LS and HS state of the Fe(II) centers, namely the crystal-field splitting and the Hund’s exchange and Coulomb interactions.

The crystal-field splitting, as mentioned above, is being fixed beforehand in the construction of the structures with given distances. In Fig. 3 we show the dependence of the crystal field splitting on the size of the FeN octahedra calculated with LAPW (magenta symbols) and NMTO (green symbols). Note the good agreement between the two calculations. The crystal field splitting values are obtained by determining the first moment of both and densities of states from non-spin polarized LAPW calculations and through construction of an Fe d only Hamiltonian in case of NMTO calculations. The atomic state diagrams for Fe 3d schematically demonstrate the relationship between crystal field splitting and the spin state ( or ). The Hund’s exchange is taken into account to some extent within the spin polarized-GGA approach.

In the extended system, two type of interactions contribute to the phase transition; the phononic excitations and the exchange interaction , due to nearest neighbor superexchange between Fe(II) centers which is typically antiferromagnetic.

The competition between all these energy scales determines the nature of the phase transition and its cooperativity, as we will discuss in section 6. The role played by the phonons in driving the LS-HS transition in the spin-crossover systems has been discussed at length in terms of elastic models [7, 8, 10, 14, 15]. In our model calculations, the phononic degrees of freedom are frozen and we investigate the role of the electronic and magnetic degrees of freedom.

## 4 Electronic structure

We performed spin resolved electronic structure calculations for these systems within DFT by considering the GGA [34] as exchange-correlation functional and the LAPW basis. In Fig. 4 we present for four of the designed structures the projection of the total density of states (DOS) on the Fe orbitals, which are responsible for the magnetism of the material. The first and second panels correspond to the structures with Å and Å respectively and show a perfect symmetry between spin up (red) and spin down (blue) density of states and therefore define a LS state (). The occupied states can be identified as the six states, while the empty states are the four states. The splitting between and states diminishes with increasing from eV to eV, respectively (see also Fig. 3). The DOS behavior completely changes as the Fe-N distance increases to Å (see the third panel of Fig. 4). Now spin up (red) and states are completely filled, and spin down (blue) states show only a partial filling with one electron. The imbalance between up and down electron numbers is which corresponds to the HS state (). This situation remains if the Fe-N distance increases further to Å, only the splitting between and diminishes.

Thus, by carefully preparing a series of model structures that correspond to the LS and HS sides of the spin crossover transition we manage to microscopically describe the LS-HS transition which occurs between the structures with Å () and Å ().

In order to quantify energetically the HS-LS spin transition, we show in Fig. 5 the total LAPW electronic energies obtained within the spin-polarized GGA (sp-GGA) approach. We note that there is a discontinuous jump between the LS (S=0) energies and the HS (S=2) energies. The relative electronic energy differences between the HS and LS systems is about which agrees with the relative energy estimates for spin crossover molecular systems [35]. Since the LS-HS phase transition occurs between structures Å and Å, we designed one more structure with Å in order to probe the sharpness of the transition. While there are indications that this structure might represent an intermediate magnetic state with per Fe(II) center, we do not include it in Fig. 5 as it is very hard to converge. This result indicates that the spin crossover transition in polymer systems occurs in a very narrow range of crystal field splittings, i.e. is very sharp, as observed experimentally [5].

## 5 Exchange interaction and magnetic properties

For the magnetic behavior of this system, we derive first from the DFT electronic calculations a Hamiltonian which describes the effective interaction between Fe(II) centers. The NMTO-based downfolding method, designed to pick up selectively the low-energy bands from the whole group of LDA/GGA bands of a compound, has been used to construct the Fe only Hamiltonian of the Fe-triazole compounds. This method has proven to provide reliable information on the low energy properties of inorganic [36, 37, 38] and organic [25] transition metal compounds. The tight-binding basis in which these Hamiltonians are constructed form the set of "effective" functions, which span the Hilbert space of the Wannier functions corresponding to the low-energy bands. Fig. 6 shows the plot of one member of such a set, namely the downfolded Wannier function corresponding to Fe . In the figure, two such Fe Wannier functions have been placed at two neighboring Fe sites. While the central part of such an effective function has the Fe symmetry, the tails of the function are shaped according to integrated out degrees of freedom in the system, like C , N , F and H . As is evident from the plot, substantial weight of these tails resides on neighboring triazole rings. The presence of these tails hints to an enhanced electronic hybridization between the adjacent Fe ions, contributing to the cooperative nature of the HS-LS transition. Out of these calculations we can estimate the various hopping matrix elements, , between the orbitals of adjacent Fe(II) centers. The values of these hopping parameters range between 1 meV to 80 meV quantifying the strength of the various interaction paths between neighboring Fe(II) orbitals.

In order to get estimates of the magnetic superexchange coupling
constants between neighboring Fe(II) centers
for the HS Fe-triazole structures, we considered two
approaches:

i) Total energies. We calculated within sp-GGA total energies of
ferromagnetic and antiferromagnetic Fe spin
configurations. Considering a spin-Hamiltonian
between nearest neighbors Fe spins , the ferromagnetic and
antiferromagnetic energies for two Fe ions in the unit cell of
the Fe-triazole are given by and
. Within sp-GGA meV for the Å structure and therefore .
ii) Perturbation theory [39] where can be
obtained in terms of the hopping parameters between Fe(II) centers
and the onsite Coulomb repulsion as . For the HS Å structure, the significant
obtained within the downfolding method is 48 meV, and for = 4-5 eV
this gives
which is very similar to the value obtained
with the difference of total energies.

The results of our model calculations can be now compared with the magnetic properties obtained on the real samples of polymeric Fe[(hyetrz)](4-chlorophenylsulfonate)HO. Variable-temperature magnetic susceptibility measurements in the temperature range K and magnetic fields T were carried out on powder samples of Fe(II) triazole using a Quantum Design SQUID magnetometer MPMS-XL. In Fig. 7 we show the susceptibility measurements where has been plotted versus . Our sample shows hysteresis at K with a width of 20 K. Since this system consists of spin Fe(II) chains with weak interchain interactions, we have analyzed the magnetic susceptibility in the frame of a spin Heisenberg chain model.

Note that the Fe(II) triazole SCP systems are Haldane chains. A Haldane gap is expected to exist between the ground state and first excited states. The reason for not observing this gap, is that the transition temperature at which the HS-LS transition happens is higher in energy than the Haldane gap energy and therefore, the system goes into the non-magnetic phase before the gap in the chain can be observed.

The susceptibility of an -site chain is given by:

(1) |

where is the Heisenberg Hamiltonian , is the Bohr magneton, is the gyromagnetic factor, is the Boltzmann constant, is temperature and is the -component of the spin on site . In the thermodynamic limit, the bulk susceptibility at high temperatures can be obtained as a series expansion in :

(2) |

For Eq. 2 compares very well to QMC data for spin chains [40].

The values obtained from our ab initio calculations ( ) are slightly larger than the extracted from the susceptibility data ( K) on the real sample, but remain in the same order of magnitude. Considering that i) we performed the calculations in a model structure and ii) the experimental measurements are affected by the quality of the samples, we can conclude that the comparison is quite good and the designed structures are reliable.

## 6 Discussion

One important topic of this work is the analysis of the various energy scales that contribute to the cooperativity of the HS-LS transition in SCP systems. In our calculations we froze the elastic degrees of freedom and concentrated on the electronic and magnetic properties for which we have quantified the corresponding parameters. A comparison with the elastic coupling constants estimated from Ising-like models [15] (in such an approach, a material with a transition temperature of K would be described by elastic interactions of K) shows that in one-dimensional Fe(II) triazole, the elastic coupling constants are of equal importance as the magnetic exchange for explaining the spin crossover transition, and the cooperativity should be understood as an interplay between elastic properties and magnetic exchange. When we cool the system from the HS state toward the HS-LS transition temperature, the elastic coupling tends to drive the system to the LS state, while the magnetic exchange tends to keep up the magnetic state for a larger temperature range (see Fig. 7). In comparison, in the heating process, the magnetic exchange is initially absent (LS) and therefore the elastic interaction (vibrational phonons) initially drives the transition which has its on-set at a higher temperature than in the cooling process. The width of the transition between cooling and heating (hysteresis) is therefore enhanced by the magnetic interaction.

A fundamental difference between the polymeric systems we are dealing with in this work, and molecular bi- (tri-, tetra-, …) nuclear Fe systems is the connectivity between the Fe(II) centers. While the molecular systems [41] form isolated clusters of Fe(II) centers and therefore there is no strong connectivity between clusters, the polymers have important nearest neighbor interactions in the thermodynamic limit. This implies that for the polymers, the magnetic superexchange is not restricted to the cluster as in the molecular systems, but rather becomes important for the nature of the HS-LS phase transition. In Ref. [42] various estimates of the magnetic J for molecular systems have been given. The values range between 4-6 K. The values we estimated for the present polymers are larger, between 11-24 K and in the energy range of the elastic constants, which indicates that the cooperativity in these SCP systems is most likely significantly enhanced by the exchange interactions.

In conclusion, this work presents an efficient route to prepare reliable model structures for microscopic investigations and provides a new interpretation about the origin of the parameters underlying traditional theoretical approaches for the polymeric spin-crossover materials. The next step of this study will be the inclusion of lattice dynamics (phonons) in the ab initio calculations which is planned in our future work.

## References

## References

- [1] Cobo S, Molnar G, Real J A, Bousseksou A 2006 Angew. Chem. 45 5786–5789
- [2] Freysz E, Montant S, Létard S, Létard J F 2005 Chem. Phys. Lett. 394 318–323 .
- [3] Létard J-F, Guionneau P, Goux-Capes L 2004 Top. Curr. Chem. 235 221–249
- [4] Garcia Y, Ksenofontov V, Gütlich P 2002 Hyperfine Interactions 139/140 543–551
- [5] Kahn O, Jay Martinez C 1998 Science 279 44–48
- [6] Gütlich P, Hauser A, Spiering H 1994 Angew. Chem. Int. Ed. 33 2024–2054
- [7] Willenbacher N, Spiering H 1988 J. Phys. C: Solid State Phys. 21 1423-1439
- [8] Spiering H, Willenbacher N 1989 J. Phys.: Condens. Matter 1 10089–10105
- [9] Spiering H, Boukheddaden K, Linares J, Varret F 2004 Phys. Rev. B 70 184106
- [10] Nishino M, Boukheddaden K, Konishi Y, Miyashita S 2007 Phys. Rev. Lett. 98 247203
- [11] Timm C, Schollwöck U 2005 Phys. Rev. B 71 224414
- [12] Timm C 2006 Phys. Rev. B 73 014423
- [13] Yoshizawa K, Miyajima H, Yamabe T 1997 J. Phys. Chem.B, 101, 4383–4389
- [14] Wajnflasz J, Pick R 1971 J. Phys.Colloq. France 32 C1-91
- [15] Boukheddaden K, Miyashita S, Nishino M 2007 Phys. Rev. B 75 094112
- [16] Boukheddaden K, Linares J, Tanasa R, Chong C 2007 J. Phys.: Condens. Matter 19 106201
- [17] Garcia Y, Niel V, Muñoz M C, Real J A 2004 Top. Curr. Chem. 233 229–257
- [18] Garcia Y, van Koningsbruggen P J, Codjovi E, Lapouyade R, Kahn O, Rabardel L 1997 J. Mater. Chem. 7 857-858
- [19] Gütlich P, van Koningsbruggen P J, Renz F 2004 Struct. Bonding 107 27–75
- [20] Garcia Y, Guionneau P, Bravic G, Chasseau D, Howard J A K, Kahn O, Ksenofontov V, Reiman S, Gütlich P 2000 Eur. J. Inorg. Chem. 1531–1538
- [21] Garcia Y, van Koningsbruggen P J, Bravic G, Chasseau D, Kahn O 2003 Eur. J. Inorg. Chem. 356–362
- [22] Garcia Y, van Koningsbruggen P J, Bravic G, Guionneau P, Chasseau D, Cascarano G L, Moscovici J, Lambert K, Michalowicz A, Kahn O 1997 Inorg. Chem. 36 6357–6365
- [23] Glaubitz C, Brüning J, Schmidt M U, unpublished results.
- [24] Jeschke H O, Salguero L A, Valentí R, Buchsbaum C, Schmidt M U, Wagner M 2007 C. R. Chim. 10, 82–88
- [25] Salguero L A, Jeschke H O, Rahaman B, Saha-Dasgupta T, Buchsbaum C, Schmidt M U, Valentí V 2007 New J. Phys. 9 26
- [26] Cambridge Structural Database, Cambridge Crystallographic Data Centre, Cambridge, England (2005).
- [27] Mayo L, Olafson B D, Goddard III W A 1990 J. Phys. Chem. 94 8897–8909
- [28] Gasteiger J, Marsili M 1980 Tetrahedron 36 3219–3228
- [29] Cerius, Version 4.9, Accelrys Inc., Cambridge, England (2003).
- [30] Blöchl P E 1994 Phys. Rev. B 50 17953–17979
- [31] Blaha P, Schwarz K, Madsen G K H, Kvasnicka D, and Luitz J, WIEN2k, An Augmented Plane Wave + Local Orbitals Program for Calculating Crystal Properties (Karlheinz Schwarz, Techn. Universität Wien, Austria), 2001. ISBN 3-9501031-1-2.
- [32] Andersen O K, Saha-Dasgupta T 2000 Phys. Rev. B 62 R16219–R16922
- [33] Andersen O K, and Jepsen O 1984 Phys. Rev. Lett. 53 2571–2574
- [34] Perdew J P, Burke S and Ernzerhof M 1996 Phys. Rev. Lett. 77 3865–3868
- [35] Paulsen H, Trautwein A X 2004 Top. Curr. Chem. 235 197–219
- [36] Valentí V, Saha-Dasgupta T, Alvarez J V, Požgajčić K and Gros C 2001 Phys. Rev. Lett. 86 5381–5384
- [37] Saha-Dasgupta T, Valentí R, Capraro F, Gros C 2005 Phys. Rev. Lett. 95 107201
- [38] Rahaman B, Jeschke H O, Valentí V, Saha-Dasgupta T 2007 Phys. Rev. B 75 024404
- [39] Anderson P W 1978Rev. Mod. Phys. 50 191–201
- [40] Yamamoto S 1996 Phys. Rev. B 53 3364–3373
- [41] Postnikov A V, Kortus J, Pederson M R 2006 phys. stat. sol. (b) 243 2533–2572
- [42] Desroches C, Pilet G, Szilágyi P A, Molnár G, Borshch S A, Bousseksou A, Parola S, Luneau D 2006 Eur. J. Inorg. Chem 2006 357–365