Combination of density functional theory and calorimetry reveals the microscopic nature of spin state switching in 1D Fe(II) spin crossover complexes

Abstract

The Gibbs free energy of spin transitions in the heptanuclear models of 1D Fe(II) spin crossover 1,2,4-triazole complexes has been estimated using DFT methods. The complexes modelled were [Fe(Htrz)₂trz]BF₄ (Htrz = 4H-1,2,4-triazole, trz = 1,2,4-triazolato), 1 the dehydrated and hydrated [Fe(NH₂trz)₃]Cl₂ (NH₂trz = 4-amino-1,2,4-triazole), 2 and 2a, and [Fe(NH₂trz)₃](NO₃)₂, 3. For each complex, the electronic energy and the vibrational energies were calculated for a heptanuclear model containing five inner Fe(II) centres in the high-spin (HS) and the low-spin (LS) states. All other possible 18 spin isomers with one to four HS centres were also modelled. Results obtained using different exchange–correlation functionals based on the B3LYP one show that each spin isomer with a particular permutation of HS and LS centres within the pentanuclear linear unit has distinctive electronic and vibrational energies. The electronic energy of each spin isomer was found to be equal to the sum of the adiabatic electronic energy of the spin transition E_ad given by the difference in energies between the LS and HS states and the strain energy H_strain. This quantity is non-zero for any spin isomer containing both LS and HS centres. Unlike E_ad, which has also been determined experimentally by calorimetric measurements, H_strain is independent of the applied functional. Calculations of the temperature dependence of the Gibbs free energy change ΔG of 19 possible spin transitions for heptanuclear model systems reveals that strain effects lead to an additional destabilisation of the spin isomers containing both LS and HS centres. The actual strain pattern depends on the chemical structure of the model molecule.

---

Introduction

Dynamic lattice properties lie at the heart of the spin crossover (SCO) phenomenon.¹ For a molecular crystal, the low energy acoustic phonons show movements of the molecules as a whole, while the optical phonons are characterized by more localized molecular vibrations that occur at higher energies. Phonons determine the key properties such as the heat capacity, vibrational entropy and energy, lattice rigidity, elastic constants and elastic interactions.² All static³ or dynamical⁴ phenomenological models aiming to explain the cooperativity of spin transition and its spatiotemporal behaviour consider these magnitudes in more or less explicit form. As Palii, Klokishner et al.⁵ state, two general trends in modelling of cooperative spin transitions can be identified:

(i) The macroscopic approach focuses on the examination of elastic interactions that change when the spins switch. These interactions result in the coupling of spin centres as they are transmitted to the lattice. One of the most effective approaches is the mechano-elastic approach, which utilises molecular dynamics (MD) and Monte Carlo (MC) methods.⁶

(ii) The microscopic approach concerns the molecular electronic structure and normal modes, while the coupling between the electronic states of the SCO centres is transmitted by phonons.

In the light of the evolution of first-principle modelling of SCO materials over the last two decades, particularly the development of density functional theory (DFT)⁷ including functional performance⁸ and modelling of structural-spin transition relationships⁹ a pertinent question emerges: how can quantum chemistry provide links between the molecular properties of the SCO centre and the relevant parameters of microscopic and macroscopic models? For example, what is the relationship between the structural, electronic and vibrational properties of the high-spin (HS) and low-spin (LS) isomers of an SCO molecule and the elastic tensor components of the macroscopic model, or the coupling parameters of the microscopic model like that used in ref. 5. This is particularly relevant in the context of molecular material science, where a key problem is determining how the chemical composition of molecules influences the mechanical properties of the material from which they are composed.

In principle, quantum state calculations for solids (periodic calculations) should result in optimised geometries, cell volumes, electronic states, energies and vibrations. This would provide a thorough thermodynamic characterization of phases with different proportions of HS centres (x_HS) and distinct spin centre distributions. According to a recent review the vast majority of calculations relating to the thermal dependence of the elastic tensor components have been conducted on simple organic molecules.¹⁰ However, the work of Erba et al.¹¹ demonstrates that a quasi-harmonic method could provide a computationally affordable approach. The parameters and their temperature dependence, as determined by X-ray methods¹² were utilised to estimate the direct interaction and self-energy contributions to the cooperativity parameter of Spiering's model.¹³ Thus, in principle, periodic calculations that yield magnitudes relevant to phenomenological models for a SCO complex are feasible. However, it must be noted that such calculations are hindered by two issues. The primary issue relates to the calculation of spin transition energy (see ref. 14 for the most recent review). The electronic energy of the spin transition, as calculated by DFT methods, depends on the exchange–correlation functional employed,¹⁵ even in the case of isolated molecules. This issue can be addressed in two ways. Firstly, a series of related systems can be compared, differing in terms of substituents or solvated molecules. Secondly, the experimental values of spin transition energies can be considered (ref. 14 and vide infra). The first approach has also been used for periodic calculations involving complexes with different anions¹⁶ and solvated molecules.¹⁷ Importantly, when considering different spatial combinations of a given number of HS and LS molecules, the dependency ceases to exist, while systems of the same multiplicity are modelled. In their study, Vela and Paulsen considered the approach of considering the different spatial distribution of HS and LS centres to ascertain the energy contribution of cooperativity for [Fe(2-pic)₃]Cl₂ (2-pic = 2-picolylamine)¹⁸ and the phenomenological interaction parameter Γ for [Fe(phen)₂(NCS)₂].¹⁹ Additionally, modelling of the SCO in the 3D Hoffman-clathrate [Fe^II(TPB){Au^I(CN)₂}]I·4H₂O·4DMF (TPB = 1,2,4,5-tetra(pyridin-4-yl)benzene)²⁰ was performed by calculating electronic energies of HS or LS “spin defects” by introducing one LS or one HS centre in a supercell containing eight SCO centres. For periodic energy calculations, the Hubbard-like parameter U is typically used to describe intra-atomic Coulomb repulsions. This is done within the framework of the DFT + U approach, which introduces an additional energy term, U, for the localised electrons, while describing the delocalised electrons using DFT functionals alone.^21–23

Another issue is the calculation of the vibrational modes, which determine the vibrational entropy that drives the SCO process. In this case, there is an inevitable conflict between achieving precise calculations and reducing computational time. It has been demonstrated that the latter increases in proportion to the size of the modelled molecules and the level of DFT tools utilised, as well as the algorithm employed for calculating Hessian matrix elements. The initial papers on periodic calculations of SCO complexes illustrate this. While modelling of the CsFe^II[Cr^III(CN)₆] Prussian blue analogue²⁴ yielded only three small imaginary frequencies (note that the electronic energies for this system were previously calculated with the GGA + U approach²⁵ using the finite-differences method to assess the Hessian matrix elements), a paper reporting a similar approach to the [Fe(phen)₂(NCS)₂] complex mentions the general problem of imaginary frequencies related to a number of soft transitions below 50 cm⁻¹.²⁶ Therefore, achieving accuracy is likely to be more challenging for molecular crystals due to their weaker intermolecular interactions compared to ionic crystals. In a subsequent report Zhang²⁷ calculated the electronic energy of spin transitions of molecular crystals of the Schiff-base ligand complex with an FeN₄O₂ coordination core. This study employed the LDA/GGA functionals within the DFT + U approach using Quantum-ESPRESSO (QE) package.²⁸ SIESTA²⁹ was used to calculate the harmonic vibrations using the finite-differences approximation (see ref. 30, where the accuracy and potential perspectives of this method are discussed). Subsequently, Vela et al.²² employed the same approach with QE for several additional SCO molecular crystals of Fe(II). They indicated that computing both analytical frequencies and frequencies obtained with the finite-difference approximation at the crystalline phase would be prohibitive. To assess the vibrational entropy, the authors used values calculated for single molecules with functionals that go beyond the LDA/GGA. This approach was further developed in subsequent ref. 17–19. Collet et al.³¹ reported the use of terahertz, infrared and Raman spectroscopy on HS [Fe(phen)₂(NCS)₂] crystals, together with periodic calculations (GGA) using the VASP package.³² Calculating vibrations using the frozen-phonon method yielded reasonable energies in the low-frequency region of 20–270 cm⁻¹. Recently, Poloni, Rodríguez-Velamazán et al.³³ reported PBE + D2 calculations of the phonon density of states (pDOS) with QE for the Fe(pz)[Pt(CN)₄] (pz = pyrazine) 3D SCO complex and its SO₂ adduct, obtaining a fit to inelastic neutron scattering (INS) data. Harmonic interatomic force constants were computed using density functional perturbation theory.³⁴ Recently Wu and coworkers³⁵ demonstrated that the PBE + U3 + D3 approach for periodic calculations allows the reproducing of the spin-transition temperatures in a series of 2D complexes. In summary, despite recent advances in periodic DFT calculations for SCO materials, these calculations still require substantial computing time and cannot be used with functionals beyond the GGA approximation. This has implications for the accuracy of low-frequency phonon mode energies, particularly in the context of molecular crystals.

Consequently, modelling the thermodynamics of multiple spatial distributions of spin centres for various spin isomers is an alternative approach for the time being. Another possibility is to use MD to model crystals involving millions of atoms.³⁶ The quality of the calculations can be improved by using DFT calculated force constants. This has resulted in a favourable alignment between the experimental and MD calculated pDOS for models of the 1D Fe(II) polynuclear SCO chains based on 1,2,4-triazole-type ligands.³⁷ In a recent report the Bousseksou group demonstrated that MD calculations (with no DFT calibration of the force constants) for the crystals of the 3D SCO compound Fe(pz)[Ni(CN)₄] successfully reproduced the thermodynamic and mechanical parameters, as well as the pDOS of both, the LS and HS phases.³⁸ Importantly, this approach enabled surface effects to be estimated, which are essential for explaining the hysteretic behaviour of the spin transition.¹³ MD has been combined with the MC approach^39,40 or with MC and ligand field molecular mechanics⁴¹ to model the spin transition energy and mechanisms in the above mentioned Fe(pz)[Pt(CN)₄] 3D SCO complex. It is reasonable to assume that this approach will provide superior modelling of the soft phonons if the lattice under investigation comprises covalently bound units rather than being a molecular crystal.

The anharmonicity of vibrations restricts the accuracy of the harmonic approximation, especially regarding the thermal lattice expansion of the system and its thermoelasticity.⁴² Another issue is that the harmonic approximation neglects intrinsic phonon anharmonicity and phonon–phonon coupling, which has consequences for the computed thermodynamic properties. These problems can be overcome using the quasi-harmonic approximation (QHA).⁴² Alternatively, optimised scale factors can be applied when calculating vibrational, harmonic and fundamental frequencies, as well as zero-point energies (ZPE).⁴³ Additionally, a composite scheme has been proposed that considers anharmonicity when calculating molecular entropy by combining DFT, semi-empirical and force-field approximations.⁴⁴

In this context, it is noteworthy that calculations of anharmonic effects for spin crossover systems are feasible with commonly used software packages, despite being computationally much more demanding than calculations of harmonic frequencies.⁴⁵ In this work we neglect anharmonic effects according to the following quotation of ref. 45: “as a conclusion, we put forward that for high precision results, one should be aware of the anharmonic effects, but as long as computational chemistry is still struggling with other larger factors like the influence of the environment and the accurate determination of the electronic energy difference between HS and LS, the anharmonicity of the vibrational modes is a minor concern”. Conversely, accounting for anharmonic effects is essential for phenomenological models of elasticity in SCO systems.^13,46

This presentation outlines the results of DFT modelling of SCO assemblies with various spin-state permutations. This approach sheds light on the differences in thermodynamic and vibrational parameters between phases with different distributions of HS and LS centres. We used 1,2,4-triazole derivatives as model compounds to simulate one-dimensional (1D) systems involving pentanuclear, hexanuclear and nonanuclear chains.⁴⁷ This enabled us to achieve a satisfactory correlation between the ⁵⁷Fe-pDOS derived theoretically and that obtained experimentally by nuclear inelastic scattering (NIS).⁴⁸ Furthermore, modelling 3D Fe(pz)[Pt(CN)₄] with a cubic cell containing 15 SCO centres yielded good agreement between the calculated ⁵⁷Fe-pDOS and the experimental pDOS obtained from nuclear inelastic scattering (NIS) measurements.⁴⁹ We also present calculations of the temperature dependence of the entropy difference ΔS and the Gibbs free energy difference ΔG of the spin transition for different spin isomers.

