Pitfalls on evaluating pair exchange interactions for modelling molecule-based magnetism

Abstract

Molecule-based magnetism is a solid-state property that results from the microscopic interaction between magnetic centres or radicals. The observed magnetic response is due to unpaired electrons whose coupling leads to a particular magnetic topology. Therefore, to understand the magnetic response of a given molecule-based magnet and reproduce the available experimental magnetic properties by means of statistical mechanics, one has to be able to determine the value of the J_AB magnetic exchange coupling between radicals. The calculation of J_AB is thus a key point for modelling molecule-based magnetism. In this Perspectives article, we will build upon our experience in modelling molecular magnetism to point out some pitfalls on evaluating J_AB couplings. Special attention must be paid to the cluster models used to evaluate J_AB, which should account for cooperative effects among J_AB interactions and also consider the environment (counterions, hydrogen bonding) of the two radicals whose interaction has to be evaluated. It will be also necessary to assess whether a DFT-based or a wavefunction-based method is best to study a given radical. Finally, in addition to model and method, the J_AB couplings have to be able to adapt to changes in the magnetic topology due to thermal fluctuations. Therefore, it is most important to appraise in which systems molecular dynamics simulations would be required. Given the large number of issues one must tackle when choosing the correct model and method to evaluate J_AB interactions for modelling magnetic properties in molecule-based materials, the “human factor” is a must to cross-examine and challenge computations before trusting any result.

---

Maria Fumanal

Dr Maria Fumanal obtained her PhD in 2015 at the University of Barcelona under the supervision of Prof. Juan Novoa and Dr Jordi Ribas-Arino. Her PhD focused on the modeling of organic π-stacked architectures for the development of magnetic switches. She joined after the University of Strasbourg to work with Dr Chantal Daniel in excited states dynamics of functional Re(I) complexes. Since 2018, she works at the École Polytechnique Fédérale de Lausanne where she was awarded a Marie Skłodowska-Curie Individual Fellowship to perform her research in the design of copolymers for singlet fission in the group of Prof. Clémence Corminboeuf.

Joaquim Jornet-Somoza

Dr Joaquim Jornet-Somoza currently works as Scientific Computing Expert at the Advanced Research Facilities (UPV/EHU, Spain) and is a guest researcher at The Max Planck Institute for the Structure and Dynamics of Matter, Hamburg (Germany). Dr Jornet-Somoza did his PhD on the Theoretical Study of Molecular Magnets under the supervision of Dr Deumal and Prof. Novoa. He was awarded with several national and international postdoctoral grants, like the Marie Curie Slwadoska Individual Fellowship in 2019, and became an expert on the ab initio study of excited states properties of condensed matter and its interaction with electromagnetic fields.

Sergi Vela

Dr Sergi Vela graduated with a chemistry degree from the University of Barcelona in 2009. He obtained his PhD in 2014 under the supervision of Prof. Juan J. Novoa and Prof. Mercè Deumal. Then, he joined the Laboratoire de Chimie Quantique at the University of Strasbourg, where he worked in the computational modelling of Spin Crossover materials with Prof. Vincent Robert. In 2018, he joined the group of Prof. Clémence Corminboeuf at EPFL under a Marie Sklodowska-Curie Fellowship, to work in the study of molecular photo-switches, and on the design and discovery of novel materials within the Marvel program.

Juan J. Novoa

Dr Juan J. Novoa obtained his PhD in 1981 at the University of Barcelona (UB). He then did postdoctoral stays with Prof. R. Carbó and Prof. M. A. Robb. Since 1997, he has been Full Professor at the Department of Materials Science and Physical Chemistry of UB. During his career he has been a visiting scientist in several advanced research laboratories (CRAY Research, Minnesota; IBM Research Laboratory in Zürich), and Visiting Professor at several universities (North Carolina State University, Univ. Bologna, Clark University, Univ. of Utah, Univ. of Sao Carlos, and Osaka Prefecture University). His current research focuses on computational studies of molecular materials.

Jordi Ribas-Arino

Dr Jordi Ribas-Arino received his PhD in 2006 with Prof. Juan J. Novoa at the University of Barcelona (UB). During his PhD, he performed two research stays at the U.S. Naval Research Laboratory, with Dr Mark R. Pederson. Thereafter, he worked as a Postdoctoral Researcher at the Ruhr-Universität-Bochum with Prof. Dominik Marx. A Humboldt Fellowship supported part of his postdoc. In 2010, he returned to UB with a “Ramón y Cajal” contract and since 2015, he has been Associate Professor at the same university. His current research interests concern the development and application of computational tools to study multifunctional molecule-based materials.

Mercè Deumal

Dr Mercè Deumal is full Professor at Secció de Química Física (Universitat de Barcelona), where she develops her tasks both as a researcher and lecturer at BSc, MSc and PhD levels. She is the local coordinator at UB of the interuniversity doctoral program on ‘Theoretical Chemistry and Computational Modeling’, and belongs to ‘Institut de Química Teòrica i Computacional’ (awarded as Unit of Excellence María de Maeztu), and ‘Xarxa de Referència en Química Teòrica i Computacional’. Her research is framed in the field of Material Science and Computational Chemistry, and focuses on the rational design of advanced multifunctional molecular materials using multiscale simulation methods, which encompass the use of ab initio quantum mechanics, molecular dynamics, Monte Carlo, statistical mechanics, machine learning and solid state calculations.

---

1. Introduction

In molecule-based magnetic materials, the microscopic interactions between magnetic centres or radicals, whose unpaired electrons couple leading to a specific magnetic topology, play a key role in defining their macroscopic magnetic properties.¹ There are two main strategies to reproduce experimental magnetic properties (e.g. magnetic susceptibility, heat capacity, magnetization) by means of statistical mechanics² depending on whether the electronic structure is explicitly considered or not. The latter approach aims at providing a numerical magnetic response for a given molecular material. A very successful approach,^3a which does not explicitly account for the crystal geometry (except in a parametric form³), is based on the use of a generalized spin Hamiltonian. Alternatively, the electronic structure-based strategy targets on capturing the microscopic complexity of the molecular material to understand, as well as reproduce, the experimental magnetic data.⁴ As computational chemists, we seek modelling molecule-based materials as a tool to rationalize their physical properties. Accordingly, within the framework of molecular magnetism, one has to be able to evaluate the magnitude of the J_AB magnetic exchange coupling between radicals. The calculation of the J_AB radical⋯radical interaction is thus one of the key points for modelling magnetism, and can be a demanding task if the number of J_AB couplings is large. At this point, anyone working in the field of molecular magnetism has wished there was a fully automated procedure to determine the magnetic response of a given molecule-based compound in a black-box fashion. Such automated code would most probably resort to the simplest model (dimer) and the least expensive method (DFT broken symmetry level)⁵ to evaluate all J_AB coupling interactions.^5b,c,6 Also the magnetic unit cell would perforce have to be found in an automated way. Model, method and magnetic unit cell would have a difficult validation since, by definition, there is no diagnosis to check whether the results obtained from a “black-box procedure” are correct or not. We must acknowledge that the selection of the right model and method to compute the J_AB magnetic coupling is sometimes very challenging. Therefore, due to the large number of issues one must tackle, the automatized selection of the model and method can be an extremely difficult task. The take-home message of this paper is not that its design is not feasible, but that there must be some careful thought behind it. This Perspective paper aims at reflecting on how molecule-based magnetism is modelled, which pitfalls could be encountered, and how to detect and prevent them.

