The quest for rationalizing the magnetism in purely organic semiquinone-bridged bisdithiazolyl molecular magnets †

Semiquinone-bridged bisdithiazolyl-based radicals (XBBO) are appealing purely organic magnetic building blocks for the synthesis of new functional materials. Remarkably, for the phenyl-derivative PhBBO, the rationalization of its magnetism becomes a proof of concept that DFT can dramatically fail to evaluate J AB magnetic interactions between purely organic radical pairs. Instead, wavefunction-based methods are required. Once J AB ’s are fully characterized, the magnetic topology of PhBBO is disclosed to consist of ferromagnetic FM p -stacks that are very weakly coupled (by FM and AFM J AB interactions). The magnetic susceptibility w T ( T ) and magnetization M ( H ) of PhBBO are then calculated using a first-principles bottom-up approach. The study of the unit cell contraction upon cooling from room temperature to zero-Kelvin is relevant to propose a suitable model for the phase transition that occurs at 4.5 K. A simplistic picture tells us that the antiparallel-aligned 1D-FM-chains convert into domains of weakly either FM- or AFM-coupled 1D-FM-chains. Accordingly, the presence of these domains may introduce geometrical spin frustration below 4.5 K.


Introduction
In an attempt to meet the electronic requirements for magnetic performance, namely, spin delocalization and orthogonal orbital overlap, heavy-heteroatom radicals have been widely studied. 1 Interest in the use of thiazyl radicals as building blocks for magnetic and conductive materials has thus grown rapidly, in part because the heavy heteroatom (sulfur) enhances intermolecular magnetic and electronic interactions. N-Alkylated pyridinebridged bisdithiazolyl radicals ( Fig. 1a) have emerged as a new versatile class of radicals for the synthesis of new functional organic materials. 2 However, they tend to pack in herringbone arrays of slipped p-stacks, which diminish the orbital overlap along the stacking direction. 2,3 The herringbone motif could be avoided using radicals based on N-alkylated pyrazine-bridged bisdithiazolyl (Fig. 1b). 4 Yet, dimerization is favored, although it can be suppressed by using the right ligands R 1 . 5 In addition to this synthetic strategy, an alternative approach consists of the isoelectronic replacement of the NR 1 moiety in 1a with a carbonyl group to result in a resonance-stabilized semiquinone-bridged bisdithiazolyl radical (denoted XBBO, see Fig. 1c). This replacement intended to be a simple isolobal exchange, 6 since the incorporation of the carbonyl, CO, group does not perturb the SOMO (see Fig. 1d) with respect to either pyridine-or pyrazine-bridged radicals due to the presence of a nodal plane. However, there is a significant mixing of the CO p*-orbital in the LUMO that results in low-lying open-shell states, whose presence, ultimately, has a major impact on both magnetic and charge-transport properties. 7 Accordingly, although this new class of radicals has a common molecular framework, the implementation of different exocyclic groups (X = H, 7c,8 F, 7b,9 Me, 10 Ph, 1a,11 Cl 1a ) is responsible for XBBO displaying a variety of structural, magnetic, conducting and electronic properties. In general, for XBBO crystals, p-stacked bis-1,2,3-dithiazolyl radicals apparently have good orbital overlap along the stacking direction. Within this context, the phenyl-substituted radical based on the semiquinonebridged skeleton (Fig. 1c with an exocyclic group X = Ph 11 ) was the first example of this new class of radicals, and is an appealing target for a computational study.
The carbonyl oxygen atom in XBBO plays an active structuremaking role, giving rise to tightly packed structures linked by short intermolecular SÁ Á ÁO contacts. These supramolecular synthons afford a rich variety of novel packing motifs, including the head-over-tail motif found for PhBBO (see Fig. 2). The experimental crystal data show that PhBBO crystallizes in the high-symmetry orthorhombic chiral space group P2 1 2 1 2 1 , whose unit cell contains two non-equivalent columns of stacked radicals organized in a head-to-tail disposition. 11 Sheets (bc-plane) of approximately coplanar radicals are laced together by a series of intermolecular S3Á Á ÁO contacts (2.858 Å) (see Fig. 2a). These bc-sheets are packed in layers to produce alternating ABABAB p-stacks (adjacent radicals are related by the 2-fold screw axis along a) (see Fig. 2b). From an experimental point of view, both the alternating p-stack motif and the paucity of close intercolumn radical-radical contacts suggested an approximately 1D electronic structure. 11 The available experimental wT(T) data at H = 1 kOe over the range 2-300 K suggested paramagnetic behavior, with strong local ferromagnetic (FM) couplings along the p-stacks. 11 These suggested strong FM interactions could be rationalized in terms of the presence of the above-mentioned head-over-tail p-stacking of PhBBO radicals. For temperatures above 100 K, the data are found to be fitted to a Curie-Weiss expression with C = 0.349 emu K mol À1 and y = +32.8 K. Within the 6-30 K temperature range, the Baker model for a Heisenberg 1D FM-coupled chain (p-stack) of S = 1 2 centers was used to fit the experimental data with a FM J chain (+29.5 cm À1 ) and an antiferromagnetic (AFM) mean field interaction zJ 0 (À2.5 cm À1 ) fitting parameters. According to ZFC-FC experiments at H = 100 Oe, a phase transition occurs at 4.5 K. It was experimentally proposed that below 4.5 K the FM chains become aligned in an antiparallel fashion to form a spincanted AFM (SC-AFM) state. The AFM ordering was observed to collapse when increasing the magnetic field and the system was experimentally suggested to undergo a spin-flop transition. 11 All the above reported experimental data motivate the need for further understanding of the electronic structure of the PhBBO material.
The standard approach for both the interpretation and the simulation of magnetic properties of materials is based on a static perspective. Accordingly, it is assumed that these properties can be obtained with a single static nuclear configuration (usually an X-ray resolved structure or, alternatively, an optimized structure). The magnetic response of the purely organic PhBBO molecular crystal is rationalized here on the basis of a first-principles bottom-up FPBU procedure. 12 One first analyzes the crystal packing to select radical pairs that might be magnetically coupled. Quantum chemistry methods are then used to evaluate all J AB magnetic couplings. The magnetically significant J AB interactions define the magnetic topology of the crystal. Full diagonalization techniques are then applied to a properly chosen magnetic model to generate the corresponding energy spectrum. Statistical mechanic tools are fed with those energy levels to simulate Fig. 1 Structures of (a) an N-Alkylated pyridine-bridged bisdithiazolyl radical, (b) an N-Alkylated pyrazine-bridged bisdithiazolyl radical, (c) a semiquinonebridged bisdithiazolyl radical XBBO and (d) the SOMO of the PhBBO radical (X = phenyl) at the B3LYP level (visualized with an isovalue of 0.02 a.u.). the macroscopic properties of interest (w, C p , M), which in turn will be compared to the available experimental data.
Summarizing, this work aims at discriminating whether the magnetic topology of PhBBO is mostly one-dimensional, with FM interactions within p-stacks and AFM interactions between p-stacks, by evaluating the magnetic coupling between radicals. The magnetic susceptibility and magnetization of PhBBO will then be calculated in order to validate the computed magnetic topology. Finally, the study of the unit cell contraction upon cooling from room temperature to 0 Kelvin targets at proposing a suitable model for the phase transition that occurs at 4.5 K.