To quantify intramolecular cooperativity, we have introduced the parameter H_coop, defined as the difference of the LS → HS spin transition electronic energy of a SCO centre with the neighbours being in LS (L) and HS (H) state,^50,47d i.e. between LLLLL → LLHLL and HHLHH → HHHHH transitions energies. H_coop was shown to be essentially independent on the applied functional. In this report we introduce a novel parameter, H_strain, which quantifies the strain present in each possible spin isomer and show that H_coop contains this quantity for two particular spin isomers. A related parameter, H_block, was recently introduced, which can be used to determine the cooperativity of the spin transition in mononuclear complexes based on the electronic energies of the LS and HS isomers within matrices that have been optimised using periodic calculations for a given spin of the switching centres.⁵¹ To quote the authors “H_block quantifies the structural quenching of a SCO molecule to remain in a given spin-state”.

The obvious cost of modelling the linear 1D SCO chain materials with a single chain of finite size is that the hysteretic character of the spin transition cannot be predicted since it is related to long-range interactions.¹³ In other words, the sharp transitions in such systems take place in 2D and 3D.⁵² Indeed, the macroscopic mean-field approach reveals that the first-order (sharp and hysteretic) transitions occur when both, the short-range (intrachain) and long-range (interchain) interactions are ferroelastic (i.e. the pure LS and HS phases are thermodynamically preferred). The two-step transitions occur with antiferroelastic (i.e. the alternating LHLHLH-like phases are thermodynamically preferred) short-range interactions and ferroelastic long-range interactions. The importance of short-range ferroelastic interactions for the first-order spin transition has been recognized.⁵³ Although the modelling of 1D SCO chains cannot account for the long-range (interchain) interactions, one may expect a fairly accurate modelling of the short-range intrachain interactions. Hence, such simulations could predict which of the following scenarios will occur for a given 1D SCO chain:

(a) The pure LS and HS phases correspond to the lowest Gibbs free energy as a function of temperature. This corresponds to the dominant ferroelastic short-range interactions and is a prerequisite for a sharp transition (enhanced by the ferroelastic long-range interactions).

(b) The alternate LHLHLH-like phase emerges as the new ground state at some temperature replacing the pure LS phase. This corresponds to the dominant antiferroelastic short-range interactions and is a prerequisite for a two-step transition (again enhanced by the ferroelastic long-range interactions).

(c) None of the above takes place, the phases with different molar fractions of HS appear to have the lowest Gibbs free energy with increasing temperature. Neither elastic nor antiferroelastic interactions dominate and either a soft or an incomplete spin transition occurs.

The current report describes the results of the DFT modelling of the oligonuclear models of the 1D SCO chains of Fe(II) complexes with 4H-R-1,2,4-triazole ligands, namely [Fe(Htrz)₂trz]BF₄ (Htrz = 4H-1,2,4-triazole, trz = 1,2,4-triazolato) (1), [Fe(NH₂trz)₃]Cl₂ (NH₂trz = 4-amino-1,2,4-triazole) (2), [Fe(NH₂trz)₃]Cl₂·2.5H₂O (2a) and [Fe(NH₂trz)₃](NO₃)₂ (3). These materials are one of the most extensively investigated SCO complexes due to their remarkable characteristics with thermal spin transitions often accompanied by wide hysteresis loops with transition temperatures located close to room temperature, leading to potential applications, e.g. as memory and sensor devices.⁵⁴ Recent work has demonstrated the additional fascinating effects of such 1D coordination polymers, such as the dependence of the spin transition on particle size⁵⁵ and mechanochemical recrystallisation,⁵⁶ as well as pronounced lattice softening of 1 upon replacing some of the Htrz ligands with NH₂trz.⁵⁷ These molecules contain the rigid bridging 1,2,4-triazole ring that distorts the geometry of the neighbours of a given spin if a centre is of different spin.⁴⁷ This effect leads to a particularly strong short-range interactions (see ref. 13) with a Γ value of more than 10³ cm⁻¹ (ref. 58) for 1.⁵⁹ Typically, in molecular crystals the direct elastic interactions between molecules does not exceed 10–100 cm⁻¹.¹³ Our previous results for 1 and [Fe(NH₂trz)₃]X₂ revealed that the DFT modelling of oligonuclear fragments provides a very good fit to the experimentally observed pDOS, both for the pure Fe(II) complex and for the Zn(II)-diluted ones.⁴¹

This report is organised as follows: firstly, we describe the modelling of four SCO materials – 1, 2, 2a (ref. 60 and 61) and 3.⁶² For each complex a heptanuclear model molecule was optimised with the B3LYP* functional. There are 20 different spin isomers with 5 switching Fe(II) centres. For each the electronic energy and the normal vibrations were calculated. The latter gave the vibrational contribution to entropy, S_vib. This allowed the calculations of the thermal dependence of the Gibbs free energy for all spin isomers. Further, we estimated the electronic energies of all spin isomers using other functionals (B3LYP and CAM-B3LYP) for the geometries obtained with B3LYP* or after optimisation. Then, for comparison, we repeated the optimisation and vibrational calculations for all systems using the B3LYP functional with dispersion correction. The analysis of the obtained electronic energies revealed that the relative (respective to that of the pure LS state) energy of each spin isomer is a sum of two factors: (a) functional dependent spin transition energy (b) strain energy due to presence of both LS and HS centres. The dependence of this energy on the applied functional is discussed. The temperature dependence of the Gibbs free energy for the B3LYP functional points towards ferroelastic interaction within the chain. The values of energy and entropy of the spin transition for 1, 2 and 2a were also determined experimentally by differential scanning calorimetry (DSC). The temperature dependence of ΔG was also derived based on the electronic spin transition energies derived from calorimetric measurements.

Results and discussion

Calculations of electronic energies and thermodynamic functions of different spin isomers of the heptanuclear models of 1, 2, 2a and 3

The electronic energies of 1, 2, 2a and 3 were calculated using the models described in the Experimental and computational methods section. Following the approach described in ref. 47, heptanuclear models were used in which the five inner Fe-atoms were coordinated with six 1,2,4-triazole-type ligands, while the terminal atoms were coordinated with three bridging ligands and three water molecules and were kept in the HS state for all spin isomers. The model of 1 is shown in Scheme 1. The optimised coordinates of all 20 spin isomers for the systems under study are given in the SI. The most important obtained structural feature is the shortening of the Fe–N bonds of the HS centres by their LS neighbours and the lengthening of the Fe–N bonds of the LS centre caused be its HS neighbours. This effect was previously modelled with DFT and detected with nuclear inelastic scattering.⁴⁷ It is also reflected in Table 1 which lists the mean Fe–N distances for the pure LS and HS models of the complexes under study compared to the spin isomers revealing the largest strain (H_strain vide infra). For further information see the list of all obtained atomic coordinates in SI.

Scheme 1 Graphic representation of the heptanuclear model used for calculations. The numbering of iron centres is given. The model molecule of 1 is shown.

Table 1 Mean Fe–N values for the LLLLL and HHHHH isomers of the complexes under study compared to those isomers revealing the highest strain (vide infra). Note the increase of the LS Fe(2)–N with HS neighbours bonds by ca. 0.2 Å compared to LS-Fe(2) ones with LS neighbours for 1, 2 and 2a and that of LS Fe(1) for LHLHL isomer of 3 (the corresponding pairs given in bold). The decrease of the Fe–N bonds for the HS centre upon introduction of the LS neighbours for 1, 2 and 2a varies from 0.2 Å for 1 to 0.32 Å for 2a. This shortening is 0.21 Å for Fe(2′) in 3 (the corresponding pairs marked blue). The distances for the negatively charge trz⁻ ligand in 1 are given in italics

		Fe(3)	Fe(2)	Fe(1)	Fe(2′)	Fe(3′)
1	LLLLL	1.995/1.999	1.985/1.989	1.979/1.989
	HHHHH	2.225/2.134	2.225/2.131
	LLHLL	1.992/1.996	2.000/1.998
2	LLLLL	1.985	1.978	1.992
	HHHHH	2.189	2.182
	HLHLH	2.166	1.997
2a	LLLLL	1.983	1.989	1.978
	HHHHH	2.167	2.180
	HLHLH	2.152	2.012
3	LLLLL	2.001	1.994	1.991	1.986	1.996
	HHHHH	2.189	2.177	2.186		2.182
	LHLHL	2.000	2.159	2.009		2.008

---

Before discussing the obtained energies and thermodynamic functions, we briefly review the intramolecular interactions between ligands and anions/solvated water. The relevant interaction patterns are discussed in detail in the SI (Fig. S5–S8). As pointed previously by Vela and Paulsen,¹⁸ the calculation of several spin isomers for different modifications of the second coordination sphere provides a wealth of structural data that will not be discussed in this paper. Here, we would like to indicate the most important patterns. The first interaction mode is the hydrogen bond between the anion or water and the C–H protons of the triazole ring. Typically, the anion or water oxygen forms a bridge between two such hydrogens of the neighbouring rings. The second one involves the hydrogen bonding between anion/water and 4-substituent of the triazole, i.e. either H or NH₂ group. The inspection of the obtained LS structures of the heptanuclear models reveals 24F⋯H contacts shorter than 3.0 Å for 1. The same distance criterion reveals 42C⋯H contacts for 2 and 25Cl⋯H and 30O⋯H contacts for 2a. The model of 3 reveals 56 contacts between ligand hydrogens and the nitrate oxygens.

This sort of interactions is likely to influence the energy of the spin transitions and affect the elasticity of the 1D chain. Notably, for 1 and 2 with only the F–H and Cl⁻–H interactions the LS → HS transition leads to the increase of the corresponding distances.

On the other hand, there is no clear effect of the spin state on the intermolecular interactions in 2a involving the O(H₂O)···HC(triazole) OH₂⋯Cl and NH₂⋯Cl interactions. For 3, the shortening of nitrate-O contacts with the triazole ring hydrogen is observed on HS to LS switching and the opposite effect occurs for nitrate-O-amine protons contacts.

With five inner centres that may be either in HS or LS state there are twenty possible spin isomers:¹

- Two with all inner Fe-centres being in either LS or HS state, that are further denoted as HHHHH and LLLLL, respectively, corresponding to x_HS of 1 and 0, respectively.

- Three with one HS centre and four LS ones (the HS defect in LS matrix, corresponding to L⁴H configuration, x_HS = 0.2), denoted as HLLLL, LHLLL and LLHLL.

- Three with one LS centre and four HS ones (the LS defect in HS matrix, corresponding to LH⁴ configuration, x_HS = 0.8), denoted as LHHHH, HLHHH and HHLHH.

- Six with three LS and two HS centres (corresponding to L³H² configuration, x_HS = 0.4), further named as “block” (LLLHH), “alternate” (LHLHL), “two pair” (LHHLL), “three” (HLLLH) and two “one pair” (LLHLH and HLLHL).

- Six with three HS and two LS centres (corresponding to L²H³ configuration, x_HS = 0.6), again named as “block” (LLHHH), “alternate” (HLHLH), “three” (LHHHL), “two pair” (HLLHH) and two “one pair” (HHLHL and LHHLH).

For each spin isomer the geometry optimisation, followed by normal modes calculation was performed with the B3LYP* functional⁷¹ and the CEP-31G basis set,⁶³ the combination used by us previously to model the 1D SCO chains under study.⁴⁷ The coordinates of all 20 model structures of spin isomers are listed in SI.

The obtained electronic energies calculated relative to those of the LLLLL models are shown in Tables 2, 3, 4 and 5 for 1, 2, 2a and 3, respectively.

Table 2 Calculated electronic energies, differences of zero-point vibrational energies and H_strain (in kJ mol⁻¹) parameters of all spin isomers of the heptanuclear model of 1 (E_el and ZPE assumed to be 0) and the derived values of H_coop and its analogue for Fe(2). The values of E_ad calculated and derived from the calorimetric results are given

	Eel				Hstrain			B3LYP D3/opt^c
	B3LYP*		B3LYP^a	CAM-B3LYP^a	B3LYP*	B3LYP^a	CAM-B3LYP^a
	Eel	ΔZPE