Usually, the working strategy that we follow to calculate macroscopic magnetic properties from microscopic magnetic interactions among radicals uses only the knowledge of X-ray crystallographic data. A four-step procedure⁷ can be followed. First, the crystal packing is analysed to select pairs of A, B radicals that might be magnetically important, in terms of radical⋯radical distances. Secondly, once all those A⋯B dimers that are univocally defined within the crystal have been identified, the microscopic J_AB radical⋯radical magnetic interactions are computed using quantum chemistry methods. Depending on the system, one can resort to a variety of either Density Functional Theory DFT⁵ or wavefunction-based (e.g. CASSCF, RASSCF, DDCI)⁸ methods. Also depending on the system, one can resort to a variety of different schemes to obtain the numerical value of J_AB magnetic coupling: localized vs. delocalized, projected vs. unprojected, etc.^5,6 Due to the molecular nature of the magnetic materials we are interested to study, the J_AB magnetic coupling between radicals is of short-range, which allows evaluating the J_AB interactions using truncated cluster models. The simplest dimer model works in most of the cases. For a pair of S = 1/2 radicals, the J_AB microscopic exchange interaction is calculated as the energy difference between open-shell singlet and triplet states.⁵ However, a dimer model is not always the best option to compute J_AB. The enlargement of the size of the cluster model (e.g. using tetramer, decamer models) is sometimes a must to account for cooperative effects among J_AB interactions not included in a radical–radical pair calculation (i.e. dimer model).⁹ The selection of the correct model must also consider the environment of the two radicals whose interaction has to be evaluated, namely counterions,¹⁰ hydrogen bonding,¹¹etc. To sum up, irrespective of the model, within this approach, the J_AB interactions are evaluated at the fixed relative positions (an X-ray resolved structure or, alternatively, an optimized structure¹²) of the selected pairs of radicals within the crystal. Since the thermal oscillations of the spin carrying units around their equilibrium positions are ignored in this type of analysis, it can be stated that the standard approach is based on a static perspective. However, for materials with dominant exchange interactions propagating through π–π labile networks,¹³ one might have to resort to a dynamic perspective, in which thermal vibrations are explicitly considered.¹⁴ Therefore, in order to get a physically correct interpretation of the magnetic response of a given molecular material one has to be able to determine whether the static perspective will suffice or, on the contrary, molecular dynamics simulations would be required.

Let us now digress to comment on the possibility of using a periodic approach rather than a cluster model to evaluate the J_AB magnetic coupling between radicals. This alternative might appeal as being an easier option since unit cells with different spin settings could be used to calculate the J_AB magnetic coupling interactions between different spin centres without having to worry about the representation of the environment. Yet performing periodic computations on the crystal structure is comparatively more difficult for two main reasons. First, the number of radicals in the unit cell (Z) can be large, which would require the evaluation of a large number of spin states that, in most cases, will become intractable (e.g. a fairly small Z = 4 in ref. 9b generates over 14 potential J_AB magnetic interactions that need to be assessed). Moreover, in most cases those spin states are not the ground state for a given spin multiplicity and, thus, the definition and – especially – the preservation of the desired spin configuration represents a challenge for solid state codes. Second, evaluation of J_AB interactions with accurate methods in solid state has an extremely large computational cost and, in many cases, both non-hybrid and hybrid DFT functionals accessible for periodic computations may have difficulties in capturing the true electronic nature of the radicals within the crystal. For these reasons, cluster models offer great advantages to calculate J_AB magnetic interactions in molecular materials at high level accurate methods. Still, the choice of the actual cluster model (dimer, tetramer, etc.) needs to be validated to ensure the reliability of the J_AB values.

Once all J_AB's have been evaluated, the magnetic topology can be next defined in terms of all computed non-negligible J_AB magnetic coupling interactions. From the magnetic topology, we select a representative magnetic model, which will enable us to consider the unpaired electrons being completely coupled antiferromagnetically (AFM), completely coupled ferromagnetically (FM), or in any other of the many states obtained as a result of solving the secular equation problem together with the energy spectra and corresponding spin quantum numbers. It is important to stress that the definition of the magnetic model is crucial since the full diagonalisation¹⁵ will be done in the space expanded by this model. One could think that the selection of the magnetic model is a question of computer power and should not be difficult to automatically build a model Hamiltonian based on the hierarchy of the J_AB magnetic couplings. However, the choice of the most appropriate magnetic model is inherently difficult in molecule-based magnets as they are usually characterised by many J_AB magnetic couplings. As a consequence, convergence with respect to the model Hamiltonians size is required, which complicates its automatic selection.

Finally, the resulting energy spectra and spin quantum numbers are used in the appropriate statistical mechanics expressions to simulate the macroscopic magnetic properties.² Here we must highlight that all simulations are always compared to the available experimental data for validation purposes. Therefore, one can calculate magnetic susceptibility, heat capacity, or magnetization, depending on the available experimental data.

In summary, the working strategy above outlined is first-principles (FP) because the J_AB magnetic interactions are computed using DFT or other wavefunction-based techniques from crystallographic data, and is bottom-up (BU) because we use microscopic data to calculate macroscopic properties. This First-Principles Bottom-Up FPBU⁷ working strategy has been applied successfully to compute the magnetic response of a variety of prototypical molecule-based magnets of different magnetic dimensionality. However, it is not just the magnetic response that is targeted, since we have been able to understand and rationalize the magnetic behaviour depending on a range of different situations such as influence of temperature,¹⁶ effect of the Madelung field,¹⁰ effect of the crystal defects,¹⁷ effect of the ligands,^11,18 magnetic dimensionality,¹⁹ bistability²⁰…. At this point, one has to realize that a too simple automated procedure would imply applying too many unsurveyed restrictions: dimer model and unrestricted DFT to compute energies and, in turn, J_AB; estimation of the magnetic unit cell upon quantitative criteria; etc. In fact, the only attempt we have come across which automatically handle all these issues is applied to neutral organic-based 1D magnets.^4d The choice of these materials greatly simplifies calculations. Yet it would most probably face serious challenges when applied to more complex materials (as the examples we discuss below). Therefore, however appraising, black-box and modelling magnetism do not necessarily hold hands. Basically, a black-box procedure based on the static evaluation of J_AB's using a dimer model with unrestricted DFT limits predictability and might only be useful for description purposes. Description purposes are obviously important. However, as computational chemists, we usually have more questions than answers, and the questions usually come out unexpectedly. Therefore using a too simple automated code is thoughtless. In fact, we will next proceed to discuss issues we must consider in terms of models and in terms of methods in order to evaluate the radical⋯radical magnetic interactions and, in the long run, simulate magnetism.

First, we will show that the dimer model cannot be taken for granted to evaluate the magnetic J_AB interaction between two radicals. In fact, the size of the cluster model to calculate J_AB must be assessed in order to be certain about it. Note that a wrong value of J_AB coupling will lead to the wrong magnetic topology which, in turn will give rise to the wrong minimal magnetic model whose eigenvalues have to be used to simulate the relevant magnetic data. This is precisely the case of Cu₂(1,4-diazacycloheptane)₂Cl₄ (CuHpCl in Fig. 1a):^9b,21 the use of the wrong model to evaluate J_AB results in a spin-ladder magnetic topology instead of a 3D network of interconnected squared-plaquettes. Secondly, we will address the choice of the right method to calculate J_AB. The calculation of a given J_AB involves the evaluation of the energy of, at least, two states. This energy evaluation can be done by means of either DFT-based or wavefunction-based methods. Since DFT calculations are less demanding, in terms of resources and computational time, it is considered as the default method to be used to calculate the energy of the states involved. Although DFT is known to provide acceptable values for J_AB coupling in many systems,^{9–11,14,16–20} we will exemplify which are the effects of choosing the wrong method with the phenylsemiquinone-bridged bisdithiazolyl (PhBBO in Fig. 1b) compound.^22,23 Finally, we will address the use of molecular dynamics simulations when studying bistable compounds whose radicals pack forming π-stacks along a given crystallographic direction and, thus, thermal fluctuations might have an impact on the magnetic topology. We have encountered that J_AB couplings have to be able to adapt to changes due to intermolecular vibrations as in the case of 1,3,5-trithia-2,4,6-triazapentalenyl (TTTA in Fig. 1c),^14a,24 which undergoes phase transition between a low temperature LT and a high temperature HT phases.

