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

This journal is © The Royal Society of Chemistry 2021 J. Mater. Chem. C, 2021, 9, 10647–10660 | 10647 Cite this: J. Mater. Chem. C, 2021, 9, 10647 Pitfalls on evaluating pair exchange interactions for modelling molecule-based magnetism Maria Fumanal, a Joaquim Jornet-Somoza, abc Sergi Vela, a Juan J. Novoa, a Jordi Ribas-Arino * a and Merce `Deumal * a 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 wavefunctionbased 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. a Dept Cie `ncia de Materials i Quı ´mica Fı ´sica and IQTCUB, Facultat de Quı ´mica, Universitat de Barcelona, Martı ´i Franque `s 1, E-08028 Barcelona, Spain. E-mail: [email protected], [email protected] b IZO-SGI SGiker, Euskal Herriko Unibertsitatea (UPV/EHU), Joxe Mari Korta Center, Av. Tolosa 72, 20018 Donostia, Euskadi, Spain c The Max Planck Institute for the Structure and Dynamics of Matter (MPSD), Bldg. 99 (CFEL) Luruper Chaussee 149, 22761 Hamburg, Germany 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 p-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 E ´cole Polytechnique Fe ´de ´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. Cle ´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. Received 8th March 2021, Accepted 17th June 2021 DOI: 10.1039/d1tc01083b rsc.li/materials-c Journal of Materials Chemistry C PERSPECTIVE Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online View Journal | View Issue 10648 | J. Mater. Chem. C, 2021, 9, 10647–10660 This journal is © The Royal Society of Chemistry 2021 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. 1 There are two main strategies to reproduce experimental magnetic properties (e.g. magnetic susceptibility, heat capacity, magnetization) by means of statistical mechanics 2 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 3 ), 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. 4 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 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. Merce `Deumal. Then,hejoinedtheLaboratoirede 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. Cle ´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. Carbo ´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 Zu ¨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-Universita ¨t-Bochum with Prof. Dominik Marx. A Humboldt Fellowship supported part of his postdoc. In 2010, he returned to UB with a ‘‘Ramo ´ny 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. Merce `Deumal Dr Merce `Deumal is full Professor at Seccio ´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 Teo `rica i Computacional’ (awarded as Unit of Excellence Marı ´a de Maeztu), and ‘Xarxa de Refere `ncia en Quı ´mica Teo `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. Perspective Journal of Materials Chemistry C Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online This journal is © The Royal Society of Chemistry 2021 J. Mater. Chem. C, 2021, 9, 10647–10660 | 10649 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) 5 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 7 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 5 or wavefunction-based (e.g. CASSCF, RASSCF, DDCI) 8 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 shortrange, 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. 5 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). 9 The selection of the correct model must also consider the environment of the two radicals whose interaction has to be evaluated, namely counterions, 10 hydrogen bonding, 11 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 12 ) 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 p–plabile networks, 13 one might have to resort to a dynamic perspective, in which thermal vibrations are explicitly considered. 14 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=4in ref. 9bgenerates 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 15 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 Journal of Materials Chemistry C Perspective Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online 10650 | J. Mater. Chem. C, 2021, 9, 10647–10660 This journal is © The Royal Society of Chemistry 2021 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. 2 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 firstprinciples (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 7 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, 16 effect of the Madelung field, 10 effect of the crystal defects, 17 effect of the ligands, 11,18 magnetic dimensionality, 19 bistability 20 .... 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 2 (1,4diazacycloheptane) 2 Cl 4 (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 squaredplaquettes. 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 p-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. 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 7 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 26 are obtained from Molcas 7.6. 27 Otherwise energies at Fig. 1 Chemical formula for (a) Cu 2 (1,4-diazacycloheptane) 2 Cl 4 , CuHpCl, (b) phenylsemiquinone-bridged bisdithiazolyl, PhBBO, and (c) 1,3,5-trithia2,4,6-triazapentalenyl, TTTA. Perspective Journal of Materials Chemistry C Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online This journal is © The Royal Society of Chemistry 2021 J. Mater. Chem. C, 2021, 9, 10647–10660 | 10651 UB3LYP 28 or range-corrected CAM-B3LYP 29 level using 6-31+G(d) or 6-311++G(d,p) basis sets 30 are computed using Gaussian09 31 for the dimer model and Orca 3.0 32 package for the tetramer model. From the general Heisenberg Hamiltonian for a pair of S= 1/2 spin centres, H ˆ=2J AB S ˆ A S ˆ 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 singledeterminant broken-symmetry (BS) solution. 5b,c Within this approximation, the expression chosen to compute the energy difference is 6d ESET¼ 2ES BS ET  1þSab2(2) where 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)). 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 wT(T), heat capacity C p (T) and magnetization M(H). 2 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 33 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  w 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 33 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. 34 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 35 within the spin unrestricted formalism (broken symmetry singlet M S = 0 state) using plane wave pseudopotentials 36 expanded at a kinetic energy cutoff of 25 Ry, together with Vanderbilt ultrasoft pseudopotentials, 37 a fictitious mass for the orbitals of 400 a.u., and G-point sampling of the Brillouin zone. The semiempirical dispersion potential introduced by Grimme, 38 in its DFT-D2 parameterization is also considered. In addition, the AIMD simulations are performed in the canonical (or NVT) ensemble using Nose ´–Hoover chain thermostats. 39 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 vibrationallyaveraged magnetic susceptibility  w vib for the HT phase at 300 K is computed by averaging the wvalue 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 wvalue 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  w vib at 300 K 14a is last compared to the w value obtained from a static analysis 20 and from experiment. 24 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 2 (1,4-diazacycloheptane) 2 Cl 4 molecular transition metal AFM complex (in short CuHpCl, see Fig. 1a). Fig. 2 Tetramer model and set of equations used to calculate J(di) magnetic interactions between CuHpCl radicals (numbered 1to 4). HS/ LS stands for high spin and low spin, respectively. Journal of Materials Chemistry C Perspective Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online 10652 | J. Mater. Chem. C, 2021, 9, 10647–10660 This journal is © The Royal Society of Chemistry 2021 Since it was synthetized in 1997, 21 CuHpCl has been the prototypical compound to quote when a magnetically isolated spinladder 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. 40 It thus follows that this copper derivative is more challenging than anticipated, as it must have many competing microscopic magnetic 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. Perspective Journal of Materials Chemistry C Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online This journal is © The Royal Society of Chemistry 2021 J. Mater. Chem. C, 2021, 9, 10647–10660 | 10653 J AB interactions. The alleged spin ladder magnetic topology of CuHpCl has been a long lasting issue. 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 beascribedtoconsistofspin-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) 41 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.4bshowsaglobalview).Noticethattheremaining non-negligible J(di) interactions are one order of magnitude smaller than d1, d5, d12 and d14 (see Table 2). Therefore, these Table 1 Magnetic J(di) interactions (in cm 1 ) 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). 21 Journal of Materials Chemistry C Perspective Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online 10654 | J. Mater. Chem. C, 2021, 9, 10647–10660 This journal is © The Royal Society of Chemistry 2021 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 w(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 w(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 spinladderora3Dmodelofinteractingspin-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 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 22 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 p-stack motifs and paucity of close inter-column radicalradical contacts (see Fig. 5a). The experimental magnetic susceptibility wT(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 1 and a mean field parameter zJ0=2.5 cm 1 . PhBBO was also found to exhibit a phase transition at 4.5 K. In this case, DFT-based methods (such as unrestricted B3LYP 28 or range-corrected CAM-B3LYP 29 )provideagooddescription 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 1 at UB3LYP, and +193.3 cm 1 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 1 at CASSCF(14,14) and 49.7 cm 1 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 wTat any given temperature (red and green lines in Fig. 5c). Contrarily, at CASSCF(14,14) level, wTis 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. 23 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 Table 2 Magnetic J(di) interactions (in cm 1 ) 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 Fig. 5 (a) ab-View of the alternating ABAB p-stacks of PhBBO with intrastacks S3S1 contacts d 1 = 3.68 Å (in black) and d 2 = 3.81 Å (in red) along the a-axis (p-stacking axis). (b) ac-View of Dimer1 along the a-axis (p-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. 22 Perspective Journal of Materials Chemistry C Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online This journal is © The Royal Society of Chemistry 2021 J. Mater. Chem. C, 2021, 9, 10647–10660 | 10655 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. 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 p-stacks along a given crystallographic direction, as in the case of 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 1 )calculatedforDimer1using theroomtemperaturecrystalstructure.Thecalculationswereperformedwith 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 nonequivalent 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 Fig. 6 Packing of the monoclinic HT (SAXPOW05) and triclinic LT (SAXPOW06) crystal structures recorded at room temperature (RT): pstacks 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 p-stacking pair. (c) Magnetic susceptibility w(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) w(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). 24 Refcodes SAXPOW01-06 given by the Cambridge Crystallographic Data Base CCDC. 42 Journal of Materials Chemistry C Perspective Open Access Article. Published on 05 July 2021. Downloaded on 1/17/2022 1:18:43 PM. This article is licensed under a Creative Commons Attribution-NonCommercial 3.0 Unported Licence. View Article Online