a Calculated for the geometry optimised with B3LYP* using B3LYP or CAM-B3LYP.b Hcoop calculated for the Fe(2) centre, equal to the difference between EHLLLLLH → EHLHLLLH and EHHLHHHH → EHHHHHH spin transition energies.c Optimised with B3LYP and D3 dispersion correction.d (Evib(HHHHH) − Evib(LLLLL))/5, i.e. the change of vibrational energy per 1 centre at the experimentally determined Tc.
LLHLL	55	−12	28	25	13	13	11	18
LHLLL	44	−11	17	15	2	2	1	8
HLLLL	40	−12	13	11.5	−2	−2	−3	4
HLLHL	83.5	−23	29	26.5	0	−1	−2	13
LHHLL	88	−23	35	31	4	5	3	16
LLHLH	94	−23	40	36	10	10	8	10
HLLLH	80	−24	26.5	23	−4	−4	−5	8
LHLHL	86	−22	33	29	2	3	1	16
LLLHH	74	−23	21	18.5	−10	−9	−10	4
HLLHH	114	−35	35	30	−11	−11	−12	8
HLHHL	127	−35	47	41.5	2	1	−1	13
HLHLH	119	−34	39	34	−6	−7	−8	12
LHHHL	121	−34	41	36	−4	−5	−6	13
LHLHH	117	−34	37.5	32.5	−8	−8	−10	11
LLHHH	118	−35	38	34	−7	−8	−8	8
HHLHH	176	−50	69.5	64.5	9	9	8	7
HLHHH	178	−50	73	68	11	12	12	6
LHHHH	180	−50	75	69	13	14	13	7
HHHHH	209	−62	76	70.5
Ead	42		15	14				16
E^calad		28.5			E^calad = ΔH^calHL (23.4) − ΔEvib (360 K)^d
Hcoop	22		21.5	19				25
Hcoop Fe(2)^b	13		14	12.5

---

Table 3 Calculated electronic energy differences, zero-point vibrational energies and H_strain (in kJ mol⁻¹) parameters of all spin isomers of the heptanuclear model of 2 (E_el and ZPE assumed to be 0) and the derived values of H_coop and its analogue for Fe(2). The values of E_ad calculated and derived from the calorimetric results are given

	Eel				Hstrain			B3LYP D3/opt^c
	B3LYP*		B3LYP^a	CAM-B3LYP^a	B3LYP*	B3LYP^a	CAM-B3LYP^a
	Eel	ΔZPE
a Calculated for the geometry optimised with B3LYP* using B3LYP or CAM-B3LYP.b Hcoop calculated for the Fe(2) centre, equal to the difference between EHLLLLLH → EHLHLLLH and EHHLHHHH → EHHHHHH spin transition energies.c Optimised with B3LYP and D3 dispersion correction.d (Evib(HHHHH) − Evib(LLLLL))/5, i.e. the change of vibrational energy per 1 centre at the experimentally determined Tc.
LLHLL	45.5	10	19	16	2.5	3	3	2
LHLLL	62	13	35	31.5	19	19	18	20
HLLLL	52.5	11	26	26	9.5	10	13	11
HLLHL	114	23	61	59	28	29	32	32
LHHLL	99	22	46	40.5	13	14	14	6
LLHLH	97	21	44	42	11	12	15	12
HLLLH	105	22	52	55	19	20	28	23
LHLHL	123	24	68.5	62.0	37	36.5	35	40
LLLHH	106	24	52	50	20	20	23	20
HLLHH	142	32	61	63	13	13	23	24
HLHHL	147	34	66.5	62	18	18.5	22	13
HLHLH	148	32	69	68	19	21	28	23
LHHHL	151	33	71	63	22	23	23	11
LHLHH	166	35	85	81	37	37	41	39
LLHHH	130	33	50	42	1	2	2	0
HHLHH	189	43	81	75	17	17	21	12
HLHHH	176	42	69	70	4	5	16	16
LHHHH	183	45	74	62	11	10	9	12
HHHHH	215	54	80	67
Ead	43		16	13
E^calad		13.5			E^calad = ΔH^calHL (8) − ^dΔEvib (340 K)
Hcoop	19.5		20	24				14
Hcoop Fe(2)^b	23		24	34.5

---

Table 4 Calculated electronic energies, zero-point vibrational energies and H_strain (in kJ mol⁻¹) parameters of all spin isomers of the heptanuclear model of 2a (E_HLLLLLH assumed to be 0) and the derived values of H_coop and its analogue for Fe(2). The values of E_ad calculated and derived from the calorimetric results are given

	Eel				Hstrain			B3LYP D3/opt^c
	B3LYP*		B3LYP^a	CAM-B3LYP^a	B3LYP*	B3LYP^a	CAM-B3LYP^a
	Eel	ΔZPE
a Calculated for the geometry optimised with B3LYP* using B3LYP or CAM-B3LYP.b Hcoop calculated for the Fe(2) centre, equal to the difference between EHLLLLLH → EHLHLLLH and EHHLHHHH → EHHHHHH spin transition energies.c Optimised with B3LYP and D3 dispersion correction.d (Evib(HHHHH) − Evib(LLLLL))/5, i.e. the change of vibrational energy per 1 centre at the experimentally determined Tc.
LLHLL	70	−14	42	38	18	17	13	17
LHLLL	60	−13	33	32	8	8	7	7
HLLLL	63	−12	36	37	11	11	12	11
HLLHL	122	−25	68	68	18	18	18	18
LHHLL	115	−26	61	57	12	11	7	9
LLHLH	132	−26	76	74	28	26	23	28
HLLLH	125	−24	71	73	21	21	23	22
LHLHL	117	−24	62	63	13	13	13	13
LLLHH	110	−25	56	58	6	6	7	6
HLLHH	172	−37	91	94	17	17	18	16
HLHHL	176	−38	95	92	21	20	17	20
HLHLH	193	−38	111	109	37	36	34	40
LHHHL	170	−35	91	94	15	16	18	2
LHLHH	167	−36	86	89	12	11	13	11
LLHHH	165	−38	84	82	10	9	7	9
HHLHH	218	−48	110	114	10	10	13	9
HLHHH	227	−50	118	117	19	18	17	19
LHHHH	209	−50	101	100	2	1	0	0
HHHHH	260	−62	125	126
Ead	52		25	25				34
E^calad		21			E^calad = ΔH^calHL (15.5) − ^dΔEvib (330 K)
Hcoop	28		21.5	19				25
Hcoop Fe(2)^b	27		14	12.5

---

Table 5 Calculated electronic energies, zero-point vibrational energies and H_strain (in kJ mol⁻¹) parameters of all spin isomers of the heptanuclear model of 3 (E_HLLLLLH assumed to be 0) and the derived values of H_coop and its analogue for Fe(2). The values of E_ad calculated and derived from the calorimetric results are given

	Eel				Hstrain			B3LYP D3/opt^c
	B3LYP*		B3LYP^a	CAM-B3LYP^a	B3LYP*	B3LYP^a	CAM-B3LYP^a
	Eel	ΔZPE
a Calculated for the geometry optimised with B3LYP* using B3LYP or CAM-B3LYP.b Hcoop calculated for the Fe(2) centre, equal to the difference between EHLLLLLH → EHLHLLLH and EHHLHHHH → EHHHHHH spin transition energies.c Optimised with B3LYP and D3 dispersion correction.d (Evib(HHHHH) − Evib(LLLLL))/5, i.e. the change of vibrational energy per 1 centre at the experimentally determined Tc.
LLHLL	52	−11	25	25	15	17	17	11
LHLLL	65	−10	37	31	27	22	22	5
HLLLL	41	−11	20	17	4	12	9	−7
HLLHL	93	−21	39	33	19	16	16	9
LHHLL	115	−22	61	55	40	38	38	5
LLHLH	85	−20	33	33	11	16	16	2
HLLLH	81	−22	27	25	6	8	8	−10
LHLHL	118	−18	64	55	44	38	38	25
LLLHH	84	−21	30	25	10	8	8	4
HLLHH	116	−31	35	30	5	4	4	18
HLHHL	131	−31	53	52	20	26	26	10
HLHLH	117	−30	39	35	6	10	10	−3
LHHHL	148	−31	67	58	36	32	32	10
LHLHH	130	28	48	42	18	17	17	35
LLHHH	139	−32	42	42	10	16	16	20
HHLHH	166	−42	58	49	18	15	15	8
HLHHH	161	−41	55	46	13	12	12	−10
LHHHH	172	−44	64	52	24	18	18	−2
HHHHH	185	−52	52	43
Ead	37		8.5	10				18
E^calad		27			E^calad = ΔH^calHL (23.1) − ^dΔEvib (330 K)
Hcoop	33		31	31				33
Hcoop Fe(2)^b	41		39	34

---

The comparison of the so defined relative energies of the spin isomers allows the following conclusions (see Tables 2–5 and Fig. S1–S4, for the E_el values):

(i) The electronic energy increases with the number of HS centres in the inner five-nuclear Fe-fragment of the heptanuclear model.

(ii) Each spin isomer containing the same ratio of HS to LS centres reveals a different electronic energy. There is hardly a pattern of stabilisation/destabilisation of a given spin isomer that is common for all systems. For example, while both “block” spin isomers LLHHH and HHHLL are of lowest energy of all L²H³ isomers (for 1, 2a and 3, respectively, see Tables 2, 4 and 5), this is not the case for 2 (cf. Table 3). The “block” LLLHH isomers are of the lowest energy for 2 and 2a but not for 1 and 3. Similarly, the “alternate” antiferroelastic LHLHL is the most destabilised among L³H² of 1 and 3 (Tables 2 and 5), but not for 2 and 2a (Tables 3 and 4).

(iii) The largest energy gap occurs between the LLLLL isomer and the one with one HS defect (LLHLL, LHLLL and HLLLL). The minimal values are 40, 45.5, 60 and 41 kJ mol⁻¹ for 1, 2, 2a and 3, respectively. This implies a high stabilisation of the LS state in the pure LS matrix. Still, there is no pattern of a particular stabilisation or destabilisation of the particular spin isomer with one HS defect. The centrosymmetric LLHLL isomer is of lowest energy for 1 and 3 (see Tables 2 and 5) while it is of highest energy for 2 and 2a (see Tables 3 and 4).

In this respect complex 2 is particular, revealing a significant energy gap between the L²H³ and LH⁴ spin isomers. The smallest energy gap is that between the LHLHH and LLHHH ones (36 kJ mol⁻¹).

(iv) The calculated electronic energies of LLLLL → HHHHH transitions per one switching centre are 42, 43, 52 and 37 kJ mol⁻¹ for 1, 2, 2a and 3, respectively. The calorimetrically obtained values (see Table 7) are respectively 23.4, 9, 12 and 23 kJ mol⁻¹. The zero-point corrected electronic energies per one center (Tables 2–5) are 29.4, 32.2, 39.6 and 26.6 kJ mol⁻¹. Thus, the calculated LLLLL → HHHHH transition energies for the 1D models of the system under study do not fit to the experimentally observed trend of the spin transition energies increasing in the 3 < 1 < 2 < 2a sequence.

Table 6 Calculated electronic energies, and H_strain (in kJ mol⁻¹) parameters of all spin isomers of heptanuclear model of 1 (E_HLLLLLH assumed to be 0) and the derived values of obtained with B3LYP and CAM-B3LYP functionals with and without dispersion correction, and with PBEh functional

	B3LYP/D3		B3LYP		CAM-B3LYP/D3		CAM-B3LYP		PBEh
	Eel	Hstrain	Eel	Hstrain	Eel	Hstrain	Eel	Hstrain	Eel	Hstrain
LLHLL	34	18	28	19	29	18	26	18	−2	19
LHLLL	24	8	18	8	20	8	15	7	−10	8
HLLLL	20	4	14	4	16	4	12	4	−17	4
HLLHL	44	13	31	12	35	12	28	11	−33	9
LHHLL	47	16	36	17	39	16	32	16	−26	16
LLHLH	42	10	29	10	33	8	26	9	−31	11
HLLLH	40	8	27	8	31	8	24	8	−33	9
LHLHL	47	16	34	15	38	15	30	14	−28	14
LLLHH	35	4	22	3	27	4	19	3	−39	3
HLLHH	56	8	36	8	43	7	32	7	−55	8
HLHHL	61	13	42	14	48	13	37	13	−49	14
HLHLH	60	12	41	12	46	11	35	11	−51	12
LHHHL	61	13	43	14	48	13	37	13	−50	13
LHLHH	59	11	39	10	46	11	34	9	−53	10
LLHHH	55	8	36	8	43	8	32	7	−51	12
HHLHH	70	7	44	6	53	6	38	5	−78	6
HLHHH	70	6	45	6	53	6	38	6	−77	7
LHHHH	71	7	46	8	54	7	40	7	−76	8
HHHHH	79		48		59		41		−105
Ead	16		10		12		8		−21
Hcoop	25		24		24		23		25