Fig. 1 Chemical formula for (a) Cu₂(1,4-diazacycloheptane)₂Cl₄, CuHpCl, (b) phenylsemiquinone-bridged bisdithiazolyl, PhBBO, and (c) 1,3,5-trithia-2,4,6-triazapentalenyl, TTTA.

2. Methodology

The standard static First-Principles Bottom-Up (FPBU) procedure is applied to study the alleged magnetically isolated cuprate spin-ladder CuHpCl system, and to rationalize the magnetism of the purely organic PhBBO semiquinone-bridged bisdithiazolyl compound. Yet it is found that the study of the LT and HT phases of TTTA requires a dynamic perspective. Hereafter both working strategies will be described.

The first-principles bottom-up FPBU⁷ procedure implies four steps, as summarized herein. First, after inspection of the crystal structure, the symmetry-unique radical pairs that are likely to be magnetically relevant are identified (using a radical⋯radical distance cutoff value between spin-carrying moieties). Second, their magnetic exchange interactions, J_AB, are computed and the magnetic topology of the crystal (i.e., the network of connectivity defined by all relevant J_AB values) is defined. When necessary, wavefunction-based multireference CASSCF/PT2 and RASSCF/PT2 energies (E^LS, E^HS)^8,25 using a DZV basis set²⁶ are obtained from Molcas 7.6.²⁷ Otherwise energies at UB3LYP²⁸ or range-corrected CAM-B3LYP²⁹ level using 6-31+G(d) or 6-311++G(d,p) basis sets³⁰ are computed using Gaussian09³¹ for the dimer model and Orca 3.0³² package for the tetramer model. From the general Heisenberg Hamiltonian for a pair of S = 1/2 spin centres,

Ĥ = −2J_ABŜ_AŜ_B, (1)

the J_AB value is defined as 2J_AB = E^S − E^T, where E^S and E^T are the energies of the open-shell singlet and triplet states, respectively, of a two-radical cluster. In DFT calculations, the energy of the singlet state can be approximated using that of the single-determinant broken-symmetry (BS) solution.^5b,c Within this approximation, the expression chosen to compute the energy difference is^6dwhere E^S_BS is the energy of the BS solution^5b,c and S_ab is the overlap integral between the magnetic a and b orbitals of the BS solution. In our case, those orbitals are localized on each of the two radicals. This leads to S_ab = 0 and to the final expression that was used to compute J_AB values:

J_AB = E^S_BS − E^T. (3)

It should be noted that in the case of the tetramer model, there are (at least) four different J(di)'s among pairs of these four radicals (see Fig. 2). Therefore, the calculation of four different J(di)'s requires the evaluation of the energy of five different spin states. Within the tetramer approach, the spin states that have been evaluated are the high spin quintuplet (HS), two low spin triplet (LS1, LS2) and two low spin singlet (LS1/4, LS1/3) states (here numbers denote radicals with spin down in Fig. 2). Specifically Fig. 2 shows the set of equations used to calculate the non-negligible J_AB interaction between CuHpCl radicals in dimers d1, d3, d5, d6, and d7 (namely, J(d1), J(d3), J(d5), J(d6), and J(d7)).

Fig. 2 Tetramer model and set of equations used to calculate J(di) magnetic interactions between CuHpCl radicals (numbered 1 to 4). HS/LS stands for high spin and low spin, respectively.

Third, the Heisenberg Hamiltonian is applied to a model space (i.e. a subset of the magnetic topology), which is designed in such a way that, ideally, the resulting set of eigenvalues reproduces those that result from the application of the Heisenberg Hamiltonian to the full infinite crystal. Finally, the resulting energies and total spin numbers are introduced into the proper statistical mechanics expressions to calculate the macroscopic properties of the system, such as the magnetic susceptibility χT(T), heat capacity C_p(T) and magnetization M(H).²

The computational scheme adopted for the study of the interplay between thermal fluctuations and magnetism in TTTA consists of three steps: (i) ab initio molecular dynamics (AIMD) simulations³³ for both LT and HT phases of TTTA; (ii) computation of J_AB values between pairs of radicals for a large number of frames along the AIMD trajectories; and lastly, (iii) calculation of the vibrationally-averaged magnetic susceptibility [small chi, Greek, macron]_vib on the basis of full diagonalizations of the Heisenberg Hamiltonian built from the previously evaluated J_AB values. Accordingly, first of all, AIMD simulations³³ are performed at 300 K for both the LT and HT phases of TTTA (ca. 10 ps, time step 4 a.u.) as implemented in the CPMD package.³⁴ Note the temperature is chosen to be 300 K because TTTA is bistable at that temperature and there is crystallographic data available for both LT and HT phases. Supercells include 32 TTTA molecules (8 stacks of 4 radicals each). DFT calculations are carried out at PBE level³⁵ within the spin unrestricted formalism (broken symmetry singlet M_S = 0 state) using plane wave pseudopotentials³⁶ expanded at a kinetic energy cutoff of 25 Ry, together with Vanderbilt ultrasoft pseudopotentials,³⁷ a fictitious mass for the orbitals of 400 a.u., and Γ-point sampling of the Brillouin zone. The semiempirical dispersion potential introduced by Grimme,³⁸ in its DFT-D2 parameterization is also considered. In addition, the AIMD simulations are performed in the canonical (or NVT) ensemble using Nosé–Hoover chain thermostats.³⁹ Periodic boundary conditions in all three directions are imposed in all AIMD simulations. Pairs of radicals are then excised from the supercell of 32 radicals after AIMD simulations every 0.97 fs. Calculations of J_AB at UB3LYP/6-31+G(d)^28,30 level are conducted using dimer models: approximately 20 000 and 60 000 J_AB evaluations are carried out for the LT and HT polymorphs, respectively. Finally, the vibrationally-averaged magnetic susceptibility [small chi, Greek, macron]_vib for the HT phase at 300 K is computed by averaging the χ value over the whole set of configurations that are used to determine the time-evolution of the J_AB values between adjacent radicals within a stack. That is to say, the χ value is computed for an overall of ca. 10 000 different molecular configurations (each configuration was collected every 0.97 fs throughout the AIMD simulations) and, then, averaged. The resulting [small chi, Greek, macron]_vib at 300 K^14a is last compared to the χ value obtained from a static analysis²⁰ and from experiment.²⁴

3. Results and discussion

In order to stress the importance of the models and methods used to evaluate the J_AB magnetic coupling between radicals, three examples will be reported in which a too simple black-box study (based on static DFT calculations performed on dimers) would have prevented us from capturing the true nature of the magnetism of those three different molecule-based magnets.

3.1. On the size of the cluster model to evaluate J_AB magnetic coupling between radicals

The first example is the Cu₂(1,4-diazacycloheptane)₂Cl₄ molecular transition metal AFM complex (in short CuHpCl, see Fig. 1a). Since it was synthetized in 1997,²¹ CuHpCl has been the prototypical compound to quote when a magnetically isolated spin-ladder system was required due to its crystal packing consisting in well-separated ladder arrangements of copper atoms (see Fig. 3a). Therefore, from direct crystal observation, the magnetic topology consists of (apparently) non-connected spin-ladders with uniform rails (see Fig. 3b). Interestingly there is a large variety of experimental data, and therefore its study is very appealing. The overall AFM behaviour of CuHpCl has been characterized by measures of magnetic susceptibility, heat capacity, magnetisation, spin gap, inelastic neutron scattering, etc. Although the magnetic topology appears to be clear, according to literature, the fittings of experimental magnetic susceptibility data to different ladder models were not conclusive.⁴⁰ It thus follows that this copper derivative is more challenging than anticipated, as it must have many competing microscopic magnetic J_AB interactions. The alleged spin ladder magnetic topology of CuHpCl has been a long lasting issue.