Methodological details
The first-principles bottom-up FPBU procedure was applied to study the magnetism of the target PhBBO system. 12 This working strategy implies the completion of 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 PhBBOÁ Á ÁPhBBO distance cutoff value of 15 Å). 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) arises. Third, the Heisenberg Hamiltonian is applied to a model space (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 wT(T) and magnetization M(H).
The evaluation of the J AB magnetic coupling between purely organic PhBBO radicals was performed using a dimer and a larger tetramer cluster models. At the DFT level, low spin (LS) broken symmetry (BS) 13 and high spin (HS) single point energy calculations were performed using Gaussian09 14 for the dimer model and Orca 3.0 15 package for the tetramer model. WFT based multireference CASSCF/PT2 and RASSCF/PT2 energies (E LS , E HS ) 16 were obtained from Molcas 7.6. 17 The evaluation of J AB interactions by DFT was performed using the 6-311++g(d,p) basis set, 18 and by WFT methods was performed using the following contraction of the ANO-RCC basis set: 19 (3s2p1d) for N, O and C atoms, (2s) for H, and (4s3p1d) for S atoms. For the latter, the Cholesky decomposition was used to reduce the computational cost associated with the calculation of the two-electron integrals. 20 When performing CASPT2 and RASPT2, the standard zero-order Hamiltonian with an IPEA 21 shift of 0.25 a.u. was applied in the definition of the diagonal Fock matrix elements of the active orbitals. The definition of the different active spaces used is given in ESI, † Section S1. All self-consistent field (SCF) energy convergences have been set up at 10 À7 a.u., thereby allowing an accuracy of 0.05 cm À1 in the evaluation of the magnetic exchange J AB values.
The unit cell contraction upon cooling from room temperature to zero Kelvin was studied employing plane wave pseudopotential calculations for the variable-cell geometry optimization. This calculation was carried out using the PBE exchangecorrelation functional 22 within the spin-unrestricted formalism, together with Vanderbilt ultrasoft pseudopotentials, 23 and G-point sampling of the Brillouin zone. In the calculation, the semiempirical dispersion potential introduced by Grimme 24 was added to the conventional Kohn-Sham DFT energy in order to properly describe the van der Waals interactions among all different PhBBO radicals (the so-called DFT-D2). Notice that good predictions for the structure and cohesive energies of molecular crystals are obtained by the use of PBE together with the D2 Grimme correction. 25 Our variable-cell optimization, in which the atomic positions and the lattice parameters are optimized simultaneously, was carried out using the QUANTUM ESPRESSO package 26 using a plane wave basis set expanded at a kinetic energy cutoff of 70 Ry. The initial atomic positions and lattice parameters for the optimization were taken from the X-ray resolved structures of PhBBO at room temperature (296 K). 11 The initial lattice parameters of the supercell are collected in ESI, † Section S2.

Results and discussion
The nature of magnetism of the PhBBO crystal: DFT failure reviewed from a WFT perspective The unit cell of the PhBBO crystal contains two non-equivalent columns of stacked radicals in a head-to-tail disposition: there is one symmetry unique intra-column p-pair (namely Dimer1 along the a-axis) and six different inter-column pairs (Pdimer1-Pdimer6) within a 15.0 Å radius cutoff (Fig. 3a). Pdimer6 and Pdimer1 pairs run along the b-axis forming zigzag chains that are connected in such a way that both dimers alternate along the c-axis (see Fig. 3b). These bc-layers are then laced together by Pdimer2, Pdimer3, Pdimers4 and Pdimer5 contacts (see Fig. 3c). The network of intermolecular interactions within the PhBBO crystal is thus particularly complex (see Fig. 3 for a complete picture of all the intra-and inter-stack radical pairs that have potential to result in non-negligible J AB magnetic couplings). Accordingly, there is a J AB magnetic interaction associated with each radicalÁ Á Áradical contact that has to be evaluated (see the color code in captions of Fig. 3 and 4 for radical pair geometries).
The well-known B3LYP hybrid functional, 27 which is acknowledged to perform well to compute J AB magnetic interactions between purely organic radicals, 28a,29,30 and the range-corrected CAM-B3LYP 31 methods are used to evaluate the J AB magnetic coupling of the previously selected seven pairs of PhBBO radicals (Fig. 4). The nuclear coordinates were extracted from X-ray data at room temperature. 11 The Hamiltonian used is Ĥ = À2 P J AB Ŝ A ÁŜ B . Thus, for a pair of PhBBO radicals showing a HS triplet state and a LS open-shell BS singlet state, J AB = E BS,LS À E HS at the DFT level within the broken symmetry, BS, approximation. 13 Table 1 shows all calculated J AB magnetic interactions using the 6-311++G(d,p) basis set. 18 The results obtained indicate that, among all pairs, only Dimer1, Pdimer1, Pdimer2 and Pdimer6 provide a non-zero J AB . By comparison with the fitted J values (+29.5, À2.5 cm À1 ), it can be inferred that the values of J AB for Dimer1 (+343.5/+193.3 cm À1 ) and Pdimer6 (+24.6/+5.9 cm À1 ) are highly overestimated at the DFT level.
An in-depth study of Dimer1 (and Pdimer6) is required in order to disentangle the disagreement between calculated J AB at the DFT level and fitted J values. While DFT calculations are known to be sometimes insufficient to properly describe inorganic molecule-based magnets, 32 our experience indicates that DFT performs correctly with purely organic radical moleculebased magnets. 29 As spin density in the PhBBO radical is delocalized over the whole p-system (see the SOMO in Fig. 1d), we focused on the spin density on the nitrogen atoms to monitor whether the radical pair is well described or not. 33 The values obtained in each case are reported in Table 2. For B3LYP and CAM-B3LYP, it can be seen that the spin densities on the N atoms obtained for the monomer are in agreement with the literature (spin densities on N1/N2 are 0.27/0.14 and 0.32/0.21, respectively). 1a,7a,c Regarding Dimer1 in its triplet state, the values again match those obtained for the monomer indicating that the spin distribution is well described. However, the spin density distribution of Dimer1 in the open-shell BS singlet differs from the latter two, since smaller values are obtained for the N2 atom. This can clearly indicate a problem in the spin density distribution of this electronic state. The problem with the open-shell BS singlet state could be related to the radical pair model used to compute the microscopic J AB magnetic interaction. The effect of the cooperativity between the closest neighboring radicals on the electronic density (and, in turn, on the J AB values) is evaluated by enlarging 34 the pair magnetic model to a tetramer cluster model to include both Dimer1 and Pdimer6 interactions. Using the tetramer model, J AB (Dimer1) coupling values are calculated to be +336.0 cm À1 and +181.2 cm À1 at B3LYP and CAM-B3LYP levels, respectively (see ESI, † Section S3, for more details). Therefore, there is no J AB improvement using a larger cluster model. Neither spin-orbit coupling nor anisotropy can be the reason for the DFT failure since a g value of 2.009 was reported for the PhBBO radical. 11 We have also evaluated the possible enhancement of the magnetic interactions due to hydrogen interactions by using larger basis sets, 35 but with the same outcome. Finally, in PhBBO there are neither counterions nor solvent to be accounted for. Therefore, so far the environment of the radicals has also been properly described. 36 Apparently there is no explanation for the failure of B3LYP or CAM-B3LYP functionals when trying to describe the magnetic interaction between Dimer1 (and Pdimer6) semiquinone-bisdithiazolyl radicals. Note that calculations done at the wB97XD level using 6-31+g(d) and 6-311++g(d,p) basis sets give +179.2 and +180.5 cm À1 as J AB values, respectively, in agreement with our previous CAM-B3LYP results, which was expected since both functionals are range separated hybrids.
We then resorted to WFT based methods, such as CASSCF and RASSCF, 16 to evaluate all J AB interactions. At the WFT level, J AB = (E LS À E HS )/2 since the Hamiltonian used is the same as for the DFT level (Ĥ = À2 P J AB Ŝ A ÁŜ B ). Several CASSCF and RASSCF spaces (including p-type orbitals and using a DZV basis set, 19 see ESI, † Section S1) were first tested at the X-ray nuclear coordinates of the PhBBO monomer (Table 3), the spin density being the criterion to decide whether the active space was adequate. 33 Indeed, the spin density of the N atoms varies depending on the active space considered. According to Table 3, CASSCF(7,7) and RASSCF(7,2,2;3,1,3) are required for a good description of the PhBBO radical (note that the RASSCF restricted active space is defined as follows: the number of electrons in all RAS space, Fig. 4 Molecular geometry of the intra-p-stack Dimer1, and inter-stack Pdimer1-6 radical pairs, for which the J AB magnetic exchange values have been evaluated. Notice that the two PhBBO radicals in Pdimer1 and Pdimer6 lay coplanar in the bc-plane. Table 1 Magnetic exchange couplings calculated for the seven pairs of radicals in the room temperature crystal structure. The calculations were performed at B3LYP and CAM-B3LYP level using the 6-311++g(d,p) basis set, and at CASSCF and RASSCF level using the DZV contraction of the ANO-RCC basis set and different active spaces (see main text and ESI, Section S1). The OÁ Á ÁO distance between each radical pair is also indicated. See the 'Methodological details' section for a description on the basis sets    the maximum number of holes in RAS1, and the maximum number of particles in RAS3; orbitals in RAS1, orbitals in RAS2, and orbitals in RAS3). Correspondingly, equivalent CASSCF and RASSCF spaces were tested for Dimer1 and Pdimer1 using two different basis sets ( Table 4). The results show that the minimal active spaces are not enough; one must resort to CASSCF (14,14), RASSCF(38,2,2;18,2,6) and RASSCF(38,2,2;16,6,4) to properly describe the spin density. Besides, the DZV to TZV basis set effect is small (see Table 4; for Dimer1: +34.6 to +36.4 cm À1 at the CASSCF level and +44.6 to +43.1 cm À1 at the RASSCF level; for Pdimer1: À0.4 to À0.5 cm À1 at the RASSCF level). Thus, DZV contraction was chosen for the calculation of the J AB interaction for all remaining radical pairs (Pdimer2-Pdimer6). For Dimer1, we also corroborated that the value of J AB changes by ca. 8% with the inclusion of dynamic correlation: J AB = +44.6 cm À1 at the RASSCF(38,2,2;18,2,6) level and J AB = +48.7 cm À1 at the RASPT2 level. Since RASPT2 calculations are very demanding in terms of computational cost and they imply a small relative improvement (8%) of the J AB values, our best options to pursue our computational study on the PhBBO crystal are the use of CASSCF (14,14), RASSCF(38,2,2;18,2,6) and RASSCF(38,2,2;16,6,4) (hereafter called CASSCF, RASSCF-1, and RASSCF-2, respectively). Table 1 shows the J AB values obtained with these active spaces for all unique radical pairs identified in the PhBBO crystal. Notice that, as suggested in the literature, 7b our calculations corroborate that crystals of PhBBO present FM interactions, as also the other members of the XBBO oxobenzobridged bisdithiazolyl family do, 7-11 and become a proof of concept that DFT can dramatically fail to evaluate J AB magnetic interactions between purely organic radical XBBO pairs. The B3LYP failure was then reviewed from a WFT perspective by comparing the TTTA (1,3,5-trithia-2,4,6-triazapentalenyl) and PhBBO crystals (see Fig. 5). TTTA has been recently shown to be well described at the B3LYP level. 28a Specifically, the high temperature phase of TTTA has a p-stacked pair of radicals whose J AB value is À143 cm À1 at the B3LYP/6-311++g(d,p) level (see Table 5 for (TTTA) 2 ). According to Table 5, for the TTTA pair, the WFT based methods require the inclusion of all p-orbitals in the active space as well as dynamic correlation to describe properly the magnetic interaction (J AB = À162 cm À1 at the CASPT2 level). Notice that this value is in agreement with the result obtained at the B3LYP level, and with the value of À135 cm À1 previously reported at the DDCI level. 37,38 For Dimer1, we have already discussed that, although it is also required to include all p-orbitals in the active space (RASSCF-1), dynamic correlation (RASPT2) hardly affects the resulting J AB value. However, in the current case, for Dimer1 there is no agreement between J AB values computed at DFT and WFT levels. The reason for B3LYP performing correctly for the TTTA radical and (TTTA) 2 but not for the PhBBO radical and Dimer1 is, basically, that the latter two show a larger multireference character than the TTTA counterparts (see ESI, † Section S4, for natural molecular orbitals and occupation numbers). Referring to Dimer1, the occupation numbers different from 2.0 vary from 1.83 to 1.94 for doubly occupied natural molecular orbitals and, consequently, they range from 0.18 to 0.06 for virtual natural molecular orbitals. Contrarily, for (TTTA) 2 , doubly occupied natural molecular orbitals have occupation numbers (other than 2.0) ranging from 1.90 to 1.94. This means that (TTTA) 2 can be adequately described using a monoreference method, such as DFT, while one has to resort to multireference methods for Dimer1, such as CASSCF and RASSCF. Indeed, this computational evidence is acknowledged in the literature as the multi-orbital effect by experimentalists themselves. 1b,c Analogous to the TTTA dimer, the exchange coupling between N-alkylated (pyridine/pyrazine)-bridged bisdithiazolyl pairs  of radicals (see Fig. 1a and b) can be evaluated using DFT-based methods since there is no multi-orbital effect involved. Indeed TTTA and N-alkylated (pyridine/pyrazine)-bridged bisdithiazolyl radicals typify conventional single orbital radicals. 1c Once more it is obvious that DFT-based methods must be used with care depending on the features of the studied system since they might be of no validity 7b,c,10 even for purely organic magnets and lead to misinterpretations.

Magnetic susceptibility vT(T) and magnetization M(H) simulations of the PhBBO crystal
The computed J AB magnetic interactions (see Table 1 for J AB values) show that at room temperature the magnetic topology is three-dimensional (3D), and consists of FM p-stacks that are very weakly coupled (both FM and AFM). We have corroborated that there are FM interactions within p-stacks (Dimer1), as suggested experimentally based on the fitting Baker model (+29.5 cm À1 ). One can safely say that the head-over-tail p-stacking of PhBBO radicals found in the crystal is responsible for remarkably strong FM interactions along the p-stacks. The main magnetic motif is thus 1D, but there are small AFM (Pdimer1) as well as FM (Pdimer2 & Pdimer6) interactions between p-stacks (see Table 1 for J AB values). Here we must stress the fact that FM interactions between p-stacks were experimentally not anticipated, since the fitted cumulative inter-stack (mean field) parameter was negative (À2.5 cm À1 ). Therefore, our results suggest that PhBBO cannot be merely viewed as a spin-canted AFM (SC-AFM) state as proposed experimentally, since the inter-stack interactions are not only AFM in nature. However, the schematic and simplest picture of the a and b spin distribution that emerges from our RASSCF-1 computed magnetic topology (see Fig. 6 and ESI, † Section S5, for full discussion on the choice of J RASSCF-1 AB values) shows an overall antiparallel alignment of FM chains. This is why it appears to be in a SC-AFM state. From Fig. 6, one realizes that the magnetic topology at room temperature introduces certain small degree of geometrical spin frustration between the radicals connected by Pdimer2 (+0.1 cm À1 ), which tend to ferromagnetically couple the FM p-stacks (+44.6 cm À1 , Dimer1) of the PhBBO crystal, otherwise antiferromagnetically coupled by Pdimer1 (À0.4 cm À1 ).
In order to explore the dimensionality required by the magnetic model for simulation purposes, a series of wT(T) calculations were performed using the J AB interactions at the RASSCF-1 level (Fig. 7). The inter-column magnetic interactions were studied using two distinct models: 2 Â 4 and 4 Â 4 models (see Fig. 7a). It is clear that the wT(T) calculated data are well converged using 2 columns, i.e. the 2 Â 4 model (Fig. 7c), since the simulated wT(T) data using either a 2 Â 4 or a 4 Â 4 model practically overlap. Also, different models were used to check the dependency of the wT(T) calculated data on the magnetic model as it is enlarged along the p-stacking direction (namely, 2 Â 4, 2 Â 6 and 2 Â 8 models in Fig. 7b). According to Fig. 7c, we can conclude that the magnetic interactions along the p-stacking direction are reasonably well Fig. 6 Schematic representation of the magnetic topology at room temperature obtained for the PhBBO crystal at the RASSCF-1 level. Color code for the significant J AB magnetic coupling: Dimer1 (dark blue, FM), Pdimer1 (orange, AFM), and Pdimer2 (green, FM). The arrows represent the spin propagation according to the J AB magnetic interaction values to provide a schematic picture of the a and b spin distribution. Fig. 7 Representation in the (a) bc-and (b) ac-planes of the magnetic exchange models employed in the computation of the wT vs. T curves shown in (c) at an applied field of 1 kOe. The experimental data are also given for comparison purposes. described at the 2 Â 8 level, i.e. using 16 PhBBO radicals. Indeed, we have confirmed that even at low temperatures the computed w(T) data agree reasonably well with the available experimental w(T) data (see ESI, † Section S6). Due to the full diagonalization code we use, we are limited to models up to N = 18 S = 1 2 radical centers for computer efficiency, the memory required being N!/[(N/2)!(N/2)!]. Although a larger model will reproduce better the numerical wT(T) values, the 2 Â 8 model is our best option. Therefore, from this point on, the 2 Â 8 magnetic model will be used in all simulations.
Let us now focus on the phase transition that the PhBBO crystal undergoes at 4.5 K. It was first suggested that, after the phase transition, crystals of PhBBO show metamagnetism, 11 which means that the antiparallel alignment of FM ordered p-stacks could be reversed by the field and, in turn, all chains would become aligned in a parallel fashion (see Fig. 8a and b for a schematic representation; notice that Fig. 8a shows a scheme of the ground spin state of Fig. 6). For a metamagnet to be detected, the magnetization as a function of the magnetic field M(H) must show an inflexion, with dM/dH(H) reaching a maximum value at the field at which the magnetic realignment occurs. Nevertheless, experimentally, no inflexion was observed in M(H), nor was any maximum in dM/dH as a function of H prior to saturation, even at T = 2 K. Thus, metamagnetism was ruled out. The available M(H), dM/dH(H) and ZFC-FC data established that AFM ordering collapses above H = 200 Oe. Accordingly, Oakley and coworkers indicated that this behavior was more consistent with a spin-flop transition (see Fig. 8a and c for a schematic representation).
In order to see which is the picture that works out better with the experimental data at the phase transition (4.5 K), we performed the optimization of the PhBBO crystal at zero Kelvin (0 K). Therefore, we will be able to discuss about the effect of the unit cell contraction upon cooling (from room temperature, RT, to 0 K). We performed variable-cell geometry optimization using periodic boundary conditions at the PBE-D2 level 22,24 (see computational details for further discussion). Notice that, although DFT is not capable of calculating the J AB value for certain PhBBO radical pairs, it is the most sensible choice to Fig. 8 (a) Schematic representation of the ground spin state obtained using the 296 K PhBBO crystal structure at the RASSCF-1 level. Schematic representation of metamagnetism: (a) to (b); and of a spin-flop transition: (a) to (c). (d) Schematic representation of the magnetic topology obtained using the 0 K PhBBO crystal structure at the RASSCF-1 level. The arrows represent the spin propagation according to the J AB magnetic interaction values to provide a schematic picture of the a and b spin distribution. Color code for the interactions: Dimer1 (dark blue, FM), Pdimer1 (orange, AFM), Pdimer2 (green, FM), and Pdimer6 (red, FM). undertake geometry optimizations, as already proved in the literature. 7b,25,28b,c The unit cell parameters slightly contract upon cooling from room temperature (X-ray data at 296 K for an orthorhombic space group P2 1 2 1 2 1 : a = 6.801(1) Å, b = 11.378(5) Å, c = 15.625(2) Å) to zero Kelvin (optimized data at 0 K: a = 6.610(1) Å, b = 11.267(5) Å, c = 15.421(1) Å) (see ESI, † Section S2). Despite the small unit cell contraction observed upon cooling, the effect on the J AB magnetic interactions is clear. Table 6 shows the J AB values that have been re-evaluated on the optimized nuclear coordinates by means of RASSCF-1/DZV. Notice that the J AB magnetic coupling for Dimer1 varies from +44.6 cm À1 at room temperature to +70.9 cm À1 at zero Kelvin. Dimer1 is thus another clear example that J AB can strongly depend on small changes in the interplanar distance between adjacent radicals within a p-stack. 28c Whereas the J AB magnetic interaction of Pdimer1 and Pdimer2 appears to be invariant upon cell contraction, the J AB of Pdimer6 is now comparable in absolute value to the magnitude of Pdimer1. It must be stressed that changes in the magnetic properties between radicals, here obtained upon cooling, can also be achieved by means of applying pressure, as has already been reported in the literature. 1c,39 Accordingly, the resulting magnetic topology at 0 K (Fig. 8d) differs significantly from the magnetic topology obtained at room temperature (Fig. 6). An analysis of the spin propagation according to the J AB values has been performed to provide a schematic and very simplistic picture of the a and b spin distribution at 0 K. Given a PhBBO radical in an A (bc-) layer, its coplanar neighboring radicals are oriented according to the AFM J(Pdimer1 opt ) and the FM J(Pdimer6 opt ) interactions. The first neighboring B (bc-) layer shows a large amount of geometrical spin frustration since it cannot satisfy either the AFM connection via Pdimer1 opt or the FM connection via Pdimer6 opt . In contrast, the second neighboring A (bc-) layer will again behave accordingly to AFM J(Pdimer1 opt ) and FM J(Pdimer6 opt ). Notice that an equivalent interpretation follows if a given PhBBO radical from a B (bc-) layer is taken as a reference. Regarding Pdimer2 opt (+0.1 cm À1 ) at 0 K, it introduces half of the geometrical spin frustration in comparison to that introduced at room temperature. It thus follows that, at the phase transition, the antiparallel-aligned 1D FM chains (Fig. 8a) convert into domains of weakly either FM-or AFM-coupled 1D FM chains (Fig. 8d).
Let us remind here that, experimentally, no inflexion was observed in M(H), nor was any maximum in dM/dH as a function of H prior to saturation. In fact, at 5 K and more clearly at 2 K, dM/dH(H) showed a minimum at zero applied field that rapidly led to a maximum value at AE10 kOe and AE20 kOe, respectively (see experimental M(H) and dM/dH(H) data in Yu et al. 11 ). Thus, in order to get a deeper insight into the phase transition that experimentally occurs at 4.5 K, M(H) simulations were performed (Fig. 9). The 2 Â 8 magnetic model was used together with both X-ray (296 K) and optimized (0 K) J AB values at the RASSCF-1 level. The magnetization M(H) simulated curves at 2, 5 and 10 K show in all cases a similar shape in comparison to experiment. Note that Table 6 Magnetic exchange couplings J AB (in cm À1 ) calculated for the seven pairs of radicals in the optimized 0 K crystal structure. The calculations were performed using RASSCF-1 active space (see main text for definition), and the DZV contraction of ANO-RCC basis set. The OÁ Á ÁO distance (in Å) between each pair of radicals is also indicated  simulated and experimental M(H) absolute values differ since we are using a finite magnetic model. However, it is dM/dH(H) that can enable us to univocally assign the true nature of the magnetism of the PhBBO purely organic crystal. Fig. 9a and b shows that the set of J AB magnetic couplings between PhBBO radicals calculated either from X-ray data or using the optimized crystal structure at 0 K cannot produce simulated dM/dH(H) data in agreement with experiment at 2 K and 5 K. The reason for this behavior is that they both share the same energy spectra and, in turn, produce a similar shape of M(H) and dM/dH(H) (see ESI, † Section S7, for details on the energy spectra). For the phase transition to occur, our study reveals that it is required that the AFM J inter-p-stack interactions have an absolute value within the range 10% to 5% of J p-stack , with J p-stack being either J X-ray,296 K (+44.6 cm À1 ) or J opt,0 K (+70.9 cm À1 ). According to our microscopic study, the J inter-p-stack estimated value would reasonably agree with the experimentally suggested mean field parameter (À2.5 cm À1 ). Indeed, we must stress that this cumulative inter-stack fitting value can be interpreted and assigned from first-principles. Thanks to these larger AFM J inter-p-stack interactions, the new energy spectrum has no degenerate lower lying energy states, i.e. it becomes an incommensurate system, which translates into a dM/dH(H) shape in agreement with experimental data (see Fig. 9c and d, and discussion in ESI, † Section S7). Apparently our calculations underestimate the inter-p-stack AFM coupling. The value of J inter-p-stack is rather small compared to the value of the intra-p-stack J AB interaction. Therefore, the quantitative evaluation of J inter-p-stack might require the inclusion of the electron correlation between the p and s subspaces (note that the RASSCF-1 space includes only the electron correlation of the p-space). In that case, the DDCI method would be more appropriate as it can consider a larger orbital space to provide highly accurate values for the singlet/triplet energy splitting. 37,38 The resulting J AB interactions might lead to a better agreement when simulating the dM/dH(H) maximum value. However, these calculations were not performed in this study since they would imply a remarkably large computational effort that will not provide new or different insights into the physical properties of the PhBBO system.

Conclusions
The p-stacking motif of crystals of purely organic PhBBO at room temperature entails a 3D magnetic topology of FM chains weakly knit together through both FM and AFM couplings, thus resulting in an overall antiparallel alignment of FM 1D p-stacked chains. Therefore, our results suggest that PhBBO cannot be merely viewed as a spin-canted AFM (SC-AFM) state as proposed experimentally, 11 since the inter-stack interactions are not only AFM in nature. Upon cooling from room temperature to zero Kelvin, our calculations uncover that there is an enhancement of p-stack and interp-stack FM interactions. At the phase transition (4.5 K), a suitable model tells us that there are domains of weakly either FM-or AFM-coupled 1D FM chains. As a result, there is a significant competition between FM chains aligning parallel (FM coupling) or antiparallel (AFM coupling) to each other. This inter-chain coupling competition clearly may lead to a certain degree of geometrical spin frustration upon phase transition at 4.5 K.
A comparison with available experimental magnetic susceptibility data has shown that, for simulation purposes, the interp-stack magnetic interactions must be taken into account to adequately calculate the magnetic susceptibility wT(T) and w(T) data. Besides, in order to reproduce the experimental magnetization M(H) and dM/dH(H) data, our study discloses that the inter-p-stack AFM interactions that couple the FM p-chains should be stronger than those computed in the present study. According to our results, the magnitude of J inter-p-stack should be (in absolute value) as large as 5-10% J p-stack magnetic coupling.
DFT dramatically fails to calculate the magnetic coupling between purely organic PhBBO radicals. Instead one has to resort to wavefunction-based methods, such as CASSCF and RASSCF. Unlike inorganic molecule-based magnets, 32 to our knowledge, this is the first time we encounter that the magnetic coupling between purely organic S = 1 2 radicals cannot be evaluated by means of DFT-based methods. An exhaustive computational study reveals that the reason behind the DFT failure turns out to be remarkably simple. Purely organic semiquinone-bridged bisdithiazolyl XBBO magnets are multireference systems due to the presence of low-lying open-shell states, unlike other bisdithiazolyl-based systems. 28,29c,30,34b This fact is indeed acknowledged in the literature as the multi-orbital effect. 1b,c It is then a must to use DFT with extreme care depending on the features of the studied system since it might be of no validity and lead to misconceptions. As a proof of concept, these limitations are illustrated for the phenyl-substituted semiquinone-bridged bisdithiazolyl PhBBO material. 11 We do believe that this conclusion transcends the specific material studied in this work, and is of general validity for the whole family of semiquinone-bridged bisdithiazolyl XBBO functional organic materials. 1a,7b,c,8-10