---

Table 7 Phase transition temperatures and thermodynamic parameters as deduced from DSC measurements and calculated (B3LYP*) vibrational entropies

1D chain	#^b	Tmax^↑ [K]	Tmax^↓ [K]	ΔT [K]	ss%^a	ΔHHL [kJ mol^−1]	ΔSHL [J mol^−1 K^−1]	ΔSvib (DFT)
a ss: switching sites determined by ^57Fe Mössbauer spectroscopy when comparing the highest and lowest temperature spectra.b Run numbering.c This work.d Ref. 69.e Ref. 64a.f 100% spin conversion assumed on the basis of magnetic and NIS data (see Experimental section for discussion).
[Fe(Htrz)2trz]BF4 (1)	1^c	387	343	44	86	23.4	64.1	50.7
								73 (390 K)
								72 (340 K)
	2^c	390	347	43	86	23.4	63.5	50.2
	3^c	391	344	47	86	24	65.3	52
[Fe(NH2trz)3]Cl2 (2)	1^c	343	335	8	^f	8.93	26.06	12.73
								52 (310 K)
[Fe(NH2trz)3]Cl2·2.5H2O (2a)	1^c	332	311	21	^f	12.22	39.38	23.05
								60 (330 K)
[Fe(NH2trz)3](NO3)2 (3)	1^d	347	314	33	92	23(1)	69.6(1)	56.2
								70 (350 K)
								56 (310 K)
[Fe(NH2trz)3]TiF6·0.5H2O	1^e	210	200	10	72	7.1	34.5	21.1
[Fe(NH2trz)3]TiF6·H2O	1^e	206	201	5	75	7.3	36	22.6
[Fe(NH2trz)3]ZrF6·0.5H2O	1^e	233	209	24	85	7.4	33.6	20.2
[Fe(NH2trz)3]ZrF6	1^e	256	220	36	85	6.7	30.5	17.1
[Fe(NH2trz)3]SnF6·0.5H2O	1^e	224	203	21	75	6.3	31.2	17.8
[Fe(NH2trz)3]SnF6·H2O	1^e	244	224	20	80	8	34.3	20.9
[Fe(NH2trz)3]TaF7·3H2O	1^e	206	197	9	72	5.9	30	17
[Fe(NH2trz)3]GeF6·H2O	1^e	213	211	2	62	6.4	30	17
[Fe(NH2trz)3]GeF6·0.5H2O	1^e	213	207	6	58	6.9	32	19

---

(v) Finally, it is of interest how much the electronic energy of the spin isomers containing the same number of HS and LS centres depends on the distribution of the centres of a given spin. For example, what is the difference between L⁴H spin isomer LLHLL, LHLLL and HLLLH. The inspection of the Tables 2–5 reveals that the largest difference of 37 kJ mol⁻¹ is obtained for LHLHL and HLLLH spin isomers of 3. For the L²H³ isomers of 2, the energy difference between LHLHH and LLHHH is 36 kJ mol⁻¹. These values are derived for the systems of the same multiplicity and are supposedly independent from the applied functional (vide infra). Consequently, one may state that the energy differences between the different pattern of spin distribution for a given number of HS and LS centres are comparable or larger than the typical spin transition energies that span the range of 2–20 kJ mol⁻¹ as observed for 1D triazole-based complexes.⁶⁴ Note that zero-point vibrational energies effectively do not depend on the spin distribution in spin isomers containing the same number of HS and LS centres (see Tables 2–5). Moreover, with the increasing number of the HS centres the ZPE decreases linearly with the number of the HS centres. For example, for 1 the ZPE is lower at −11/-12 kJ mol⁻¹ for L⁴H isomers compared to the fully LS model. The corresponding values for L³H², L²H³, LH⁴ and full HS spin isomers are respectively −23, −34 to −35, −50 and −62 kJ mol⁻¹.

In the next step, the cooperativity parameter H_coop was derived for the central Fe(1) and for Fe(2) (Scheme 1). As defined in ref. 50, H_coop corresponds to the difference of calculated spin transition energies within the matrix of LS and HS neighbours. Hence, H_coop for Fe(1) is equal to [(E(LLHLL) − E(LLLLL)) − (E(HHHHH) − E(HHLHH))], while H_coop for Fe(2) is equal to [(E(LHLLL) − E(LLLLL)) − (E(HHHHH) − E(HLHHH))]. H_coop for Fe(1) results to 22, 19.5, 28 and 33 kJ mol⁻¹ for 1, 2, 2a and 3 respectively. For Fe(2) the corresponding values are 13, 23, 27 and 41 kJ mol⁻¹ (Tables 2–5). As stated above (ref. 50) these values reflect the amount of the deformation of the FeN₆-core by the matrix of a different spin for the LS and HS centres. Another parameter is the stress due to neighbourhood of both LS and HS centre, being the net result of deformation (strain) of all centres from their regular HS and LS geometry. To derive it, one needs to assume that the calculated electronic energy E_el of a given spin isomer relative to the LLLLL one is the sum of two contributions: (i) the electronic energy (dependent on the exchange–correlation functional) of a single-centre spin transition (E_ad). The latter is independent on the position of the centre in the chain. (ii) The stress due to geometric strain in a given spin isomer. Hence, the relative electronic energy of a given isomer containing nHS centres and 5-nLS centres is:

E_el(H_nL_5−n) = nE_ad + H_strain (1)

where H_strain is the cumulative stress due to the deformation of all HS and all LS centres from their geometries in the pure HS and LS phases, respectively. It becomes zero for the pure LS and HS phases.

To derive H_strain for our models, we estimated E_ad to be 1/5 of the energy difference between LLLLL and HHHHH spin isomers, assuming a negligible stress in the heptanuclear model molecule for the LS and the HS state (i.e. for the LLLLL and HHHHH systems). We obtain 42, 43, 52 and 37 kJ mol⁻¹ for 1, 2, 2a and 3 respectively. With these numbers we derive H_strain for all spin isomers that are listed in Tables 2–5. The NH₂trz complexes 2, 2a and 3 display a similar pattern of H_strain being the largest for the L³H² and L²H³ spin isomers with maximal values of +44 kJ mol⁻¹ for the “alternate” LHLHL spin isomer of 3. The other spin isomers L³H² and L²H³ display a very low H_strain with that for the “block” LLHHH isomer of 2 being only 1 kJ mol⁻¹. On the other hand, the H_strain values obtained for 1 reveal an opposite trend. Particularly for L²H³ and L³H² spin isomers H_strain is negative, indicating that strain induces a spin stabilizing effect. The structural reason for this effect is not immediately clear, yet with B3LYP/D3 optimisation only the positive values of H_strain are obtained (see Table 2).

It is of interest whether the values of H_strain are reproducible using other functionals since any functional chosen leads to different absolute electronic energy values.¹⁵ Therefore, for all structures optimised with B3LYP*, the calculations of the electronic energies were performed also with CAM-B3LYP⁶⁵ and B3LYP⁶⁶ using the CEP-31G basis set (Tables 2–5). The calculated values of E_ad decrease from 42–43 kJ mol⁻¹ to 16–13 kJ mol⁻¹ on changing from B3LYP* to B3LYP/CAM-B3LYP for 1 and 2, from 52 to 25 kJ mol⁻¹ for 2a and from 37 to 10 and 8.5 kJ mol⁻¹, respectively for 3. Nevertheless, the energy sequence as well as H_coop and H_strain are comparable for all three functionals. There are some exceptions for 2 calculated with CAM-B3LYP, in line with previous findings⁵¹ but H_coop and H_strain of 2 are equivalent. Assuming a constant value of E_ad the cooperativity parameter H_coop can by expressed for the central Fe(1) as H_coop (Fe(1)) = H_strain(LLHLL) − H_strain(HHLHH). Importantly, for the LLLLL → LLHLL transition the strain is introduced, while for the HHLHH → HHHHH one it is released. Noteworthy, the energy of elastic amplification observed after photo-induced LS → HS transition⁶⁷ can be also expressed in term of H_coop and H_strain.²

Based on the electronic energies of spin isomers collected in Tables 2–5 several other parameters characterizing the short-range interactions may be derived (see Table S1). These results are presented in SI. By and large, the following conclusions may be drawn:

(i) The presence of HS neighbours generally favours the LS to HS transition four all four modelled systems.

(ii) This effect weakens with the distances between the HS centre and the switching centre only when the HS centre can be considered as a defect in a LS matrix. For example, the LLLLL to HLLLL and HLLLL to HLLLH transition display practically no difference for all modelled molecules, while the LHHHL to LHHHH and LHHHH to HHHHH transition energies show large differences (e.g. from −12 to 30 kJ mol⁻¹ for 2a and 1, respectively, 29 and 41, and 59 and 29 kJ mol⁻¹). The exception is 2, for which these two transitions reveal the same energies. This suggests that the elastic effect of the defect has a higher range in the more elastic matrix of the HS centres. This shows that the short-range interactions affect centres at a distance of ca. 1.6 nm in the HS matrix. On the other hand, there is no difference in the calculated spin transition energies of LLLLL → LLLLH and LLLLH → HLLLH transitions.

(iii) The obtained parameters (see Table S1) for 2 are somewhat lower than for the other systems while 3 is exceptional in showing practically no zero values. Considering that NH₂trz–Cl (with no additional ligand–water interaction) contacts are supposed to be the weakest while those of NO₃–NH₂trz type are supposed to be the strongest, we conclude that the strong anion–ligand interactions may increase the transfer of the elastic distortion.

In next step, we investigated how H_strain differs if the optimisations and energy calculations are performed with other functionals. In this case we carried out the modelling with B3LYP functional including dispersion correction. The obtained values of H_strain and H_coop are shown in Tables 2–5 the calculated electronic energies for 1 given in Table 6, while those for 2, 2a and 3 are given in, together with the ZPE corrections SI (Tables S3–S5). The values of ZPE corrections for 1 obtained with different functionals are collected in Table S6.

The comparison of H_strain obtained for the structures optimised with B3LYP* and B3LYP/D3 reveals a surprising pattern (see Tables 2–5). Nearly the same values were derived for 2a (with exception of the LHHHL isomer). Fairly similar values were found for 2, the most significant differences found for LHHHL (11 and 22 kJ mol⁻¹ for B3LYP/D3 and B3LYP*, respectively) and HLHHH (16 and 4 kJ mol⁻¹, respectively). Noteworthy, the latter two values obtained with B3LYP/D3 are very close to those obtained with CAM-B3LYP for the geometry obtained with B3LYP* (see Table 3). On the other hand, the H_strain parameters derived with B3LYP/D3 for 1, show a different pattern than those derived with B3LYP*, particularly for the L³H² and L²H³ spin isomers where the negative values obtained with B3LYP* turn positive for B3LYP/D3 (see Table 2). For 3 this discrepancy is also observed, again with three values of H_strain revealing the negative values for B3LYP/D3 (see Table 5). On the other hand, the values of H_coop seem to be closer for the two functionals yielding 22 and 25 kJ mol⁻¹ for B3LYP* and B3LYP/D3 for 1, and respectively 19.5 and 14 kJ mol⁻¹ for 2, 28 and 25 kJ mol⁻¹ for 2a and 33 kJ mol⁻¹ for both functionals for 3. That means that the H_strain parameters for LLHLL and HHLHH are less dependent on the choice of the functional.