Fig. 3 (a) [101] view of crystal packing of CuHpCl consisting of well-separated ladder arrangements of Cu atoms (highlighted as red arrows). (b) View of an isolated ladder showing intra-ladder first nearest neighbours nn pairs of radicals (d1–d5). Coloured lines have been added between Cu atoms to distinguish different radical⋯radical pairs. Evaluation of J(d3) (in green) using as a cluster model: (c) a bare d3 dimer, (d) a tetramer-based model that explicitly accounts for the radical d3 pair under study and the point charges of its d1 counterparts [dimer/d3-d1PC], (e) an eight radical model consisting on a two radical d3 pair embedded in six radicals represented by point charges [dimer/d3-d1PC]-4PC, and (f) a tetramer. Colour code: C (black), H (pale pink), N (light blue), Cu (blue) and Cl (green). PCs represented as shadowed atoms in blue.

For CuHpCl, all possible intra-ladder nearest-neighbour pairs of radicals were selected, namely d1, d2, d4 as rungs and d3, d5 as rails (see Fig. 3b). Once selected, all five corresponding J_AB magnetic couplings were evaluated using a dimer model. The only magnetically significant radical pairs were d1 (rung), and d3, d5 (rails) (see Table 1). Accordingly, the magnetic topology could be ascribed to consist of spin-ladders. Strikingly, the J_AB magnetic interaction for d3 and d5 (the rail interactions) has opposite sign. This magnetic behaviour is apparently at odds with their geometry, which would suggest that the spin-ladder should have uniform rails. Therefore, to be certain about the numerical value of J_AB for d3 and d5 pairs of radicals, the importance of the surrounding radicals was assessed by using different size cluster models with or without point charges (PC)⁴¹ representing the environment of the radicals whose magnetic strength was evaluated. For instance, regarding d3, we evaluated J(d3) using a dimer with its first nearest-neighbour replaced by point charges (namely dimer/d3-d1PC), using the previous model embedded in the point charges of its second nearest neighbour radicals (namely, [dimer/d3-d1PC]-4PC), and using a tetramer model (see Fig. 3c–f). The results show that with a tetramer model both rail interactions are AFM but not uniform at all (see Table 1). It is thus exceedingly important to stress that if the environment is not well described the coupling interactions are meaningless. For CuHpCl, it is also fundamental to realize that for a good description of each Cu monomer one has to account explicitly for the d1 pair, which becomes the magnetic building block of this compound and corresponds to a ferromagnetic FM interaction (see Fig. 4a in red, and Table 2). Each two of these dinuclear FM building blocks are then connected antiferromagnetically by d5 pairs of radicals (see Fig. 4a in green). It thus follows that the final 3D magnetic topology results from weaker AFM interactions (see Fig. 4a in orange; Fig. 4b shows a global view). Notice that the remaining non-negligible J(di) interactions are one order of magnitude smaller than d1, d5, d12 and d14 (see Table 2). Therefore, these magnetic couplings will not be important for simulation purposes. The resulting magnetic topology can be pictured as a 3D network of interacting squared plaquette magnetic building blocks. The simulation of the χ(T) magnetic susceptibility using J_AB magnetic interactions extracted from dimer and tetramer models further supports the failure of the dimer model since the numerical error involved in χ(T) when using the J_AB calculated at dimer level is larger than 200% compared to experiment, irrespective of the magnetic model used being an isolated spin-ladder or a 3D model of interacting spin-ladders (see full symbols in Fig. 4c). Yet the simulation using a 3D magnetic model of interacting plaquettes parameterized with J_AB's obtained from a 4-radical tetramer cluster model certainly agrees with the experimental data (see empty symbols in Fig. 4c). We can thus safely conclude that for CuHpCl the dimer model fails to evaluate J_AB magnetic couplings. Further we can also come to the conclusion that the resulting 3D magnetic topology is by far more complex than the alleged spin-ladder hinted from direct observation of the crystal packing (compare Fig. 3b and 4b).^9b

Table 1 Magnetic J(di) interactions (in cm⁻¹) for d1, d3 and d5 pairs of radicals (see Fig. 3b) using different cluster models, namely dimer, dimer with point charges ([dimer/di-d1PC] and [dimer/di-d1PC] – 4PC) and tetramer (see Fig. 3c–f for models to evaluate J(d3)). Note that X-ray crystallographic data at 4 K is used

Model	J(d1)	J(d3)	J(d5)
Dimer	+3.13	+2.11	−3.12
[dimer/di-d1PC]	+3.13	+0.08	−4.31
[dimer/di-d1PC] – 4PC	+2.02	+0.28	−3.68
Tetramer	+2.30	−0.37	−3.88

---

Fig. 4 (a) Schematic representation of two interacting squared-plaquette magnetic motifs. (b) Magnetic topology of CuHpCl, where only copper atoms are represented. Colour code: intra-plaquette AFM J(d5) green, and FM J(d1) red; inter-plaquettes AFM J(d12) orange and AFM J(d14) yellow. (c) Magnetic susceptibility curves for (i) J_AB values obtained with the dimer approach using one isolated spin-ladder magnetic model (full blue circles) and a 3D magnetic model consisting of three interacting ladder models (full red circles), and (ii) J_AB values obtained by the tetramer approach using a 3D magnetic model that includes four interacting plaquette-based magnetic building blocks (empty red circles). Note that experimental data is given for comparison purposes (full black circles).²¹

Table 2 Magnetic J(di) interactions (in cm⁻¹) between radicals using a tetramer model using X-ray crystallographic data at 4 K. J(di) classified according to crystal packing motif

J(di)/“intra”	J(d1)	J(d3)	J(d5)
	+2.30	−0.37	−3.88

---

J(di)/“inter”	J(d6)	J(d7)	J(d10)	J(d11)	J(d12)	J(d14)
	−0.29	−0.12	+0.22	−0.31	−1.05	−1.38

---

3.2. Assessing the adequacy of the method to evaluate J_AB magnetic coupling between radicals

The main objection in the CuHpCl case laid on the size of the dimer model, i.e. the dimer model was too small to evaluate correctly the J_AB microscopic interaction between two radicals. Changing gear, we will focus on the methods rather than on the models to evaluate J_AB magnetic coupling between radicals. The choice of the right method is also crucial and has to be assessed while studying the molecule-based material of interest. We would like to mention that CuHpCl was studied at UB3LYP/6-31+G(d) level,^28,30 which our experience grants as appropriate for a large variety of systems.^{9–11,14,16–20} We will now resort to the PhBBO magnet²² that belongs to the purely organic semiquinone-bridged bisdithiazolyl family of compounds (see Fig. 1b). In short, it has a 1D electronic structure due to the presence of π-stack motifs and paucity of close inter-column radical⋯radical contacts (see Fig. 5a). The experimental magnetic susceptibility χT(T) data suggests paramagnetic behaviour, with strong local FM interactions (see black line in Fig. 5c). This data was fitted to the Baker model for Heisenberg 1D FM chain of S = 1/2 centers resulting in J = +29.5 cm⁻¹ and a mean field parameter zJ′ = −2.5 cm⁻¹. PhBBO was also found to exhibit a phase transition at 4.5 K.

Fig. 5 (a) ab-View of the alternating ABAB π-stacks of PhBBO with intrastacks S3⋯S1 contacts d₁ = 3.68 Å (in black) and d₂ = 3.81 Å (in red) along the a-axis (π-stacking axis). (b) ac-View of Dimer1 along the a-axis (π-stacking axis). (c) Magnetic susceptibility as a function of temperature using J_AB calculated at unrestricted B3LYP (in red), CAM-B3LYP (in green), CASSCF(14,14) (in blue) and RASSCF(38,2,2;16,6,4) (in orange). Note that the line in black corresponds to experimental data.²²