To get a further insight into this dependence we performed optimisations with B3LYP and CAM-B3LYP with and without dispersion correction and additionally for PBEh functional for 1 as an example of another hybrid functional, not related to B3LYP. The results are given in Table 6. Additionally, the CAM-B3LYP/D3 and PBEh modelling was done for 2, the results are shown in Table S3. The results show that within these functionals quite consistent values of H_strain are obtained. The largest differences of H_strain are obtained for the LLHHH spin isomer, with all B3LYP and CAM-B3LYP (with or without dispersion) functionals yielding the values of 7–8 kJ mol⁻¹ while PBEh yielded 12 kJ mol⁻¹. The differences between results with and without dispersion correction are negligible. Again, the obtained H_coop values lie all in the range 23–25 kJ mol⁻¹. Further calculations of H_coop with other functionals (TPSS and TPSSh) show quite constant values for each modelled system (see Table S2). The data obtained for 2 with B3LYP/D3, CAM-B3LYP/D3 and PBEh are also quite close (see Table S3), the exceptions being “two pair” LHHLL and HLLHH spin isomers (H_strain for LHHLL equal to 6, 4 and 13 kJ mol⁻¹) for B3LYP/D3, CAM-B3LYP/D3 and PBEh, respectively, the corresponding values for HLLHH being 24, 5 and 3 kJ mol⁻¹. To summarize, the values of H_strain for different spin isomers are comparable if for the same geometry (obtained with B3LYP*) the energies are calculated with B3LYP and CAM-B3LYP. Whether that obtained with B3LYP* (optimisation and energy calculations) differ from that obtained for the optimisation and energy calculations with B3LYP, CAM-B3LYP or PBEh depends on the model molecule. Within the last series of functionals they match well. The reason for the different values obtained with B3LYP* is not immediately clear, yet the comparison of the selected Fe–N distances suggests that the obtained values of H_strain correlate with the ability of different functionals to give the HS Fe–N bonds lengths that differ less than 0.01 Å (see Table S7).

Thermodynamics of spin transition

The temperature dependencies of the Gibbs free energy changes for the LLLLL → L_5−nH_n spin transitions have been computed for the 19 spin isomers. Our calculations are based on the approach described in ref. 42, i.e. according to the following formula:

ΔG(T) = ΔE_el + ΔH_vib(T) − T (ΔS_vib(T) + ΔS_mag + ΔS_mix + ΔS_conf) (2)

where E_el denotes the calculated electronic energies, n the number of HS centres in the inner chain of the five irons, ΔS_mag = R[ln(2S + 1)_LSⁿ_HS⁽⁵⁻ⁿ⁾/ln(2S + 1)_LS⁵] (ref. 1a). The latter is the entropy change due to change of multiplicity on going from the pure LS state to a given spin isomer LⁿHS⁵⁻ⁿ. In simple form, it is equal to nR ln 5. ΔS_mix is the mixing entropy change given by the formula: ΔS_mix = −R(x_HS ln x_HS + (1 − x_H)ln(1 − x_HS)). Note that x_HS in our case is 0.2 for L⁴H, 0.4 for L³H², 0.6 for L²H³ and 0.8 for LH⁴ spin isomers. ΔS_conf reflects the probability of formation of a given isomer. It is zero for the centrosymmetric isomers and R ln 2 for those bearing no inversion centre. It may be given by the formula ΔS_conf = R ln(nk), were nk denotes the number of molecular permutations corresponding to a given spin isomer. For example, there are two ways the isomer HLLLL could be formed, but only one for LLHLL one. E_vib(T) is the vibronic energy and ΔS_vib(T) denotes the change of vibronic entropy upon the spin transition. Both, E_vib(T) and ΔS_vib(T) have been obtained by DFT based normal mode calculations.

Additionally, the ΔG(T) dependencies were also obtained using the E_el values calculated with B3LYP and CAM-B3LYP and all other parameters obtained with B3LYP* (see further discussion). For these calculations the geometries optimised with B3LYP* were used for the point-energy calculations of energies (followed by the stability check, vide infra). In this way the only parameter that is changed is the electronic energy characteristic for a given functional.

The plots of ΔG(T) calculated with the B3LYP* functional of all spin isomers of 1, 2, 2a, and 3 are shown in Fig. 1. The respective calculated plots of vibrational entropy obtained for B3LYP* functional are shown in Fig. 2.

Fig. 1 Calculated temperature dependence of ΔG of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear models of 1, 2, 2a and 3 obtained with B3LYP* functional and CEP-31G basis set.

Fig. 2 Calculated temperature dependence of ΔS_vib of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear model of 1, 2, 2a and 3 obtained with B3LYP* functional.

There are distinct values of Gibbs free energy and vibration entropy for each spin isomer of a given multiplicity. This means that electronic and vibrational energy, Gibbs free energy and entropy are dependent on the distribution of the LS and HS centres within the pentanuclear fragment of the chain. Accordingly, the different distribution leads to different level of strain extorted by the proximity of the neighbours of different spin. The so obtained ΔG(T) curves show that the combination of strain and entropic effects leads to a situation where the thermodynamic stability of a given isomer is not linearly dependent on the number of HS centres.

This is particularly pronounced for 2, where the LLHHL and LHLHL isomers are up to 15 kJ mol⁻¹ less stable than the LLHHH and HHLLH isomers below 100 K (Fig. 2).

According to Fig. 2, the entropy changes with respect to the LLLLL isomer scales with the number of HS centres. Somewhat less obvious, however, is the observation that the higher the multiplicity, the stronger is the dependence of the vibrational entropy on temperature. While for all systems, particularly for 1, ΔS_vib for model molecules with three HS centres reaches (blue lines in Fig. 4) its maximal value at 250–300 K. For the systems with more than three HS centres ΔS_vib still grows at higher temperatures: The calculated value of ΔS_vib for the HHHHH isomer of 2 grows at ca. 20 J K⁻¹ mol⁻¹ between 300 and 400 K. For 1, the vibrational entropy curves for spin isomers containing the same number of HS and LS centres are very close to each other. Those for 2, 2a and 3 deviate more.

Electronic energy tuning. ΔG(T) calculated with B3LYP and CAM-B3LYP for the geometries and entropy obtained with B3LYP

Actually, the obtained temperature dependencies of ΔG, calculated with the B3LYP* functional do not reveal the properties of the model systems, namely the abrupt spin transition, i.e. switching between fully HS and fully LS chains. Instead, differences of the Gibbs free energy values of all possible spin isomers tend to decrease with temperature. The fully LS system is predicted to have the lowest Gibbs free energy up to 400 K. This effect is likely to be due to the propensity of the B3LYP* functional to both lower the energy of the low-spin state and higher for the HS one, thus increasing the spin gap ΔE_HL = E_HS − E_LS (E_el in notation used here).¹⁵ Generally, the computed electronic energies of spin transition are highly dependent on correlation-exchange functionals.^14,15

On the other hand, the differences seem not to be that large in what concerns the vibrational properties and hence the vibrational energy and entropy. Note the importance of the accuracy of calculating ΔS_vib.⁶⁸ For the critical temperature ΔS_vib = 33 ± 13 J K⁻¹ mol⁻¹ was obtained for [Fe(phen)₂(NCS)₂] using different basis sets. This issue is discussed in detail in SI where we present the results of vibrational entropy calculation obtained with the B3LYP/D3 optimised structures of all spin isomers of 1–3. For 1, the results used with CAM-B3LYP/D3 and CAM-B3LYP and B3LYP with no dispersion correction are given. For 2 the results with CAM-B3LYP/D3 are also given. All calculated harmonic frequencies are given together with xyz coordinates in SI.

Thus, we first performed the calculations of the ΔG(T) dependencies using the vibrational entropy and energy values obtained with B3LYP*, while the electronic energies were computed with the B3LYP functional. The latter is known to give lower energies for the HS isomers of Fe(II). With CAM-B3LYP a further fine-tuning of the calculated electronic energy is provided. In each case, the geometries obtained with B3LYP* were applied. In the next step, ΔG(T) curves are given for the models optimised with B3LYP/D3.

We have shown above that the differences of the electronic energies within the series of spin isomers of the same multiplicity, that reflect the strain effects leading to intramolecular cooperativity are overall quantitatively retained, independently on the applied functional based on the B3LYP one. Yet, the relative stability of the isomers containing a given number of HS centres differs according to the relation CAM-B3LYP ≥ B3LYP > B3LYP*. This gives rise to a higher weight of entropic effects in the predicted values of ΔG.

The plots of ΔG(T) calculated with B3LYP and CAM-B3LYP of all spin isomers of 1, 2, 2a and 3 are shown in Fig. 3 and 4.

Fig. 3 Calculated temperature dependence of ΔG of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear model of 1, 2, 2a and 3 obtained with CAM-B3LYP functional, using the geometries and vibrational entropies and energies obtained for the B3LYP* functional.

Fig. 4 Calculated temperature dependence of ΔG of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear model of 1, 2, 2a and 3 obtained with B3LYP functional, using the geometries and vibrational entropies and energies obtained for the B3LYP* functional.

The inspection of respective figures allows the following conclusions: (i) with the lower value of E_ad the entropic effects may lead either to (i) a predicted change from the LLLLL ground state to the HHHHH spin isomer for 2, (ii) to HHHHH being the ground state in case of 3, (iii) a predicted change from the LLLLL ground state to the HHHLL isomer for 2a, and (iv) switching between HHHLL and HHHHH isomers in case of 1. For both 1 and 2, both the B3LYP and the CAM-B3LYP functionals predict the purely HS state as the ground state at temperatures above 150–200 K. Fig. 3 and 4 show that the entropy lowers the Gibbs free energy of the pure HS phase with increasing temperature. For 1, for which the strain could have negative values (Table 2) the same effect is obtained while ΔS_vib is significantly larger than for other complexes under study (Fig. 2).

(ii) With decreasing E_ad the strain effects may lead to a higher stabilisation of spin isomers having a majority number of HS centres. This is well seen for 2 for which the L²H³ spin isomers are predicted to be more stable than the L³H² ones when E_ad is calculated with B3LYP and CAM-B3LYP functionals. For 3, this is seen even with B3LYP*. In all cases this effect occurs at low temperatures at which the entropic effects are small.

In the next step, we calculated the ΔG(T) dependencies for all four systems under study using the previously mentioned calculations using B3LYP functional with the dispersion correction D3. Thus, we can compare the results obtained with B3LYP* with those obtained for the B3LYP functional giving another values of E_ad with the vibronic properties calculated also with B3LYP. The calculated ΔG(T) curves for the model molecules optimised with this method are given in Fig. 5.

Fig. 5 Calculated temperature dependence of ΔG of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear model of 1, 2, 2a and 3 obtained with B3LYP/D3 functional.

The comparison of the results shown in Fig. 5 with those in Fig. 4, i.e. these obtained with B3LYP/D3 with those calculated with B3LYP* for the geometry obtained with B3LYP* reveals that the ΔG(T) curves for 1 and 2 seems similar with some shift of the temperature at which the HHHHH turns to be ground state to higher temperature for B3LYP/D3. For 2a and 3 the differences are more pronounced. This effect may be related to the different values of E_ad within the two approaches. The data in Tables 2–5 show that E_ad are quite close for 1 (15 and 16 kJ mol⁻¹ for B3LYP* and B3LYP/D3 geometries, respectively) and for 2 (16 and 17 kJ mol⁻¹, respectively). For 2a and 3 they are respectively 25 and 34 kJ mol⁻¹, and 8.5 and 18 kJ mol⁻¹. The similar pattern was obtained for the optimisation with B3LYP with no dispersion correction and CAM-B3LYP (with and without dispersion correction) for 1 and for CAM-B3LYP/D3 (see SI). In summary, our study shows that for the two modelled systems (1 and 2) E_ad values obtained with B3LYP and CAM-B3LYP functionals (and geometries from B3LYP*) all HS state turns to be the ground state above 200 K. The B3LYP/D3 approach predicts all HS state to be the ground state of 1 above 100 K and over 200 K for 2. With the B3LYP* geometry, for 2 the pure LS state turns to be the ground state up to 100 K with both CAM-B3LYP and B3LYP obtained energies for 3, the HS ground state is predicted in the whole temperature range. The B3LYP/D3 calculations for 3 predict the all HS state over 300 K.

Thus, our approach suggests that the short-range intra-chain interactions lead to a ferroelastic character of the 1D SCO chains investigated in this study. An experimental determination of E_ad is possible with calorimetry. Therefore, we decided to determine experimentally SCO transition enthalpies in order to model the SCO behaviour with model molecules of 1–3.

Calorimetric data as input for DFT modelling

The temperature dependence of the heat capacity, C_p, of 1, 2 and 2a was recorded by DSC. (see chapter Thermal analysis and Experimental and computational methods for details) to derive phase transition temperatures on warming and cooling, hysteresis width, both energy and entropy associated to the spin transition, ΔH_HL and ΔS_HL as well as ΔS_vib (Table 7). Reported data fall in the range of previously reported [Fe(NH₂trz)₃]X₂ complexes with different anions, including 3 (Table 7).

The obtained ΔH_HL values for 1, 2 and 2a and for 3 (ref. 69) provide the sum of the electronic energy (i.e. E_ad), vibrational energy and the intermolecular energy of the spin transition. We assume that for each spin isomer ΔG regarding the LLLLL state per one Fe centre is given by:

ΔG(T) = n(ΔH_HL − ΔH_vib(T_c)/5) + ΔH_vib(T) + H_strain − T(ΔS_vib(T) + ΔS_mag + ΔS_mix + ΔS_conf) (3)

where n is the number of the HS centres within the pentanuclear chain. ΔS_vib and H_strain denote the vibrational entropy and interaction parameters, respectively (the latter listed in Table 2) obtained with B3LYP*. ΔH_HL is the calorimetrically determined energy of the all-LS to the all-HS transition. Note that on a contrary to eqn (2) the first term of eqn (3) corresponds to the sum of the electronic and vibronic energies. The former is independent of temperature. An additional contribution to ΔG comes from stress for those spin isomers containing the centres of both spin states. Thus, E_ad can be determined by summing up calorimetrical ΔH_HL and DFT calculated ΔH_vib(T_c).

For the LS → HS transition ΔH_vib is negative as the energy of the stretching vibrations is lowered upon this transition. Therefore, E_ad is larger than ΔH_HL. The so estimated values of E_ad are listed in Tables 2–5. For all systems our estimated E_ad values are significantly lower than the ones calculated with B3LYP*, while for the B3LYP/D3 optimised structures the relation between these values varies between the systems. Thus, for 2 the calorimetrically and B3LYP/D3 calculated values of E_ad are close to each other (13.5 vs. 16 and 17 kJ mol⁻¹). For 1, and 3, the values derived from calorimetric data are 9–15.5 kJ mol⁻¹ higher than those calculated with the above two functionals. For 2a, on the other hand the calculated E_ad of 34 kJ mol⁻¹ is higher than that determined calorimetrically (21 kJ mol⁻¹). The so derived ΔG(T) curves are shown in Fig. 6:

Fig. 6 Calculated temperature dependence of ΔG of the spin transition from the LLLLL spin isomers to all spin isomers of heptanuclear model of 1, 2, 2a and 3 obtained with the E_ad value estimated from the calorimetric data and H_strain for each spin isomer (see eqn (3)), using the geometries, H_strain values and vibrational entropies and energies obtained for the B3LYP* and B3LYP/D3.

(i) For 1, calculated dependencies point towards the full LS isomer to be the ground state up to ca. 20–210 K using the vibrational data obtained with both B3LYP* and B3LYP/D3. Around this temperature a crossing point of ΔG for all spin isomers is seen full HS isomer becomes that of the lowest Gibbs free energy at 315 K and above 400 K with B3LYP* and B3LYP/D3 vibrational data, respectively.

(ii) For 2, the fully LS isomer is predicted to be the ground state up to about 100 K with the vibrational data obtained with B3LYP*. Then the fully HS turns to be the one with lowest Gibbs free energy. The B3LYP/D3 predicts this to happen at ca. 65 K. For 2a, similarly full low spin switching temperatures of 65 and 85 K, respectively are predicted.

(iii) For 3, full HS isomer is predicted to be the ground state at 350 K if calculated with B3LYP/D3 derived vibrational data and at 350 K for the B3LYP* derived ones.

In all cases, this composite approach predicts the fully LS and HS chains to be the ground state depending on temperature. Thus our combined experimental and theoretical approach confirms that the heptanuclear 1D chain fragments of Fe(II) complexes with 1,2,4-triazole-based ligands have a ferroelastic character and switch between the full LS and the full HS state on temperature change. This ferroelasticity is a result of a strong increase of the entropic effects for the HS chains with increasing temperature and the additional strain occurring for the spin isomer containing both spin centres.

The above results allow the estimation of the effect of how the temperature dependencies of the free enthalpies of a spin isomer depend on the adiabatic electronic energies of the spin transition. Fig. 7 shows the dependencies obtained for 1 using eqn (3) with E_ad = ΔH_HL − ΔH_vib(T_c) as parameter. The diagrams show a very strong dependence of the predicted ground state on E_ad. For example, the change from 20 kJ mol⁻¹ to 10 kJ mol⁻¹ shift the temperature at the which the fully HS state exhibits the lowest Gibbs free energy from ca. 320 K to ca. 270 K. Note that the difference of 10 kJ mol⁻¹ may be compared to 12 kJ mol⁻¹ difference of the conformer energies for ethane,⁷⁰ thus corresponding to a very subtle changes of the molecular geometry. This effect is shown more clearly in Fig. 8 in which the relative free enthalpies of the full HHHHH model systems are shown as function of the temperature for different values of E. Last, but not least, the obtained temperature dependencies of ΔG obtained both on the basis of calculation and composite approach using the calorimetrically determined spin transition energies allow the estimation of the spin transition temperatures between fully LS and fully HS chains. The results are given in Table 8.

Fig. 7 Temperature dependencies of the Gibbs free energy of HHHHHH spin isomer of 1 for different values of ΔE_ad.

Fig. 8 Temperature dependencies of the free enthalpy of HHHHH (relative LLLLL one) spin isomer of 1 for different values of ΔE_ad.

Table 8 Transition temperatures estimated on the basis of the ΔG(T) diagrams from Fig. 2–6

Model molecule	1	2	2a	3
a Composite approach using the calorimetry data and vibrational data from DFT.
Tc experimental (K)	390↑	343↑	332↑	347↑
	345↓	335↓	311↓	311↓
Method
B3LYP*	Over 400 K	Over 400 K	Over 400 K	Over 400 K
B3LYP/geoB3LYP*	190 K	190 K	305 K	10 K
CAM-B3LYP/geoB3LYP*	160 K	150 K	320 K	120 K
B3LYP/D3	160 K	260 K	310 K	120 K
Cal/vib B3LYP*^a	315 K	100 K	65 K	350 K
Cal/vib B3LYP/D3^a	400 K	75 K	85 K	315 K

---

Interestingly, in spite of the obvious simplifications of the models (isolated molecules of finite size with no energetic and vibrational effects of crystal packing, within the harmonic approximation) give for systems 1 and 3 the reasonable estimation of the transition temperatures when the calorimetrically obtained values of spin transition energies are used.

Thermal analyses

[Fe(Htrz)₂trz]BF₄ (1) was studied by DSC on both warming and cooling modes over the temperature range 300–425 K at a 10 K min⁻¹. An endothermic peak is observed on warming at T_max^↑ = 387 K and an exothermic peak is recorded on cooling at T_max^↓ = 343 K (Fig. 9), in good agreement with the transition temperatures provided by Kröber et al.⁵⁹

Fig. 9 DSC measurements for [Fe(Htrz)₂trz]BF₄ (1) over the 300–423 K temperature range, at a scan rate of 10 K min⁻¹, in cooling and warming modes, for three consecutive runs.

The energy and entropy associated with the spin transition were evaluated by considering the fraction of switching sites evaluated from Mössbauer spectroscopy (86%). This correction was not done in the original DSC report on this compound, which explains why corrected values are lower (compared to ΔH_HL = 27 kJ mol⁻¹).⁵³

TGA was undertaken under N_2(g) (100 mL min⁻¹) from room temperature to 873 K with a heating rate of 5 K min⁻¹ to determine water content (Fig. 10). The ca. 10.8% mass loss below 390 K was attributed to the loss of 2.5 water molecules, corresponding to the formula [Fe(NH₂trz)₃]Cl₂·2.5H₂O. Thermogravimetric measurements were also performed in air (100 mL min⁻¹, 5 K min⁻¹) over the temperature range of 298 K to 433 K. The mass loss before 420 K also corresponds to the loss of 2.5 water molecules (Fig. 10) no mass change was observed with decreasing temperature, indicating that the sample did not rehydrate in air, thus affording compound 2.

Fig. 10 TGA curves of [Fe(NH₂trz)₃]Cl₂·2.5H₂O (2a) under N_2(g) (top) and under air (bottom).

DSC measurements were performed under a N_2(g) atmosphere with a heating rate of 5 K min⁻¹ in the temperature range from 223 K to 433 K. Calorimetric data of 2a were first analysed within the temperature range of 223 K to 433 K, revealing two distinct endothermic peaks. A sharp peak at T_max^↑ = 332 K corresponding to the SCO of 2a from the LS state to the HS state (see colour change in Fig. 11), while a broader peak at 370 K is associated with the dehydration process (Fig. 11).

Fig. 11 (Top) DSC curves of [Fe(NH₂trz)₃]Cl₂·2.5H₂O from 223 K to 433 K showing two peaks associated to SCO and to the removal of water molecules (red line) and [Fe(NH₂trz)₃]Cl₂ during cooling (blue line) and heating modes (green line). (Middle) DSC curves of [Fe(NH₂trz)₃]Cl₂·2.5H₂O during cooling (blue line) and heating (red line) modes. (Bottom) Color change of [Fe(NH₂trz)₃]Cl₂·2.5H₂O from pink (303 K) to white (332 K).

The warming mode was proceeded until complete dehydration observed at 433 K, according to thermogravimetry (Fig. 10). Subsequently, the dehydrated sample 2 was cooled to 223 K, during which an exothermic peak at T_max^↓ = 335 K was observed, indicating the transition from the HS state back to the LS state. During reheating, the transition from LS to HS was detected at 343 K, resulting in a hysteresis loop of 8 K width for 2. The energy and entropy changes were determined as ΔH = 8.93 kJ mol⁻¹ and ΔS = 26.06 J K⁻¹ mol⁻¹, respectively. DSC characterization of 2a was carried out in the range of 273 K to 343 K, below the dehydration temperature (Fig. 10) to allow the missing thermodynamic parameters to be determined. Correspondingly, endothermic and exothermic peaks corresponding to the SCO were detected at 332 K and 311 K, respectively. This resulted in the formation of a hysteresis loop of 21 K width, which is much broader than that observed for 2. Such an increase in hysteresis width may be due to an increase of intermolecular interactions involving non coordinated water molecules with the NH₂trz ligand and/or chlorine atoms in this system, as observed in the crystal structure of [Cu(NH₂trz)₃](NO₃)₂·H₂O.⁷¹ The changes in energy and entropy were estimated as ΔH = 12.22 kJ mol⁻¹ and ΔS = 39.38 J mol⁻¹ K⁻¹, respectively.

Experimental and computational methods

Syntheses. [Fe(Htrz)trz]BF₄ (1) was prepared at room temperature as a pink powder using an ethanol/water mixture of Fe(BF₄)₂·6H₂O mixed with three equivalents of Htrz.⁵³

[Fe(NH₂trz)₃]Cl₂·2.5H₂O (2a) was synthesized by adapting the reported synthesis.⁵⁴ Dissolve 2 mmol (400 mg) of FeCl₂·4H₂O along with a small amount of ascorbic acid in water (10 mL). In a separate container, dissolve 20 mmol (1.68 g) of NH₂trz in ethanol (15 mL). Stir the combined solutions at room temperature for 15 h. A white precipitate formed that gradually turned purple. The resulting precipitate was washed with ethanol and dried in air.

[Fe(NH₂trz)₃]Cl₂ (2) was prepared in situ; warming 2a above at a temperature where full dehydration occurs (420 K), accurately determined by thermogravimetric analysis (TGA), see below.

Thermal analyses. DSC experiments were carried out on 1 in a He_(g) atmosphere using a PerkinElmer DSC Pyris 1 instrument equipped with a cryostat and operating down to 98 K. Aluminum capsules were loaded with 21.9 mg of sample and sealed. The heating and cooling rates were fixed at 10 K min⁻¹. Temperatures and enthalpies were calibrated over the temperature range of interest (300–425 K) using the solid–liquid transitions of pure indium (99.99%). Enthalpy data were obtained by integration of the peaks using the PYRISTM DSC software 7.0 in specific heat C_p (J mol⁻¹ K⁻¹) format. The transition temperatures were derived by considering the maximum (T_max) of the thermal anomalies.