In this case, DFT-based methods (such as unrestricted B3LYP²⁸ or range-corrected CAM-B3LYP²⁹) provide a good description of the spin density of an isolated radical (see monomer in Table 3). However, using a dimer model to evaluate the magnetic interaction between two radicals, DFT-based methods are able to correctly describe the spin density for the triplet state but fail to describe the open-shell singlet state (see Dimer1 in Fig. 5b and Table 3 for spin density). As a result, the FM interaction that we want to evaluate is overestimated at DFT level (+343.5 cm⁻¹ at UB3LYP, and +193.3 cm⁻¹ at CAM-B3LYP, see Table 4). We must then resort to wavefunction-based methods (such as CASSCF and RASSCF) in order to correctly describe the radicals spin density and, in turn, compute adequately the J_AB coupling interaction between them (see Tables 3 and 4). The J_AB magnetic interaction results to be +34.6 cm⁻¹ at CASSCF(14,14) and 49.7 cm⁻¹ at RASSCF(38,2,2;16,6,4) levels. Calculations at CASPT2 level are helpful to assess that dynamic correlation does not necessarily need to be considered to adequately evaluate J_AB magnetic interactions (see Table 4). The magnetic susceptibility is then calculated at all levels of theory. Clearly, unrestricted both B3LYP and CAM-B3LYP overestimate the value of χT at any given temperature (red and green lines in Fig. 5c). Contrarily, at CASSCF(14,14) level, χT is underestimated (blue line in Fig. 5c). One has to resort to RASSCF(38,2,2;16,6,4) to reproduce correctly the experimental data (orange line in Fig. 5c). For PhBBO, it can thus be concluded that DFT methodology dramatically fails.²³ Instead wavefunction-based methods are here required. The reason for the DFT failure turns out to be simple. It is because semiquinone-bridged bisdithiazolyl compounds are multireference systems due to the presence of low-lying open-shell states (what is called multi-orbital effect by experimentalists). This is thus the second example in which a black-box study based on standard DFT methods would have led to an incorrect interpretation, and to quantitatively-wrong values for the J_AB pair interactions of PhBBO.

Table 3 Mulliken spin densities of N atoms calculated for the doublet ground state of the PhBBO monomer, and for the triplet and broken symmetry (BS) open-shell singlet states of Dimer1 (see Fig. 5b for ac-view of geometry). The value of the spin density of the N atoms of the second radical within Dimer1 is indicated in parenthesis. The calculations were performed at unrestricted B3LYP and CAM-B3LYP levels using the 6-311++G(d,p) basis set

Method	System	Spin state	Spin density
			N1	N2
B3LYP	Monomer	Doublet	+0.27	+0.14
	Dimer1	Triplet	+0.27 (+0.27)	+0.13 (+0.13)
J(Dimer1) = +343.5 cm^−1		Singlet BS	+0.27 (−0.26)	+0.06 (−0.05)
CAM-B3LYP	Monomer	Doublet	+0.32	+0.21
	Dimer1	Triplet	+0.30 (+0.30)	+0.21 (+0.21)
J(Dimer1) = +193.3 cm^−1		Singlet BS	+0.25 (−0.24)	+0.16 (−0.16)

---

Table 4 Magnetic exchange couplings (in cm⁻¹) calculated for Dimer1 using the room temperature crystal structure. The calculations were performed with unrestricted B3LYP and CAM-B3LYP using 6-311++G(d,p) basis set, and with CASSCF and RASSCF using the DZV contraction of ANO-RCC basis set (note that active spaces are detailed). Also Mulliken spin densities of the two non-equivalent N atoms of Dimer1 at their triplet state are given

Basis	Method	J AB/cm^−1	Spin density
			N1	N2
6-311++g(d,p)	B3LYP	343.5	+0.27	+0.13
	CAM-B3LYP	193.3	+0.30	+0.21
DZV	CASSCF(14,14)	34.6	+0.30	+0.15
	RASSCF(38,2,2;16,6,4)	49.7	+0.29	+0.15
	RASPT2	48.7

---

3.3. Addressing thermal fluctuations to evaluate J_AB magnetic coupling between radicals

Finally, in addition to model and method, the J_AB couplings have to be able to adapt to changes in the magnetic topology due to thermal fluctuations. We have encountered that molecular dynamics simulations are necessary when studying bistable compounds whose radicals pack forming π-stacks along a given crystallographic direction, as in the case of 1,3,5-trithia-2,4,6-triazapentalenyl (TTTA, see Fig. 1c).¹⁴ Note that bistability is the ability of a material to present two stable phases that can both exist within a given range of temperatures but above and below that range only one or the other phase exists. TTTA was chosen because the neutral radical TTTA is the most studied molecule among dithiazolyl DTA radicals. Therefore, TTTA has become a prototypical molecule-based bistable material whose low temperature (LT, see Fig. 6a) phase is diamagnetic and high temperature (HT, see Fig. 6b) phase is paramagnetic.

Fig. 6 Packing of the monoclinic HT (SAXPOW05) and triclinic LT (SAXPOW06) crystal structures recorded at room temperature (RT): π stacks in the (a) LT (viewed along c) and (b) HT (viewed along b) phases. See inset values of J_AB magnetic couplings for (a) eclipsed and slipped pairs and (b) regular π-stacking pair. (c) Magnetic susceptibility χ(T) for LT and HT phases using structures determined at RT SAXPOW06 () and 250 K SAXPOW03 ()/RT SAXPOW05 () employing a 16-radical magnetic model. The experimental curve showing the thermal hysteresis present in the range 210–320 K, where HT and LT structures can both exist, is also given in black. All calculations were done at UB3LYP/Aug-cc-pVTZ level. (d) χ(T) plot with the approximate positions on hysteresis loop of each of the six published X-ray data sets of TTTA crystal (LT: recorded at 150 K and RT, SAXPOW01 and 06; HT: recorded at 225 K, 250 K, RT, and 310 K, SAXPOW04, -03, -05 and -02).²⁴ Refcodes SAXPOW01-06 given by the Cambridge Crystallographic Data Base CCDC.⁴²

Analysis of the crystallographic data at room temperature shows that the radicals pack forming π-stacks of eclipsed and slipped dimers in the LT phase (Fig. 6a) and regular π-stacks in the HT phase (Fig. 6b). The evaluation of J_AB magnetic coupling between radicals indicates that at LT the magnetic topology consists in AFM chains of strongly AFM dimers (−1755 cm⁻¹, see Table 5). In fact, it is because dimers have such strong AFM coupling that the LT phase is magnetically silent (see SAXPOW06 data in Fig. 6c). Contrarily, the HT phase shows a regular AFM chain topology (−135.6 cm⁻¹, see Table 5), which results in a paramagnetic response (see SAXPOW05 data in Fig. 6c). Since there are X-ray data at 4 different temperatures for the HT phase and at 2 different temperatures for the LT phase (see Fig. 6d), we wondered which would be the effect (if any) of the temperature in the evaluation of J_AB. New calculations showed that the magnetic topology was preserved as AFM π-stacks of TTTA radicals (see Table 5), irrespective of the temperature at which the crystal data was characterized. Note that basis set was also tested, and corroborated to not influence in the final J_AB calculated. The LT phase always behaves as magnetically silent due to the very large AFM interaction between radicals (see bottom simulation in Fig. 6c). In contrast, for HT, depending on the temperature, a shift in the value of the magnetic susceptibility was observed (see upper simulation in green and blue in Fig. 6c). That is, the dependence of J_AB magnetic coupling on the interplanar distance between radicals, which decreases as temperature lowers, causes the magnetic response to also be temperature-dependent. At this point, it is clear that the importance of thermal fluctuations must be quantified. Therefore, molecular dynamics simulations were carried out to be able to account for them in order to properly study TTTA.

Table 5 Values of J_AB interactions (in cm⁻¹) for the unique radical⋯radical pairs along the π-stacking for the HT and LT phases at UB3LYP level using 6-31+G(d) and Aug-cc-pVTZ basis sets. Note that the temperature at which J_AB is evaluated is also indicated (together with the CCDC code)

LT pair	Basis set	J AB/cm^−1	HT pair	Basis set	J AB/cm^−1
150 K SAXPOW01	6-31+G(d)	−2020.1	225 K SAXPOW04	6-31+G(d)	−175.0
	aug-cc-pVTZ	−2261.1		aug-cc-pVTZ	−199.9