Thermogravimetric analyses were undertaken on a Mettler Toledo TGA/SDTA 851e analyzer under N_2(g) (100 mL min⁻¹) from room temperature to 873 K with a heating rate of 5 K min⁻¹.

Differential scanning calorimetry (DSC) measurements on 2 and 2a were performed on a Mettler Toledo DSC 3 using a N_2(g) atmosphere with a heating rate of 5 K min⁻¹ in the temperature range from 223 K to 433 K. Aluminium standard 40 mL capsules were loaded with 5.70 mg of sample and sealed. Temperatures and enthalpies were calibrated using pure indium and zinc. The specific heat capacity (C_p) was determined using STARe software.

The models. Most of the modelling was performed with the heptanuclear models represented in Scheme 1. For 1 and 2 they were based on the heptameric models reported for some of the spin isomers in ref. 47b and 47c, respectively. The model 2a was built up after all spin isomers of 2 were optimised and then 12 molecules of water were introduced to form the CH⋯O(H₂)···HC bridges between each two parallel arranged triazole ligands (see pdb files in SI). The chlorine anions were shifted on the rim of the chain in starting structures providing the stoichiometry of Fe₅(H₂O)₁₂ for the inner switching iron centres in the heptanuclear units. This stoichiometry is between Fe(H₂O)₂ reported in literature⁶² and the Fe(H₂O)_2.5 derived from thermogravimetric data in this paper (vide supra). The model for 3 was based on the crystal structure of the LS isomer.⁵⁶ While the calorimetric data for this complex are available for the anhydrous complex, the calculation were performed without water molecules. All system were modelled as the neutral molecules, with 2, 2a and 3 containing 14 anions, while the 1 one with 8 anions (there are six negatively charged triazole ligands in the structure of 1). Th structure of the HHHHH spin isomer of 3 was obtained on the basis of the optimised Zn₇ analogue starting from the LLLLL structure with the further optimisation of the HHHHH model. As described previously⁴⁷ all spin isomers containing both HS and LS centres in the inner pentanuclear unit were modelled on the basis of the optimised systems in which the HS Fe(II) were replaced with Zn(II). The Zn centres were subsequently replaced with Fe(II), the stability of the wave function was then checked with the Gaussian stable = opt option and the molecule was further optimised for all-Fe centres. Such an approach yielded all spin isomers under consideration. The H_coop parameter was additionally estimated using the nonanuclear and undecanuclear models of four isomers, two with fully HS or LS inner Fe centres and two with one Fe centres of the other spin (i.e. the analogues of LLLLL, HHHHH, LLHLL and HHLHH, respectively). The model molecules were built up on the basis of the optimised heptanuclear models. In order to decrease the computational time the terminal HS centres with three coordinated waters were replaced with the Zn ones.

DFT methods. As stated previously, all optimisations and normal mode calculations were firstly performed with B3LYP* functional⁶³ using Gaussian 16 package,⁷² Then the calculations were performed with B3LYP⁶⁶ functional using the dispersion correction. For 1 and 2 the further optimization were performed with CAM-B3LYP⁶⁷ functionals and the CEP-31G basis.⁷³ Additionally, the PBEh functional was used for the above two systems.⁷⁴ In order to check the dependence of the obtained H_coop parameters the relevant four isomers of 1, 2 and 2a were optimised with B3LYP, CAM-B3LYP and TPSS⁷⁵ functionals with Grimme's D3 dispersion correction.⁷⁶ The CEP-31G basis set was used for the former and the tzvp⁷⁷ for the latter one. Also, the TPSSh functional^76,78 was used together with the tzvp basis set. We used the default settings of Gaussian 16, optimising the internal coordinates with Berny algorithm⁷³ (see https://gaussian.com/opt/ for details). The further details concerning the threshold values for optimisation, integral grid, the dispersion functions, etc. are given in SI.

Conclusions

Four different 1D chain spin crossover Fe(II) complexes with 1,2,4-triazole based ligands exhibiting sharp or hysteretic spin transitions have been modelled by DFT methods using heptanuclear model fragments. In each case 20 possible spin isomers were calculated varying the permutation of spin centres within the inner pentanuclear fragment. The B3LYP*/CEP-31G calculations revealed that each spin isomer has a distinct electronic and vibrational energy. Further energy calculations with the B3LYP and CAM-B3LYP functionals for the geometry of the spin isomers optimised with B3LYP* exhibit that each spin transition energy from the fully LS state to a given spin isomer may be considered as the sum of the functional-dependent quantity E_ad (derived from the calculated spin transition energy from the fully LS to the fully HS state) and the functional independent H_strain. The latter is the result of the strain in the particular spin isomer due to short-range interactions between HS and LS centres. Calculations of the temperature dependence of the Gibbs free energy change on the spin transition from the LLLLL isomer to all 19 spin isomers were performed. These results point towards ferroelastic properties of the isolated 1D chains under study suggesting a direct switch between the fully LS and HS isomers at a sufficiently low value of E_ad. It was also found that short-range interactions within the chain affect even centres that are 1.6 nm apart (i.e. with 3 Fe(II) centres in between) in an HS matrix with LS defects. The corresponding interactions in the LS matrix occur at significantly lower distances (1.1 nm for 3). The comparison between obtained results and the experimental finding may be summarized in three points: (i) all four systems under study reveal the sharp transitions indicating the ferroelastic character of both short- and long-range interactions. Our results show that the ferroelastic properties are present in the relatively short fragments of chains of these molecules. Additionally, the most destabilised spin isomers seem to be those with antiferroelastic structure. (ii) Also for complexes 2a and 3 the effect of elongation of the bonds of the LS centres by the HS neighbours and the shortening of the HS centres by the LS neighbours is predicted with DFT methods. This effect was previously predicted and experimentally observed with nuclear inelastic scattering for 1 and 2.

(iii) The composite approach combining the DFT calculated H_strain parameters and vibrational properties with the calorimetrically determined spin transition energies gives the reasonable estimate of the experimentally determined spin transition energies for two of the four studied compounds.

We believe that the above results may be helpful in the further development of the phenomenological models.

Author contributions

JAW, HP and VS conceived the paper, MD, XL and YG performed and evaluated the calorimetric measurements, JAW performed the DFT calculations, KG, TH and JAW calculated the thermodynamic functions, XL synthesized the complexes, JAW, YG, HP and VS wrote the paper.

Conflicts of interest

There are no conflicts to declare.

Data availability

The calculated relative energies are given in text. In SI the calculated optimised structures are given as xyz file and harmonic frequencies for each modelled systems are included.

Also the computational results and details, and details of thermochemical measurements are included. See DOI: https://doi.org/10.1039/d5ra03472h.

Acknowledgements

V. S. and J. A. W. acknowledge the support by the Deutsche Forschungsgemeinschaft (DFG) through CRC/TRR173, “Spin + X”, project A4. V. S. also acknowledges support by the German Ministry of Research (BMBF) under 05K22UK1. J. A. W. and V. S. are grateful to Allianz für Hochleistungsrechnen Rheinland–Pfalz (AHRP) for providing CPU-time within the project RPTU-SPINPLUSXA4. Y. G. thank the Fonds de la Recherche Scientifique-FNRS (PDR T.0095.21, CDR J.0064.23, EQP U.N027.24). X. Li was supported by a fellowship from the China Scholarship Council (202108310050). We would like to thank the reviewers or for their valuable comments, which helped improve the manuscript.

Notes and references