T room SAXPOW06	6-31+G(d)	−1755.0	250 K SAXPOW04	6-31+G(d)	−160.6
	aug-cc-pVTZ	−1967.3		aug-cc-pVTZ	−183.8
			T room SAXPOW05	6-31+G(d)	−135.6
				aug-cc-pVTZ	−155.6
			310 K SAXPOW02	6-31+G(d)	−134.2
				aug-cc-pVTZ	−154.1

---

The time resolved evolution of the distance between centroids of adjacent TTTA radicals in one π-stack of the LT phase at 300 K shows that the arrangement of eclipsed dimers and slipped dimers is retained throughout the dynamics (see Fig. 7a). In contrast, the dynamics of the HT polymorph at 300 K shows that a given TTTA radical continually exchanges adjacent neighbour (upper or lower) with which it forms an eclipsed dimer (see Fig. 7b). In other words, our results indicate that each regular π-stacking of radicals in the HT phase is the resulting average of a fast pair exchange dynamics, characterized by a rapid interconversion between two degenerate distorted stacks (see Fig. 7c), which we have named pair exchange dynamics. This dynamic picture at the structural level leads to wondering whether it makes sense to use the common static approach to interpret the magnetism of TTTA, since each geometrical arrangement along the dynamics must involve a different magnetic coupling. Using selected geometries directly extracted from the molecular dynamics trajectories, the corresponding J_AB values of both polymorphs were evaluated (see Fig. 8a and b). Clearly the values of J_AB magnetic couplings feature remarkable large-amplitude fluctuations because the structure of the material undergoes large fluctuations over time. For LT, all evaluated J_AB are very large and negative along the entire trajectory (Fig. 8a), and lead to the overall diamagnetic behaviour of this polymorph. For HT, the effect of the thermal fluctuations is more important since the values of J_AB range from FM to highly AFM (see Fig. 8b and c for J_AB values in two random snapshots at 2.5 and 4.5 ps). Note that 1 ps involves the evaluation of 1000 J_AB interactions between pairs of radicals, i.e. the calculation of 1000 triplet and 1000 open-shell singlet states. In fact, if each set of these J_AB fluctuations is used to calculate their corresponding value of the magnetic susceptibility at 300 K and, then, the vibrationally-averaged magnetic susceptibility^14a is calculated, a substantial improvement is observed compared to the static data (see Fig. 8d). Specifically, one must compare an experimental χ(T) value of 4.5 × 10⁻⁴ emu mol⁻¹ to the computed magnetic susceptibility value using the static FPBU approach of 6.7 × 10⁻⁴ emu mol⁻¹ and the corresponding computed vibrationally average value of 5.2 × 10⁻⁴ emu mol⁻¹. As mentioned, the improvement obtained using a dynamical perspective is evident. Yet a diagnostic criterion is still missing and, thus, it is most important to assess in which systems molecular dynamics simulations are required.

Fig. 7 Time-resolved evolution of the distance between centroids of adjacent TTTA radicals in one of the columns (see d₁ to d₃ in a 4 radicals π-stack cluster model) of (a) LT-300 and (b) HT-300. (c) Regular π-stack of the HT polymorph (centre) results from the dynamic interconversion between two distorted stacks (left- and right-hand side). Examples of each distorted π-stack as extracted from the AIMD HT simulation at 300 K.

Fig. 8 Time-resolved fluctuations of (a) J(d₁) and J(d₃) for the LT polymorph, and (b) J(d₁), J(d₂) and J(d₃) for the HT polymorph. Each of these J(d_n) values corresponds to the pair of radicals marked with d_n in Fig. 7. (c) J_AB values and corresponding distances in two random snapshots at 2.5 and 4.5 ps. (d) Comparison between the values of the magnetic susceptibility calculated using the static FPBU approach (red) and its vibrationally averaged counterpart (green) at 300 K. Notice that the experimental temperature dependence of the magnetic susceptibility for TTTA on cooling (empty symbols) and on heating (full symbols) is given as a reference.

Overall, insight into the J_AB coupling between magnetic sites is crucial to understand the properties of magnetic molecule-based materials. Quantum chemical calculations are thus essential to shed light on the nature and mechanism of the magnetic coupling, namely whether it is through space/through ligand, direct exchange/kinetic exchange, spin polarization/spin delocalization, etc. Understanding and determining the nature and magnitude of the microscopic interactions and their dependence on electronic and structural factors is thus exceedingly important to make real advances in the fields of experimental characterization and modelling.

4. Conclusions

Albeit highly pursued, the total automatization of the modelling of the magnetism of molecular materials using their crystallographic data as only input faces significant challenges. A sufficiently complex automated black-box procedure should be able to address (i) whether the size of the cluster models to calculate the J_AB magnetic coupling is the most adequate, (ii) which is the most correct method to be used in the J_AB calculation (DFT vs. wavefunction-based, static vs. dynamic approach) and, finally, (iii) (although not discussed here) whether the chosen magnetic models that will be used for full diagonalisation purposes are the most suitable. Three examples have been carefully chosen to emphasize the importance of controlling the evaluation of the J_AB magnetic coupling to model and rationalize magnetism in molecule-based materials.

In the first case, a dimer-based automated procedure would have prevented uncovering the 3D nature of the magnetic topology of CuHpCl. Let us stress that this 3D magnetic topology is not hinted by J_AB values extracted from a dimer model nor could be inferred from direct observation of the crystal. In the second case, an overestimation of the FM interaction between radicals would have been obtained for PhBBO at the default unrestricted DFT level. Instead, wavefunction-based methods are required to adequately describe the presence of low-lying open-shell states in these PhBBO semiquinone-bridged bisdithiazolyl compounds. Finally, it has been shown that for compounds with π–π radical labile packing, such as TTTA, thermal fluctuations entail large amplitude variations in the value of J_AB's, which have significant effects in the magnetic response of the material. Effects that would have gone undetected under a total automatization of the modelling of the magnetism of TTTA using a static approach.

To sum up, in order to evaluate correctly J_AB exchange interactions, we need to control the size of the cluster models, and the correct method to be used. In addition to model and method, the J_AB couplings have to be able to adapt to changes in the magnetic topology due to thermal fluctuations. Therefore, building upon experience, in the field of modelling molecule-based magnetic systems, the “human know-how” cannot yet be replaced by a simple automated black-box procedure. Our take-home message is not that the design of an automated strategy is not attainable, but that needs to be incorporated in a sufficiently complex workflow. Considering the recent advances in Chemo-Informatics and Machine-Learning tools at complementing automatized workflows relevant to the task (e.g. oxidation state recognition,⁴³ active-space selection⁴⁴), we anticipate that these techniques will soon be able to ease part of the current working strategy.

Conflicts of interest

There are no conflicts to declare.

Note added after first publication

This article replaces the version published on 5^th July 2021, which contained errors in Tables 3 and 5.

Acknowledgements

MD, JRA, and JJN acknowledge financial support from MINECO (CTQ2017-87773-P/AEI/FEDER, UE), Spanish Structures Excellence María de Maeztu program (MDM-2017-0767), and Catalan DURSI (2017SGR348).

References

1. (a) E. Coronado, Nat. Rev. Mater., 2020, 5, 87; (b) A. E. Thorarinsdottir and T. D. Harris, Chem. Rev., 2020, 16, 8716; (c) Molecular Magnetic Materials: Concepts and Applications, ed. B. Sieklucka and D. Pinkowicz, Wiley-VCH, 2017; (d) S. M. Winter, S. Hill and R. T. Oakley, J. Am. Chem. Soc., 2015, 137, 3720; (e) J. S. Miller, Mater. Today, 2014, 16, 224; (f) I. Ratera and J. Veciana, Chem. Soc. Rev., 2012, 41, 303; (g) Molecular Materials, ed. D. W. Bruce, D. O’Hare and R. I. Walton, Wiley, 2011; (h) P. M. Lahti, Adv. Phys. Org. Chem., 2011, 45, 93; (i) Stable Radicals: Fundamentals and Applied Aspects of Odd-Electron Compounds, ed. R. Hicks, Wiley, New York, 2010.
2. (a) R. L. Carlin, Magnetochemistry, Springer-Verlag, Berlin Heidelberg, 1986; (b) O. Kahn, Molecular Magnetism, VCH Publishers, New York, 1993; (c) R. Boca, Theoretical Foundations of Molecular Magnetism, Current Methods in Inorganic Chemistry, Elsevier Science, Amsterdam, 1999, vol. 1.
3. (a) J. J. Borras-Almenar, J. M. Clemente-Juan, E. Coronado and B. S. Tsukerblat, J. Comput. Chem., 2001, 22, 985; (b) W. Duffy and K. P. Barr, Phys. Rev., 1968, 165, 647; (c) K. M. Diederix, H. W. J. Blote, J. P. Groen, T. O. Klaassen and N. J. Poulis, Phys. Rev. B: Condens. Matter Mater. Phys., 1979, 19, 420; (d) J. W. Hall, W. E. Marsh, R. R. Weller and W. E. Hatfield, Inorg. Chem., 1981, 20, 1033; (e) W. E. Hatfield, J. Appl. Phys., 1981, 52, 1985; (f) E. Coronado, M. Drillon, A. Fuertes, D. Beltran, A. Mosset and J. Galy, J. Am. Chem. Soc., 1986, 108, 900; (g) T. Barnes, E. Dagotto, J. Riera and E. S. Swanson, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 3196; (h) J. J. Borras-Almenar, E. Coronado, J. Curely, R. Georges and J. C. Gianduzzo, Inorg. Chem., 1994, 33, 5171.
4. (a) J. J. Novoa, M. Deumal and J. Jornet-Somoza, Chem. Soc. Rev., 2011, 40, 3182; (b) N. A. G. Bandeira, D. Maynau, V. Robert and B. Le Guennic, Inorg. Chem., 2013, 52, 7980; (c) J. P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf and N. Guihery, Chem. Rev., 2014, 114, 429; (d) T. Dresselhaus, S. Eusterwiemann, D. R. Matuschek, C. G. Daniliuc, O. Janka, R. Pottgen, A. Studera and J. Neugebauer, Phys. Chem. Chem. Phys., 2016, 18, 28262; (e) C. de Graaf and R. Broer, Magnetic Interactions in Molecules and Solids, Springer, 2016; (f) G. David, F. Wennmohs, F. Neese and N. Ferre, Inorg. Chem., 2018, 57, 12769.
5. (a) R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules, Oxford University Press, New York, 1989; (b) L. Noodleman, J. Chem. Phys., 1981, 74, 5737; (c) L. Noodleman and E. R. Davidson, Chem. Phys., 1986, 109, 131.
6. For a discussion on unrestricted DFT with broken symmetry, see: (a) L. Noodlemann and D. A. Case, Adv. Inorg. Chem., 1992, 38, 423; (b) L. Noodlemann, C. Y. Peng, D. A. Case and J. M. Mouesca, Coord. Chem. Rev., 1995, 144, 199; (c) E. Ruiz, P. Alemany, S. Alvarez and J. Cano, J. Am. Chem. Soc., 1997, 119, 1297; (d) R. Caballol, O. Castell, F. Illas, I. D. R. Moreira and J. P. Malrieu, J. Chem. Phys. A, 1997, 101, 7860; (e) E. Ruiz, J. Cano, S. Alvarez and P. Alemany, J. Comput. Chem., 1999, 20, 1391; (f) H. Nagao, M. Nishino, Y. Shigeta, T. Soda, Y. Kitagawa, T. Onishi, Y. Yoshioka and K. Yamaguchi, Coord. Chem. Rev., 2000, 198, 265; (g) J.-M. Mouesca, J. Chem. Phys., 2000, 113, 10505; (h) F. Illas, I. D. R. Moreira, C. de Graaf and V. Barone, Theor. Chem. Acc., 2000, 104, 265; (i) F. Neese, Coord. Chem. Rev., 2009, 253, 526.
7. M. Deumal, M. J. Bearpark, J. J. Novoa and M. A. Robb, J. Phys. Chem. A, 2002, 106, 1299.
8. (a) T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic-structure Theory, Wiley, Chichester, 2000; (b) C. de Graaf and R. Broer, Magnetic Interactions in Molecules and Solids, Springer, 2016.
9. (a) S. Vela, A. Sopena, J. Ribas-Arino, J. J. Novoa and M. Deumal, Chem. – Eur. J., 2014, 20, 7083; (b) J. Jornet-Somoza, F. Cosi, M. Fumanal and M. Deumal, Dalton Trans., 2021, 50, 1754.
10. (a) M. Deumal, G. Giorgi, M. A. Robb, M. M. Turnbull, C. P. Landee and J. J. Novoa, Eur. J. Inorg. Chem., 2005, 4697; (b) A. Shapiro, C. P. Landee, M. M. Turnbull, J. Jornet, M. Deumal, J. J. Novoa, M. A. Robb and W. Lewis, J. Am. Chem. Soc., 2007, 129, 952.
11. (a) S. Vela, J. Jornet-Somoza, M. M. Turnbull, R. Feyerherm, J. J. Novoa and M. Deumal, Inorg. Chem., 2013, 52, 12923; (b) J. C. Monroe, M. A. Carvajal, M. Deumal, C. P. Landee, M. Rademeyer and M. M. Turnbull, Inorg. Chem., 2020, 59, 6319.
12. (a) M. Fumanal, S. Vela, J. J. Novoa and J. Ribas-Arino, Chem. Commun., 2015, 51, 15776; (b) M. Fumanal, F. Mota, J. J. Novoa and J. Ribas-Arino, J. Am. Chem. Soc., 2015, 137, 12843.
13. (a) J. L. Brusso, O. P. Clements, R. C. Haddon, M. E. Itkis, A. A. Leitch, R. T. Oakley, R. W. Reed and J. F. Richardson, J. Am. Chem. Soc., 2004, 126, 8256; (b) J. L. Brusso, O. P. Clements, R. C. Haddon, M. E. Itkis, A. A. Leitch, R. T. Oakley, R. W. Reed and J. F. Richardson, J. Am. Chem. Soc., 2004, 126, 14692; (c) A. Alberola, D. J. Eisler, L. Harvey and J. M. Rawson, CrystEngComm, 2011, 13, 1794; (d) W. Fujita, K. Kikuchi and K. Awaga, Angew. Chem., Int. Ed., 2008, 47, 9480; (e) M. Nakano, R. Kishi, K. Yoneda, Y. Inoue, T. Inui, Y. Shigeta, T. Kubo and B. Champaigne, J. Phys. Chem. A, 2011, 115, 8767; (f) R. C. Haddon, A. Sarkar, S. K. Pal, X. Chi, M. E. Itkis and F. S. Tham, J. Am. Chem. Soc., 2008, 130, 13683; (g) K. Lekin, H. Phan, S. M. Winter, J. W. L. Wong, A. A. Leitch, D. Laniel, W. Yong, R. A. Secco, J. S. Tse, S. Desgreniers, P. A. Dube, M. Shatruk and R. T. Oakley, J. Am. Chem. Soc., 2014, 136, 8050.
14. (a) S. Vela, M. Deumal, M. Shiga, J. J. Novoa and J. Ribas-Arino, Chem. Sci., 2015, 6, 2371; (b) M. Fumanal, S. Vela, J. Ribas-Ariño and J. J. Novoa, Chem. – Asian J., 2014, 9, 3612.
15. All energy levels are essential in the computation of magnetic susceptibility data. Unsuccessful attempts to calculate magnetic susceptibility data: (i) used a finite-size cluster approach to compute the energy levels of only the four lowest-lying spin states for binuclear complexes ( K. Fink, C. Wang and V. Staemmler, Inorg. Chem., 1999, 38, 3847 ), and (ii) solved the Heisenberg Hamiltonian using Lanczos techniques on large (20–32 spins) clusters to obtain the lowest/highest energy levels and then assumed an analytical density of states for the middle eigenvalues ( C. S. Hellberg, W. E. Pickett, L. L. Boyer, H. T. Stokes and M. J. Mehl, J. Phys. Soc. Jpn., 1999, 68, 3489 ).
16. (a) M. Deumal, M. J. Bearpark, M. A. Robb, Y. Pontillon and J. J. Novoa, Chem. – Eur. J., 2004, 10, 6422; (b) J. Jornet-Somoza, M. Deumal, M. M. Turnbull and J. J. Novoa, Polyhedron, 2009, 28, 1965; (c) M. Deumal, J. M. Rawson, A. E. Goeta, J. A. K. Howard, R. C. B. Copley, M. A. Robb and J. J. Novoa, Chem. – Eur. J., 2010, 16, 2741.
17. J. Jornet-Somoza, M. Deumal, M. A. Robb, C. P. Landee, M. M. Turnbull, R. Feyerhem and J. J. Novoa, Inorg. Chem., 2010, 49, 1750.
18. (a) J. Jornet-Somoza, N. Codina-Castillo, M. Deumal, F. Mota and J. J. Novoa, Inorg. Chem., 2012, 51, 6315; (b) Ref. 11a.
19. S. Vela, M. Deumal, M. M. Turnbull and J. J. Novoa, Polyhedron, 2013, 52, 699.
20. C. S. Clarke, J. Jornet-Somoza, F. Mota, J. J. Novoa and M. Deumal, J. Am. Chem. Soc., 2010, 132, 17817.
21. (a) B. Chiari, O. Piovesana, T. Tarantelli and P. F. Zanazzi, Inorg. Chem., 1990, 29, 1172; (b) M. B. Stone, Y. Chen, J. Rittner, H. Yardimci, D. H. Reich, C. Broholm, D. V. Ferraris and T. Lectka, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65.
22. X. Yu, A. Mailman, P. A. Dube, A. Assouda and R. T. Oakley, Chem. Commun., 2011, 47, 4655.
23. M. Fumanal and M. Deumal, Phys. Chem. Chem. Phys., 2016, 18, 20738.
24. (a) W. Fujita and K. Awaga, Science, 1999, 286, 261; (b) G. D. McManus, J. M. Rawson, N. Feeder, J. van Dujin, E. J. L. McInnes, J. J. Novoa, R. Burriel, F. Palacio and P. Oliete, J. Mater. Chem., 2001, 11, 1992; (c) W. Fujita, K. Awaga, H. Matsuzaki and H. Okamoto, Phys. Rev. B: Condens. Matter Mater. Phys., 2002, 65, 064434.
25. P.-A. Malmqvist, A. Rendell and B. O. Roos, J. Phys. Chem., 1990, 94, 5477.
26. DZV contraction of the ANO-RCC basis set using (3s2p1d) for N, O and C atoms, (2s) for H, and (4s3p1d) for S atoms: (a) P.-O. Widmark, P.-A. Malmqvist and B. O. Roos, Theor. Chim. Acta, 1990, 77, 291; (b) B. O. Roos, R. Lindh, P.-A. Malmqvist, V. Veryazov and P.-O. Widmark, J. Phys. Chem. A, 2004, 108, 2851.
27. F. Aquilante, L. De Vico, N. Ferre, G. Ghigo, P.-A. Malmqvist, P. Neogrady, T. B. Pedersen, M. Pitonak, M. Reiher, B. O. Roos, L. Serrano-Andres, M. Urban, V. Veryazov and R. Lindh, J. Comput. Chem., 2010, 31, 224.
28. (a) A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098; (b) A. D. Becke, J. Chem. Phys., 1993, 98, 5648; (c) C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785.
29. T. Yanai, D. Tew and N. Handy, Chem. Phys. Lett., 2004, 393, 51.
30. (a) P. C. Hariharan and J. A. Pople, Theor. Chim. Acta, 1973, 28, 213; (b) M. M. Francl, W. J. Petro, W. J. Hehre, J. S. Binkley, M. S. Gordon, D. J. DeFrees and J. A. Pople, J. Chem. Phys., 1982, 77, 365.
31. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. Bloino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery Jr, J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V. N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski and D. J. Fox, Gaussian 09, Revision D.01, Gaussian, Inc., Wallingford, CT, 2009.
32. F. Neese, ORCA-An ab initio, Density Functional and Semi-empirical Program Package, Version 3.0, Max Planck Institute for Bioinorganic Chamistry, Mullheim, 2012.
33. (a) D. Marx and J. Hutter, Ab Initio Molecular Dynamics: Basic Theory and Advanced Methods, Cambridge University Press, 2009; (b) R. Car and M. Parrinello, Phys. Rev. Lett., 1985, 55, 2471.
34. CPMD, Version 3.13.2; copyright IBM Corp., 1990–2009, copyright MPI fur Festkorperforschung Stuttgart, 1997–2001; http://www.cpmd.org.
35. (a) J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865; (b) J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
36. E. G. Parr and W. Yang, Density Functional Theory, Oxford University Press, New York, 1989.
37. D. Vanderbilt, Phys. Rev. B: Condens. Matter Mater. Phys., 1990, 41, 7892.
38. S. Grimme, J. Comput. Chem., 2006, 27, 1787.
39. G. J. Martyna, M. L. Klein and M. Tuckerman, J. Chem. Phys., 1992, 97, 2635.
40. (a) G. Chaboussant, M.-H. Julien, Y. Fagot-Revurat, M. Hanson, L. P. Lévy, C. Berthier, M. Horvatic and O. Piovesana, Eur. Phys. J. B, 1998, 6, 167; (b) G. Chaboussant, M.-H. Julien, Y. Fagot-Revurat, L. P. Levy, C. Berthier, M. Horvatic and O. Piovesana, Phys. Rev. Lett., 1997, 79, 925; (c) G. Chaboussant, Y. Fagot-Revurat, M.-H. Julien, M. E. Hanson, C. Berthier, M. Horvatic, L. P. Levy and O. Piovesana, Phys. Rev. Lett., 1998, 80, 2713; (d) G. Chaboussant, M.-H. Julien, Y. Fagot-Revurat, H. Mayaffre, M. Horvatic, L. P. Levy, C. Berthier and O. Piovesana, Phys. B, 2000, 280, 315; (e) P. R. Hammar and D. Reich, J. Appl. Phys., 1996, 79, 5392; (f) P. R. Hammar, D. H. Reich, C. Broholm and F. Trouw, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 7846; (g) X. Wang and L. Yu, Phys. Rev. Lett., 2000, 84, 5399; (h) W. Zheng, R. R. P. Singh and J. Oitmaa, Phys. Rev. B: Condens. Matter Mater. Phys., 1997, 55, 8052.
41. Point charges (PC) used are fitted to reproduce the electrostatic potential of the isolated cations according to Merz–Singht–Kollman scheme: (a) B. H. Besler, K. M. Merz Jr. and P. A. Kollman, J. Comput. Chem., 1990, 11, 431; (b) U. C. Singh and P. A. Kollman, J. Comput. Chem., 1984, 5, 129.
42. CCDC: (a) F. H. Allen, Acta Crystallogr., Sect. B: Struct. Sci., 2002, B58, 380; (b) F. H. Allen and W. D. S. Motherwell, Acta Crystallogr., Sect. B: Struct. Sci., 2002, B58, 407.
43. K. M. Jablonka, D. Ongari, S. M. Moosavi and B. Smit, 2020. https://chemrxiv.org/articles/preprint/Using_Collective_Knowledge_to_Assign_Oxidation_States/11604129/1.
44. (a) C. J. Stein and M. Reiher, J. Chem. Theory Comput., 2016, 12, 1760–1771; (b) C. J. Stein and M. Reiher, J. Comput. Chem., 2019, 40, 2216.

---