1. (a) P. Gütlich and H. A. Goodwin, Top. Curr. Chem., 2003, 233–235; (b) M. H. Halcrow, Spin-Crossover Materials, John Wiley & Sons, 2013.
2. G. Molnár, M. Mikolasek, K. Ridier, A. Fahs, W. Nicolazzi and A. Bousseksou, Ann. Phys., 2019, 531, 1900076.
3. J. Pavlik and R. Boča, Eur. J. Inorg. Chem., 2013, 697 and references therein..
4. F. Varret, K. Boukheddaden, E. Codjovi, C. Enanescu and J. Linares, Top. Curr. Chem., 2004, 234, 199 and references therein..
5. A. Palii, S. Ostrovsky, O. Reu, B. Tsukerblatt, S. Descurtins, S.-X. Liu and S. Klokishner, J. Chem. Phys., 2015, 143, 084502 and references therein..
6. A.-I. Popa, L. Stolieriu and C. Enachescu, J. Appl. Phys., 2021, 129, 131101 and references therein..
7. D. C. Ashley and E. Jakubikova, Coord. Chem. Rev., 2017, 337, 97 and references therein..
8. (a) S.-G. Coca and E. Ruiz, Inorg. Chem., 2024, 63, 13338; (b) S. Rajpurohit, V. Vennelakanti and H. J. Kulik, J. Chem. Phys. A., 2024, 128(41), 9082.
9. (a) F. Daumann, G. Hörner, P. Köhler and B. Weber, Cryst. Growth Des., 2025, 25, 4565; (b) Y.-Y. Wu, Z.-Y. li, S. Peng, Z.-Y. Zhang, H.-M. Cheng, H. Su, W.-Q. Houm, F.-L. Yang, S.-Q. Wu, O. Sato, J.-W. Dai, W. Li and X.-H. Bu, J. Am. Chem. Soc., 2024, 146, 8206.
10. P. R. Spackman, A. Grosjean, S. P. Thomas, D. P. Karothu, P. Naumov and M. A. Spackman, Angew. Chem., Int. Ed., 2022, 134, e202110716 and references therein..
11. J. Maul, D. Ongari, S. M. Moosavi, B. Schmidt and A. Erba, J. Phys. Chem. Lett., 2020, 11, 8543.
12. (a) L. Wiehl, G. Kiel, C.-P. Köhler, H. Spiering and P. Gütlich, Inorg. Chem., 1986, 25, 1565; (b) J. Jung, J. Schmitt, A. Hauser, L. Wiehl, K. Knorr, H. Spiering and P. Gütlich, Z. Phys. B, 1996, 100, 523.
13. H. Spiering, Top. Curr. Chem., 2004, 235, 171.
14. M. Radon, Phys. Chem. Chem. Phys., 2023, 25, 30800.
15. (a) S. R. Mortensen and K. P. Kepp, J. Phys. Chem. A, 2015, 119, 4041; (b) D. Vidal, J. Cirera and J. Ribas-Arino, Dalton Trans., 2021, 50, 17635; (c) H. Paulsen, V. Schünemann and J. A. Wolny, Eur. J. Inorg. Chem., 2013, 628; (d) H. Paulsen and A. X. Trautwein, Top. Curr. Chem., 2004, 235, 197; (e) S. Rajpurohit, V. Vennelekanti and H. J. Kulik, J. Phys. Chem. A, 2024, 128, 9082; (f) G. Paveliuc and L. M. Lawson Daku, J. Phys. Chem. A, 2024, 128, 8404; (g) V. Vennelakanti, M. G. Taylor, A. Nandy, C. Duan and H. J. Kulik, J. Chem. Phys., 2023, 159, 024120; (h) J. Cirera, M. Via-Nadal and E. Ruiz, Inorg. Chem., 2018, 57, 14097.
16. G. Lemercier, N. Bréfuel, S. Shova, J. A. Wolny, F. Dahan, M. Verelst, H. Paulsen, A. X. Trautwein and J.-P. Tuchagues, Chem.–Eur. J., 2006, 12, 742.
17. S. Vela and H. Paulsen, Dalton Trans., 2019, 48, 1237.
18. S. Vela and H. Paulsen, Inorg. Chem., 2018, 57, 9478.
19. (a) H. Paulsen, Magnetochemistry, 2016, 2, 14; (b) L. Kreuzburg, C. G. Hübner and H. Paulsen, Materials, 2017, 10, 172.
20. S. G. Wu, L.-F. Wang, Z.-Y. Ruan, S.-N. Du, S. Gómez-Coca, Z.-P. Ni, E. Ruiz, X.-M. Chen and M.-L. Tong, J. Am. Chem. Soc., 2022, 144, 14888.
21. Y. Zhang, J. Chem. Phys., 2014, 141, 214703 and references therein..
22. S. Vela, M. Fumanal, J. Rinas-Arino and V. Robert, Phys. Chem. Chem. Phys., 2015, 17, 16306.
23. M. Ohlrich and B. J. Powell, J. Chem. Phys., 2020, 153, 104107.
24. D. S. Middlemis, D. Portinari, C. P. Grey, C. A. Morrison and C. C. Wilson, Phys. Rev., 2010, B81, 184410.
25. J. Wojdeł, I. d. P. R. Moreira and F. Illas, J. Chem. Phys., 2009, 130, 014702.
26. T. Bučko, J. Hafner, S. Lebègue and J. Ángyán, Phys. Chem. Chem. Phys., 2012, 14, 5389.
27. Y. Zhang, J. Chem. Phys., 2014, 141, 214703.
28. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni and I. Dabo, J. Phys.:Condens. Matter, 2009, 21, 395502.
29. E. Artacho, D. Sánchez-Portal, P. Ordejón, A. García and J. M. Soler, Phys. Status Solidi B, 1999, 215, 809.
30. K.-Y. Liu, J. Liu and J. M. Herbert, J. Comp. Chem., 2017, 38, 1678.
31. E. Collet, G. Azzolina, T. Ichii, L. Guerin, R. Bertoni, A. Moréac, M. Cammarata, N. Daro, G. Chastanet, J. Kubicki, K. Tanaka and S. F. Matar, Eur. Phys. J. B, 2019, 92, 12.
32. G. Kresse and J. Furthmüller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169.
33. Á. Fernández-Blanco, L. Piñeiro-López, M. Jiménez-Ruiz, S. Rols, J.-A. Real, J. A. Rodríguez-Velamazán and R. Poloni, J. Chem. Phys. C., 2022, 126, 8090.
34. S. Baroni, S. de Gironcoli, A. Dal Corso and P. Giannozzi, Rev. Mod. Phys., 2001, 73, 515.
35. Y.-Y. Wu, Z.-Y. Li, S. Peng, Z. Y. Zhang, H.-M. Cheng, H. Su, W.-Q. Houm, F.-L. Yang, S.-Q. Wu, O. Sato, J.-W. Dai, W. Lim and X.- H Bu, J. Am. Chem. Soc., 2024, 146(12), 8206.
36. (a) A. P. Thompson, H. Metin Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, Comp. Phys. Comm., 2022, 271, 108171; (b) J. C. Phillips, D. J. Hardy, J. D. C. Maia, J. E. Stone, J. V. Ribeiro, R. C. Bernardi, R. Buch, G. Fiorin, J. Henin, W. Jiang, R. McGreevy, M. C. R. Melo, B. K. Radak, R. D. Skeel, A. Singharoy, Y. Wang, B. Roux, A. Aksimentiev, Z. Luthey-Schulten, L. V. Kale, K. Schulten, C. Chipot and E. Tajkhorshid, J. Chem. Phys., 2020, 153, 044130.
37. (a) R. Meyer, C. Mücksch, J. A. Wolny, V. Schünemann and H. M. Urbassek, Chem. Phys. Lett., 2019, 733, 136666; (b) J. Herz, R. Meyer, J. A. Wolny, V. Schünemann and H. M. Urbassek, Appl. Phys. A: Mater. Sci. Process., 2023, 129, 345.
38. S. Mi, A. Fahs, G. Molnar, W. Nicolazzi and A. Bousseksou, Chem. Phys. Lett., 2023, 811, 140232.
39. C. Huy Pham and F. Paesani, J. Phys. Chem. Lett., 2016, 7, 4022.
40. C. Huy Pham, J. Cireara and F. Paesani, J. Am. Chem. Soc., 2016, 138, 6123.
41. (a) R. J. Deeth, Coordination Chemistry, Springer, Berlin Heidelberg, 1995, vol. 82, p. 1; (b) R. J. Deeth, Coord. Chem. Rev., 2001, 212, 4; (c) R. J. Deeth, A. Anastasi, C. Diedrich and K. Randell, Coord. Chem. Rev., 2009, 253, 795.
42. (a) A. Erba, J. Maul, M. Ferrabone, P. Carbonnière, M. Rérat and R. Dovesi, J. Chem. Theory Comput., 2019, 15(6), 3755 and references therein;; (b) G. Comas-Vilà and P. Salvador, Inorg. Chem., 2025, 64, 15347 and references therein..
43. I. M. Alecu, J. Zheng, Y. Zhao, D. G. Truhlar, J. Chem. Theory Comput. 2010, 6, 9, 3755 and references therein. See https://comp.chem.umn.edu/freqscale/index.html for the database of these factors.
44. P. Pracht and S. Grimme, Chem. Sci., 2021, 12, 6551.
45. J. Wu, C. Sousa and C. de Graaf, Magnetochemistry, 2019, 5(3), 49.
46. O. Markin, I. Gudyma and K. Boukheddaden, Phys. Rev. B, 2024, 110, 174416.
47. (a) S. Rackwitz, J. A. Wolny, K. Muffler, K. Achterhold, R. Rüffer, Y. Garcia, R. Diller and V. Schünemann, Phys. Chem. Chem. Phys., 2012, 14, 14650; (b) J. A. Wolny, I. Faus, J. Marx, R. Rüffer, A. I. Chumakov, K. Schlage, H.-C. Wille and V. Schünemann, Magnetochemistry, 2016, 2, 19; (c) K. Jenni, L. Scherthan, I. Faus, J. Marx, C. Strohm, H.-C. Wille, P. Würtz, V. Schünemann and J. A. Wolny, Phys. Chem. Chem. Phys., 2017, 19, 18880; (d) J. A. Wolny, T. Hochdörffer, S. Sadashivaiah, H. Auerbach, K. Jenni, L. Scherthan, A.-M. Li, C. von Malotki, H.-C. Wille, E. Rentschler and V. Schünemann, J. Phys.: Condens. Matter, 2021, 33, 034004.
48. H. Paulsen, V. Schünemann, A. X. Trautwein and H. Winkler, Coord. Chem. Rev., 2005, 247, 25.
49. T. Hochdörffer, A. I. Chumakov, H.-C. Wille, V. Schünemann and J. A. Wolny, Dalton Trans., 2019, 48, 15625.
50. S. Rackwitz, W. Klopper, V. Schünemann and J. A. Wolny, Phys. Chem. Chem. Phys., 2013, 15, 15450.
51. A. Moneo-Corcuera, D. Nieto-Castro, J. Cirera, V. Gómez, J. Sanjosé -Orduna, C. Casadevall, G. Molnar, A. Bousseksou, T. Parella, J.-M. Martánez-Agudo, J. Lloret-Fillol, M. H. Perez-Temprano, E. Ruiz and J.-R. Galán-Mascarós, Chem, 2023, 9, 1.
52. K. Boukheddaden, S. Miyashita and M. Nishino, Phys. Rev. B:Condens. Matter Mater. Phys., 2007, 75, 094112.
53. M. Nishino and S. Miyashita, Phys. Rev. B:Condens. Matter Mater. Phys., 2013, 88, 014108.
54. M. C. Muñoz and J. A. Real, in Spin-Crossover Materials, ed. M. A. Halcrow, John Wiley & Sons, 2013.
55. S. A. Siddiqui, O. Domanov, E. Schafler, J. Vejpravova and H. Shiozawa, J. Mat. Chem. C, 2021, 9, 1079.
56. D. Nieto-Castro, F. A. Garcés-Pineda, A. Moreno-Corcuera, B. Pato-Doldan, F. Gispert-Guirado, J. Benet-Buchholz and J. R. Galan-Mascarós, Inorg. Chem., 2020, 59, 7953.
57. M. Pedrahita-Bello, K. Ridier, M. Mikolasek, G. Molnár, W. Nicolazzi, L. Salmon and A. Bousseksou, Chem. Commun., 2019, 55, 4769.
58. J. Linarès, H. Spiering and F. Varret, Eur. J. Phys. B, 1999, 10, 271.
59. J. Kröber, J. P. Audière, R. Claude, E. Codjovi, O. Khan, J. G. Haasnoot, F. Grolière, F. Jay, A. Bousseksou, J. Linarès, F. Varret and A. Gonthier-Vassal, Chem. Mater., 1994, 6, 1404.
60. L. G. Lavrenova, O. G. Shakirova, V. N. Ikorskii, V. A. Varnek, L. A. Sheludyakova and S. V. Larionov, Russ. J. Coord. Chem., 2003, 29, 22–27.
61. L. G. Lavrenova, V. N. Ikorskii, V. A. Varnek, I. M. Oglezneva and S. V. Larionov, Koord. Khim., 1986, 12, 207.
62. A. Grosjean, N. Daro, B. Kauffmann, A. Kaiba, J.- F- Létard and Ph. Guionneau, Chem. Commun., 2011, 47, 12382.
63. (a) W. J. Stevens, H. Basch and J. Krauss, J. Phys. Chem., 1984, 81, 6026; (b) W. J. Stevens, M. Krauss, H. Basch and P. G. Jasien, Can. J. Chem., 1992, 70, 612; (c) T. R. Cundari and W. J. Stevens, J. Chem. Phys., 1993, 98, 5555.
64. (a) M. M. Dîrtu, A. Rotaru, D. Gillard, J. Linares, E. Codjovi, B. Tinant and Y. Garcia, Inorg. Chem., 2009, 48, 7838; (b) O. Roubeau, M. Castro, R. Burriel, J. G. Haasnoot and J. Reedijk, J. Phys. Chem. B, 2011, 115, 3003; (c) G. Berezovskii, M. Bushuev, D. Pishchur and L. Lavrenova, J. Therm. Anal. Calorim., 2008, 93, 999.
65. T. Yanai, D. Tew and N. Handy, Chem. Phys. Lett., 2004, 393, 51.
66. A. D. Becke, J. Chem. Phys., 1993, 98, 5648.
67. R. Bertoni, E. Collet, H. Cailleau, M.-L. Boillot, A. Tissot, J. Laisney, C. Enachescu and M. Lorenc, Phys. Chem. Chem. Phys., 2019, 21, 6606.
68. K. L. Ronayne, H. Paulsen, A. Höfer, A. C. Dennis, J. A. Wolny, A. I. Chumakov, V. Schüneman, H. Winkler, H. Spiering, A. Bousseksou, P. Gütlich, A. X. Trautwein and J. J. McGarvey, Phys. Chem. Chem. Phys., 2006, 4685.
69. M. M. Dîrtu, C. Neuhausen, A. D. Naik, A. Rotaru, L. Spinu and Y. Garcia, Inorg. Chem., 2010, 49, 5723.
70. H. Butenschön, K. P. C. Volhardt and N. E. Schore, Organische Chemie, Wiley-VCH, 2001, p. 90.
71. (a) M. Reiher, O. Salomon and B. A. Hess, Theor. Chim. Acta, 2001, 107, 49; (b) O. Salomon, M. Reiher and B. A. Hess, J. Chem. Phys., 2002, 117, 4729.
72. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Calmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, B. Foresman and D. J. Fox, Gaussian 16, Revision A.03, Gaussian, Inc., Wallingford CT, 2016.
73. (a) W. J. Stevens, H. Basch and M. Krauss, J. Phys. Chem., 1984, 81, 6026; (b) W. J. Stevens, M. Krauss, H. Basch H and P. G. Jasien, Can. J. Chem., 1992, 70, 612; (c) T. R. Cundari and W. J. Stevens, J. Chem. Phys., 1993, 98, 5555.
74. M. Ernzerhof and J. P. Perdew, J. Chem. Phys., 1998, 109, 3313.
75. J. M. Tao, J. P. Perdew, V. N. Staroverov and G. E. Scuseria, Phys. Rev. Lett., 2003, 91, 146401.
76. S. Grimme, J. Antony, S. Ehrlich and H. Krieg, J. Chem. Phys., 2010, 132, 154104.
77. (a) A. Schaefer, H. Horn and R. Ahlrichs, J. Chem. Phys., 1992, 97, 2571; (b) A. Schaefer, C. Huber and R. Ahlrichs, J. Chem. Phys., 1994, 100, 5829.
78. V. N. Staroverov, G. E. Scuseria, J. Tao and J. P. Perdew, J. Chem. Phys., 2003, 119, 12129.

---

Footnote

† Dedicated to the memory of Prof. Dr hab. Cristian Enachescu (1976–2024) who devoted his career to the computation of spin crossover systems.

---