# Magnetic exchange coupling in Cu dimers studied with modern multireference methods and broken-symmetry coupled cluster theory

Gurjot Singh, Stefani Gamboa, Maylis Orio, Dimitrios Pantazis, Michael Roemelt

► To cite this version:

Gurjot Singh, Stefani Gamboa, Maylis Orio, Dimitrios Pantazis, Michael Roemelt. Magnetic exchange coupling in Cu dimers studied with modern multireference methods and broken-symmetry coupled cluster theory. Theoretical Chemistry Accounts: Theory, Computation, and Modeling, 2021, 140 (10), pp.139. 10.1007/s00214-021-02830-0. hal-03862126

HAL Id: hal-03862126

https://hal.science/hal-03862126

Submitted on 22 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Magnetic exchange coupling in Cu dimers studied with modern multireference methods and broken-symmetry coupled cluster theory

Gurjot Singh<fnref n="1" /> · Stefani Gamboa<fnref n="2" /> · Maylis Orio<fnref n="2" /> · Dimitrios A. Pantazis<fnref n="3" /> · Michael Roemelt<fnref n="1" />

Received: 2 July 2021 / Accepted: 19 August 2021 / Published online: 20 September 2021
© Springer-Verlag GmbH Germany, part of Springer Nature 2021

## Abstract

Spin-state energetics of exchange-coupled copper complexes pose a persistent challenge for applied quantum chemistry. Here, we provide a comprehensive comparison of all available theoretical approaches to the problem of exchange coupling in two antiferromagnetically coupled bis-μ-hydroxo Cu(II) dimers. The evaluated methods include multireference methods based on the density matrix renormalization group (DMRG), multireference methods that incorporate dynamic electron correlation either perturbatively, such as the N-electron valence state perturbation theory, or variationally, such as the difference-dedicated configuration interaction. In addition, we contrast the multireference results with those obtained using broken-symmetry approaches that utilize either density functional theory or, as demonstrated here for the first time in such systems, a local implementation of coupled cluster theory. The results show that the spin-state energetics of these copper dimers are dominated by dynamic electron correlation and represent an impossible challenge for multireference methods that rely on brute-force expansion of the active space to recover correlation energy. Therefore, DMRG-based methods even at the limit of their applicability cannot describe quantitatively the antiferromagnetic exchange coupling in these dimers, in contrast to dinuclear complexes of earlier transition metal ions. The convergence of the broken-symmetry coupled cluster approach is studied and shown to be a limiting factor for the practical application of the method. The advantages and disadvantages of all approaches are discussed, and recommendations are made for future developments.

**Keywords** Exchange coupling · Copper · Multireference methods · Broken symmetry · Coupled cluster

## 1 Introduction

Local spins in inorganic and bioinorganic systems with multiple open-shell transition metals couple to a multitude of total spin states whose energy is a function of the total spin S. This phenomenon is called magnetic exchange coupling

and is usually described by the Heisenberg–Dirac–van Vleck Hamiltonian [1–4]

$$ \hat{H}_{HDVV} = -2 \sum_{I<j} J_{II} \hat{S}_I \hat{S}_J \quad (1) $$

$\hat{H}_{HDVV}$ is an effective Hamiltonian that acts only on a fictitious set of spin states without any reference to the spatial wavefunction. Instead, the coupling of each pair of spins $\hat{S}_I$ and $\hat{S}_J$ is associated with a coupling parameter $J_{II}$. Negative values for $J_{II}$ correspond to antiferromagnetic coupling, i.e., the lowest total spin state is energetically favored. In contrast, positive values describe a situation in which the highest total spin state is energetically favored, hence ferromagnetic coupling. While magnetic exchange coupling has notable effects on spectroscopic and chemical properties of polynuclear transition metal compounds, the accurate prediction of coupling constants by electronic structure calculations remains a great challenge for quantum chemistry.

The problem of magnetic coupling in polynuclear transition metal clusters is intimately related to the

<footnote>Maylis Orio
maylis.orio@univ-amu.fr</footnote>

<footnote>Dimitrios A. Pantazis
dimitrios.pantazis@kofo.mpg.de</footnote>

<footnote>Michael Roemelt
michael.roemelt@theochem.rub.de</footnote>

<footnote n="1">¹ Lehrstuhl Für Theoretische Chemie, Ruhr-Universität Bochum, 44780 Bochum, Germany</footnote>

<footnote n="2">² Institut des Sciences Moléculaires de Marseille, Aix Marseille Université, CNRS, Centrale Marseille, ISM2 UMR 7313, 13397 Marseille, France</footnote>

<footnote n="3">³ Max-Planck Institut Für Kohlenforschung, Kaiser-Wilhelm Platz 1, 45470 Mülheim an der Ruhr, Germany</footnote>
multideterminantal character of the involved electronic states and the concomitant static electron correlation. For example, a binuclear copper(II) system with local spins of $S = 1/2$ is characterized by two total spin states arising from the same ground-state electron configuration, and only one of them can be represented by a single determinant. Therefore, a straightforward approach to the theoretical description of this phenomenon must properly account for this aspect of the electronic structure, as provided by multireference (MR) methods.

The most commonly used methods of this family of electronic structure methods are the complete active space self-consistent field (CASSCF) and the complete active space configuration interaction (CASCI). A key feature of both methods is the division of the molecular orbital space into a set of internal orbitals that are always doubly occupied, a set of virtual orbitals that are always empty and the set of active orbitals that can have any occupation between 0 and 2. More precisely, the CASCI and CASSCF wavefunctions are linear combinations of all possible active space configuration state functions, thus solving the full configuration interaction (FCI) problem in this subset of orbitals. Although CASCI and CASSCF provide the required flexibility to yield a qualitatively correct wavefunction (given the active space is chosen to be adequately large), a quantitatively correct description of magnetic coupling can only be achieved when dynamic electron correlation is taken into account [5–7]. Well-established approaches to accomplish this crucial task are based on perturbation theory, configuration interaction or coupled cluster theory [8–11].

A common concern for all multireference methods is their high computational cost, arising from the underlying CASSCF calculation and the treatment of dynamic electron correlation. As the active space size increases, both aspects of an MR calculation quickly become unfeasible. Although modern implementations allow for CASSCF calculations with up to ca. 18 active electrons and orbitals (up to 22 electrons and orbitals have even been reached with massively parallel implementations), the treatment of dynamic electron correlation causes even more severe restrictions on the active space size. Therefore, a number of previous studies of magnetic couplings relied on a minimal active space that comprises only the partially occupied metal d-orbitals [5, 12–14]. Since this active space does not capture all of the essential physics of magnetic exchange coupling, sophisticated and computationally demanding approaches to dynamic electron correlation such as the difference-dedicated CI (DDCI) are needed to obtain qualitatively and quantitatively accurate results [15, 16]. Consequently, this approach is limited to small systems with few unpaired electrons.

In recent years, a different strategy to compute magnetic exchange couplings based on large active spaces has been explored. This conceptual shift was triggered by the

development of techniques that allow for large active spaces in CASSCF-type calculations on the order of 20–50 orbitals at considerably reduced computational cost. Among these techniques, the density matrix renormalization group (DMRG) is the most prominent [17–19]. Relying on a so-called matrix-product state wavefunction, the DMRG optimizes only a polynomial number of parameters to solve a problem whose complexity grows exponentially with increasing system size. A number of pilot studies of Fe, Cr and Mn clusters showcased the feasibility of DMRG calculations for the prediction of magnetic exchange coupling constants [20–23]. However, previous work by some of us revealed that for an oxo-bridged mixed-valence Mn dimer truly quantitative accuracy can only be reached when dynamic electron correlation is treated on top of a relatively large active space that comprises all Mn d-orbitals and all occupied valence orbitals of the oxo-bridges [6].

On the opposite end of computational complexity, a widely employed approach to the description of magnetic exchange coupling involves the broken-symmetry (BS) approach, particularly in combination with unrestricted density functional theory (DFT) as introduced by Noodleman [24, 25]. Instead of a multideterminantal wavefunction, the BS approach relies on a single determinant to indirectly deduce the energy of the antiferromagnetically coupled state. For example, in the aforementioned Cu(II) dimer one Cu center would be associated with positive spin density, while the other Cu center would feature negative spin density. As a consequence, the BS determinant does not feature the correct spin symmetry. Nevertheless, projection schemes allow for an estimation of magnetic exchange coupling constants based on the energies of the high-spin determinant and the BS determinant as well as the corresponding expectation values for the total spin $\langle S^2 \rangle$. Compared to the multideterminantal approaches discussed above, BS-DFT yields magnetic exchange coupling constants at the cost of two self-consistent field (SCF) calculations, thus making it feasible even for large systems with more than 100 atoms. Yet the BS-DFT approach is problematic because the properties of spin states other than the highest-spin state cannot be obtained directly, and because in practice the usage of different functionals and spin projection schemes leads to different exchange coupling values in an often unsystematic way. The dependence of BS-DFT results on the employed functional originates from the different treatment of dynamic electron correlation which, owing to the unknown form of the exact Hohenberg–Kohn functional, cannot be systematically improved [26]. Recently, it was suggested that instead one may employ the same broken-symmetry approach but using the well-established coupled cluster theory for the treatment of dynamic electron correlation [27]. This is a potentially powerful approach that remains to be explored. In the present work, we assess the performance of all four
<figure><img src="image_2.png" /><figcaption>Fig. 1 The two dinuclear CuII complexes studied in the present work</figcaption></figure>

existing approaches to the exchange coupling problem, using the two archetypical antiferromagnetically coupled bis-μ-hydroxo-bridged copper dimers **1** and **2** depicted in Fig. 1.

**2 Computational details**

This work reports on the results of several sets of electronic structure calculations: geometry optimizations on the DFT level of theory and single point energy evaluations using DFT, DMRG-CI, [28] DMRG-SCF, [29, 30] strongly contracted DMRG-NEVPT2, [31–33] DDCI, [15, 16] and broken-symmetry coupled cluster theory [27]. The reported DMRG and DMRG-NEVPT2 calculations were conducted with the MOLBLOCK program, while all other calculations were performed with the ORCA program package [6, 34, 35].

To facilitate comparisons between theory and experiment, we optimized the X-ray structure of complexes **1** (ID: HTMCUP) and **2** (ID: CUXJEK) while constraining the positions of all heavy atoms to their experimentally derived coordinates [36, 37]. Thus, only the positions of hydrogen atoms were relaxed since these are not reliably determined from X-ray diffraction experiments. Geometry optimizations were conducted using the TPSS functional with the def2-TZVP basis set and taking advantage of the RI approximation with the appropriate auxiliary Coulomb fitting basis sets [38–42]. Subsequent DFT single point calculations were carried out with the same basis set and a large variety of functionals (see below). Increased integration grids and tight SCF convergence criteria were used throughout. It should be noted that the usage of Grimme’s dispersion correction with Becke–Johnson damping (D3BJ) during the geometry optimization did not result in a significant change of the obtained results (see Supporting Information).

The reported high-spin and broken-symmetry coupled cluster (CC) calculations incorporated single and double excitations (CCSD). Owing to the size of **1** and **2**, the solution of the canonical CCSD equations is intractable. To reduce the associated computational cost, the corresponding equations were solved in their local pair-natural orbital (LPNO-CCSD) form [43–46]. All reported coupled cluster

calculations relied on unrestricted Kohn–Sham reference determinants that were generated with the B3LYP functional. [47, 48] If not stated otherwise, the def2-TZVP basis set was used in CC calculations.

DDCI calculations incorporated excitations with up to three degrees of freedom (DDCI3) employed the def2-SVP basis set as the maximum applicable basis for the present systems [15, 16, 39]. Selection thresholds have an important effect on the predicted spin-state energetics, and these will be discussed in the appropriate section where the results are presented.

DMRG-SCF, DMRG-CI and DMRG-NEVPT2 calculations were performed using the MOLBLOCK in-house code [6, 34]. In the current context, MOLBLOCK employs the BLOCK DMRG code developed by Chan and co-workers as approximate Full-CI solver [17, 28, 49]. Accordingly, all reported DMRG-NEVPT2 calculations use the implementation reported in reference [32]. In the following, active spaces are labeled (n,m) to denote n electrons in m orbitals. Starting orbitals for CASSCF calculations with small active spaces ((2,2) and (14,8)) were obtained from DFT calculations using the B3LYP functional [47, 48]. During these DFT calculations, the RIJCOSX approximation with the def2/J auxiliary basis was applied [50–52]. Multireference calculations with larger active space are based on the CASSCF(14,8) orbitals as described in detail below. Wavefunction calculations used the RI approximation where applicable with the def2-TZVP/C auxiliary basis set (def2-SVP/C for DDCI calculations) [53]. The state-averaged approach was followed for orbital optimization in CASSCF and DMRG-SCF calculations.

**3 Results and discussion**

**3.1 A. Description of the Cu dimers**

The structures of **1** and **2** display the same coordination core $[Cu_2(\mu-OH)_2(L)_2]^{2+}$ with two Cu(II) bridging two hydroxo ligands (Fig. 1 and Figure S1). Complex **1** features a square planar geometry with each copper center coordinating two nitrogen atoms of the N,N,N',N''-tetramethylethylenediamine ligand [36]. Complex **2** exhibits a square pyramidal geometry with the metal being now coordinated to three nitrogen atoms of the N,N',N''-trimethyl-1,4,7-triazacyclononane ligand [37]. The Cu–O-Cu angles are 101.6° and 100.1° for complexes **1** and **2**, respectively (see Table S1 for detailed structural parameters). Coordination distances around the Cu(II) centers are in the expected range: For complex **1**, Cu–O distances are 1.931 Å and 1.897 Å, and Cu–N distances are 1.996 Å and 2.033 Å, while for complex **2** Cu–O distances are 1.938 Å and 1.935 Å, and equatorial Cu–N distances are 2.084 Å and 2.060 Å, with an elongated
apical Cu–N distance of 2.238 Å. The experimentally fit-
ted exchange coupling constants *J* for both complexes are
antiferromagnetic, *J*(1) = −180 cm⁻¹ and *J*(2) = −45 cm⁻¹.

**3.2 B. Broken-symmetry density functional theory**

The accuracy of BS-DFT for the calculation of *J* has notably
increased with the use of modern hybrid functionals, and
nowadays, its applicability expands to diverse and complex
systems. However, many examples show that the optimal
functional is not transferable from system to system. For
example, it was observed that B3LYP [47, 48] can give
satisfactory results for copper [54–57] and cobalt [58, 59]
complexes, but is not the best functional for chromium [60],
manganese [61] and iron [62] dimers. In these latter cases,
TPSS0 [63, 64], TPSSh [65] and M06 [66], respectively, are
the best performing functionals to describe the experimental
results. There have also been intense debates [67–69] regard-
ing the correct way of calculating the exchange coupling
constant from the high-spin and BS energies, since three dif-
ferent schemes can be employed within the BS-DFT frame-
work (eqs. (2)–(5)).[70] These formalisms were described
in works from Noodleman's (Eqs. (2) and (5)) [24], Ruiz's
(Eq. (3)) [71] and Yamaguchi's (Eq. (4)) [72] groups:

$$
J_1 = - \frac{E_{HS} - E_{BS}}{S_{max}^2} \qquad (2)
$$

$$
J_2 = - \frac{E_{HS} - E_{BS}}{S_{max}(S_{max} + 1)} \qquad (3)
$$

$$
J_3 = - \frac{E_{HS} - E_{BS}}{\langle S^2 \rangle_{HS} - \langle S^2 \rangle_{BS}} \qquad (4)
$$

$$
J_4 = \frac{E_{HS} - E_{BS}}{(1 + S_{ab})^2} \qquad (5)
$$

where E<sub>HS</sub> and E<sub>BS</sub> are the high-spin and broken-symmetry
energies, respectively, S<sub>ab</sub> is the orbital overlap, S<sub>max</sub> is the
total spin of the high-spin state, and ⟨S<sup>2</sup>⟩<sub>HS</sub> and ⟨S<sup>2</sup>⟩<sub>BS</sub> are
the spin expectation values for the high-spin and broken-
symmetry states, respectively.

Noodleman's method (eq. (2)), also referred to as the
spin-projected approach, computes the exchange coupling
constant *J* in the weak limit of overlap of the magnetic orbit-
als [26]. Alternatively, Ruiz and co-workers developed the
non-projected approach (eq. (3)) assuming a strong coupling
limit, arguing that the broken-symmetry energy is a good
approximation of the low-spin state of the system [67, 71].
Yamaguchi proposed a more general spin-projected approach
(eq. (4)) that adapts and covers all situations, from weak-to-
strong coupling [72]. Similarly to the latter, Noodleman also

proposed a formulation based on the overlap of the magnetic
orbitals accounting for all coupling regimes (Eq. 5).

Discrepancies have been observed in the use of the equa-
tions to calculate *J* values without a clear justification of
the choice. Regarding dinuclear copper and heterobimetallic
copper complexes, Santiago et al. [73] favored computation
of exchange coupling constants using eq. (2). On the con-
trary, Luo et al. [74] and Ruiz et al. [75, 76] described dinu-
clear copper systems with eq. (3). Finally, Luo et al. [57],
Reis et al.[77] and Comba [55] used eq. (4) in their studies.
Moreover, Simões et al. [78] observed a better quantitative
agreement between theoretical and experimental results
when using eqs. (2) and (4), while eq. (3) highly under-
estimated the values of the calculated exchange coupling
constants. Systematic inferences are not possible because
different spin projection approaches and different choices of
functional can have effects of similar magnitude and interact
in ways that can be unpredictable.

To assess the performance of BS-DFT for predicting the
exchange coupling constant *J* of the selected copper com-
plexes **1** and **2**, we present a series of calculations using
the generalized-gradient-approximation (GGA) functionals
BLYP, [47, 79] PBE, [80, 81] and BP86 [82], the meta-
GGA TPSS, [38] hybrid GGA B3LYP, [47, 48] B1LYP, [48]
PBE0 [83] and B3PW91 [48, 80] and meta-hybrid GGA
functionals TPSSh [65] and TPSS0 [63, 64]. Tables 1 and 2
display the computed *J*<sub>1</sub>, *J*<sub>2</sub>, *J*<sub>3</sub> and *J*<sub>4</sub> for complex **1** and **2**,
respectively, along with the absolute percentage deviation
(APD) calculated with respect to the experimental results.
Using any of Eqs. (2–5), we predict negative *J* values for the
two complexes, in line with the experimentally established
antiferromagnetic coupling between the two Cu(II) centers,
leading to a spin singlet (*S*=0) ground state.

Prior to a detailed discussion of the numerical perfor-
mance of the various tested functionals in the current con-
text, it is instructive to examine the nature of antiferromag-
netic coupling between the metal centers. To this end, the
corresponding orbitals presented in Fig. 2 (here, the TPSSh
orbitals are shown as representative case) were constructed
according to established procedures and their shape ana-
lyzed. For both complexes, the magnetic orbitals are the 3
d<sub>x<sup>2</sup>-y<sup>2</sup></sub> Cu orbitals in σ* combination with ligand p orbitals.
Exchange is mediated mostly by overlap over the bridging
hydroxyl groups (2p<sub>x</sub> and 2p<sub>y</sub> orbitals of the bridging O
atoms). According to the detailed results listed in Tables 1
and 2, each individual functional provides greater overlap
integrals for complex **1** than for complex **2**, in line with the
stronger antiferromagnetic coupling describing complex **1**
compared to complex **2**. The magnitude of the overlap inte-
gral is correlated with the nature of the functional, specifi-
cally with the percentage of exact exchange. This in turn is
correlated with the energetic separation between the BS and
<caption>**Table 1** BS-DFT results for complex 1. DFT-calculated spatial overlap ($S_{ab}$) of corresponding orbitals in the broken-symmetry state, exchange coupling constants (J, cm⁻¹) and absolute percentage deviation (APD) for the 4 different spin projection formalisms. The experimental J value is –180 cm⁻¹</caption>

<table><caption>**Table 2** BS-DFT results for complex 2. DFT-calculated spatial overlap ($S_{ab}$) of corresponding orbitals in the broken-symmetry state, exchange coupling constants (J, cm⁻¹) and absolute percentage deviation (APD) for the 4 different spin projection formalisms. The experimental J value is –45 cm⁻¹</caption><thead><tr><th>Functionals</th><th>% of HF exchange</th><th>S<sub>ab</sub></th><th>J<sub>1</sub></th><th>J<sub>2</sub></th><th>J<sub>3</sub></th><th>J<sub>4</sub></th><th>APD J<sub>1</sub></th><th>APD J<sub>2</sub></th><th>APD J<sub>3</sub></th><th>APD J<sub>4</sub></th></tr></thead><tbody><tr><td>TPSSh</td><td>10</td><td>0.219</td><td>-336</td><td>-168</td><td>-321</td><td>-336</td><td>87</td><td>7</td><td>78</td><td>87</td></tr><tr><td>B3LYP</td><td>20</td><td>0.164</td><td>-221</td><td>-111</td><td>-215</td><td>-221</td><td>23</td><td>38</td><td>19</td><td>23</td></tr><tr><td>B3PW91</td><td>20</td><td>0.163</td><td>-220</td><td>-110</td><td>-214</td><td>-220</td><td>22</td><td>39</td><td>19</td><td>22</td></tr><tr><td>B1LYP</td><td>25</td><td>0.136</td><td>-167</td><td>-84</td><td>-164</td><td>-167</td><td>7</td><td>53</td><td>9</td><td>7</td></tr><tr><td>PBE0</td><td>25</td><td>0.134</td><td>-166</td><td>-83</td><td>-163</td><td>-166</td><td>8</td><td>54</td><td>9</td><td>8</td></tr><tr><td>TPSS0</td><td>25</td><td>0.124</td><td>-150</td><td>-75</td><td>-148</td><td>-150</td><td>17</td><td>58</td><td>18</td><td>17</td></tr></tbody></table><table><thead><tr><th>Functionals</th><th>% of HF exchange</th><th>S<sub>ab</sub></th><th>J<sub>1</sub></th><th>J<sub>2</sub></th><th>J<sub>3</sub></th><th>J<sub>4</sub></th><th>APD J<sub>1</sub></th><th>APD J<sub>2</sub></th><th>APD J<sub>3</sub></th><th>APD J<sub>4</sub></th></tr></thead><tbody><tr><td>BLYP</td><td>0</td><td>0.441</td><td>-740</td><td>-370</td><td>-619</td><td>-740</td><td>311</td><td>106</td><td>244</td><td>311</td></tr><tr><td>PBE</td><td>0</td><td>0.436</td><td>-725</td><td>-363</td><td>-609</td><td>-725</td><td>303</td><td>102</td><td>238</td><td>303</td></tr><tr><td>BP86</td><td>0</td><td>0.433</td><td>-723</td><td>-361</td><td>-608</td><td>-723</td><td>302</td><td>101</td><td>238</td><td>302</td></tr><tr><td>TPSS</td><td>0</td><td>0.371</td><td>-627</td><td>-314</td><td>-551</td><td>-627</td><td>248</td><td>74</td><td>206</td><td>248</td></tr><tr><td>TPSSh</td><td>10</td><td>0.219</td><td>-336</td><td>-168</td><td>-321</td><td>-336</td><td>87</td><td>7</td><td>78</td><td>87</td></tr><tr><td>B3LYP</td><td>20</td><td>0.164</td><td>-221</td><td>-111</td><td>-215</td><td>-221</td><td>23</td><td>38</td><td>19</td><td>23</td></tr><tr><td>B3PW91</td><td>20</td><td>0.163</td><td>-220</td><td>-110</td><td>-214</td><td>-220</td><td>22</td><td>39</td><td>19</td><td>22</td></tr><tr><td>B1LYP</td><td>25</td><td>0.136</td><td>-167</td><td>-84</td><td>-164</td><td>-167</td><td>7</td><td>53</td><td>9</td><td>7</td></tr><tr><td>PBE0</td><td>25</td><td>0.134</td><td>-166</td><td>-83</td><td>-163</td><td>-166</td><td>8</td><td>54</td><td>9</td><td>8</td></tr><tr><td>TPSS0</td><td>25</td><td>0.124</td><td>-150</td><td>-75</td><td>-148</td><td>-150</td><td>17</td><td>58</td><td>18</td><td>17</td></tr></tbody></table>

<table><thead><tr><th>Functionals</th><th>% of HF exchange</th><th>S<sub>ab</sub></th><th>J<sub>1</sub></th><th>J<sub>2</sub></th><th>J<sub>3</sub></th><th>J<sub>4</sub></th><th>APD J<sub>1</sub></th><th>APD J<sub>2</sub></th><th>APD J<sub>3</sub></th><th>APD J<sub>4</sub></th></tr></thead><tbody><tr><td>TPSSh</td><td>10</td><td>0.148</td><td>-93</td><td>-47</td><td>-91</td><td>-93</td><td>107</td><td>4</td><td>102</td><td>107</td></tr><tr><td>B3LYP</td><td>20</td><td>0.111</td><td>-39</td><td>-19</td><td>-38</td><td>-39</td><td>13</td><td>58</td><td>16</td><td>13</td></tr><tr><td>B3PW91</td><td>20</td><td>0.109</td><td>-36</td><td>-18</td><td>-36</td><td>-36</td><td>20</td><td>60</td><td>20</td><td>20</td></tr><tr><td>B1LYP</td><td>25</td><td>0.091</td><td>-18</td><td>-9</td><td>-18</td><td>-18</td><td>67</td><td>80</td><td>60</td><td>60</td></tr><tr><td>PBE0</td><td>25</td><td>0.089</td><td>-16</td><td>-8</td><td>-16</td><td>-16</td><td>64</td><td>82</td><td>64</td><td>64</td></tr><tr><td>TPSS0</td><td>25</td><td>0.082</td><td>-15</td><td>-7</td><td>-15</td><td>-15</td><td>67</td><td>84</td><td>67</td><td>67</td></tr></tbody></table><table><thead><tr><th>Functionals</th><th>% of HF exchange</th><th>S<sub>ab</sub></th><th>J<sub>1</sub></th><th>J<sub>2</sub></th><th>J<sub>3</sub></th><th>J<sub>4</sub></th><th>APD J<sub>1</sub></th><th>APD J<sub>2</sub></th><th>APD J<sub>3</sub></th><th>APD J<sub>4</sub></th></tr></thead><tbody><tr><td>BLYP</td><td>0</td><td>0.301</td><td>-261</td><td>-131</td><td>-240</td><td>-261</td><td>480</td><td>191</td><td>433</td><td>480</td></tr><tr><td>PBE</td><td>0</td><td>0.294</td><td>-247</td><td>-123</td><td>-227</td><td>-247</td><td>449</td><td>173</td><td>404</td><td>448</td></tr><tr><td>BP86</td><td>0</td><td>0.293</td><td>-246</td><td>-123</td><td>-227</td><td>-246</td><td>447</td><td>173</td><td>404</td><td>448</td></tr><tr><td>TPSS</td><td>0</td><td>0.252</td><td>-219</td><td>-110</td><td>-206</td><td>-219</td><td>387</td><td>144</td><td>358</td><td>387</td></tr><tr><td>TPSSh</td><td>10</td><td>0.148</td><td>-93</td><td>-47</td><td>-91</td><td>-93</td><td>107</td><td>4</td><td>102</td><td>107</td></tr><tr><td>B3LYP</td><td>20</td><td>0.111</td><td>-39</td><td>-19</td><td>-38</td><td>-39</td><td>13</td><td>58</td><td>16</td><td>13</td></tr><tr><td>B3PW91</td><td>20</td><td>0.109</td><td>-36</td><td>-18</td><td>-36</td><td>-36</td><td>20</td><td>60</td><td>20</td><td>20</td></tr><tr><td>B1LYP</td><td>25</td><td>0.091</td><td>-18</td><td>-9</td><td>-18</td><td>-18</td><td>67</td><td>80</td><td>60</td><td>60</td></tr><tr><td>PBE0</td><td>25</td><td>0.089</td><td>-16</td><td>-8</td><td>-16</td><td>-16</td><td>64</td><td>82</td><td>64</td><td>64</td></tr><tr><td>TPSS0</td><td>25</td><td>0.082</td><td>-15</td><td>-7</td><td>-15</td><td>-15</td><td>67</td><td>84</td><td>67</td><td>67</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Corresponding orbitals and spatial overlap integrals $S_{ab}$ from TPSSh calculations describing the broken-symmetry states of complex 1 (top) and complex 2 (bottom)</figcaption></figure>

high-spin states, which affects the derived exchange coupling constant.

As expected from past experience [84], our calculations follow the trend that hybrid functionals are more accurate in predicting *J* values than GGA or meta-GGA functionals, which largely overestimate the magnitude of the coupling. It is noted that the choice of equation for the *J* is crucial in determining the final value. Although B1LYP and PBE0 perform well for complex 1 when considering *J*₁ and *J*₃, the results using these functionals appear disappointing for complex 2. A combination that appears to work well for both complexes is the TPSSh functional with the *J*₂ formalism. This is associated with the lowest APD overall (see also Figure S2 for a graphical representation of all APDs). These results highlight the strong dependence of the BS-DFT approach on both, the functional and the formalism used to calculate the exchange couplings. Therefore, DFT can provide acceptable agreement with experiment at an empirical level, i.e., as long as an adequate screening of functionals is performed on the studied systems. Clearly, however, the wide range of results that can be obtained means that the predictive ability of the approach is limited, at least for the present systems, and that it is hard to identify or quantify the role of error cancellation. To overcome the fundamental weakness of the BS-DFT approach, more advanced electronic structure methods need to be employed and the applicability of such methods will be considered in the following.

### 3.3 C. CASSCF- and DMRG-based approaches

In contrast to Kohn–Sham DFT, ab initio multireference methods do not rely on any parameters or procedures of semi-empirical nature but constitute a systematically improvable approach to the calculation of coupling constants through direct calculation of spin-state energies. Besides the basis set size, the size of the *N*-electron Hilbert space in which the wavefunction is expanded is the only factor that determines the accuracy of a given calculation. As discussed above, one way to approach the required numerical accuracy for the present problem is to utilize the DMRG to cover static correlation effects within a large set of active orbitals, while all remaining dynamic electron correlation effects are treated by means of second-order perturbation theory, i.e., NEVPT2 [6, 22, 23]. Importantly, this strategy allows one to analyze the nature of antiferromagnetic coupling through a gradual extension of the active orbital space in various ways. By permuting combinations of groups of chemically equivalent (or similar) orbitals within reasonable perimeters, the most important pathway for the observed magnetic interaction can be identified [6, 7, 60]. Alternatively, one may analyze the entanglement within the active space to gain information about the coupling pathway [85, 86]. It should be noted at this point that in many cases where the experimental exchange coupling constant is unknown, it quickly becomes impractical to increase the active space size or the level of theory used to describe dynamic electron correlation until convergence of the computed exchange couplings is achieved.
<figure><img src="image_4.png" /><figcaption>**Fig. 3** Definition of orbital groups for complex 1, used in the construction of different active spaces for multireference calculations. Groups A to D correspond to sets of ligand-based and copper-based occupied orbitals</figcaption></figure>

In the present work, we have defined 7 groups of orbitals labeled A through G, shown in Figs. 3 and 4 to compose different active spaces employed in a series of multireference calculations. Group A consists of 4 MOs that are dominated by in-plane N-2p contributions, thus contributing to the Cu–N bonding, while group B comprises 6 MOs of mostly O-2p character. The MOs in Group C are essentially linear combinations of the “doubly occupied” Cu 3d orbitals (xy, xz, yz and z²), whereas group D includes linear combinations of the formally “singly occupied” Cu 3 dx²-y² orbitals. Finally, groups E, F and G (Fig. 4) correspond to “second-shell” counterparts of groups B, C and D in the sense that they have mostly Cu 4d and O 3p character.

As pointed out previously, optimizing the orbitals according to the variational principle is imperative to obtain accurate results. However, some orbitals that are expected to play

a prominent role in facilitating magnetic couplings feature occupation numbers close to 0 or 2 rendering the orbital optimization procedure with these orbitals in the active space difficult if not impossible. Therefore, we have omitted orbitals from the active space during orbital optimization if they were found to be close to doubly occupied or unoccupied after a considerable number of optimization cycles. In particular, orbital groups E, F and G could not be fully optimized regardless of the active space composition. Instead they were localized by means of the Pipek–Mezey localization scheme in the spirit of the “split localization” approach previously used in large-scale multireference calculations [28]. In the discussion below, we will refer to the combination of orbital optimization and virtual orbital localization as DMRG-SCF/CI approach. To avoid numerical biases toward any spin state, all orbital optimizations were conducted in
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Definition of orbital groups for complex **1**, used in the construction of different active spaces for multireference calculations. Groups E to G correspond to sets of unoccupied orbitals with Cu 4d and O 3p character</figcaption></figure>

a state-averaged fashion with respect to the first singlet and triplet state [6].

Tables **3** and **4** present the calculated magnetic coupling constants *J* for various active spaces using the bare CASSCF or DMRG-SCF/CI methods and their respective NEVPT2 implementations. The composition of the corresponding active spaces in terms of orbital groups A through G is given in the respective second column. In general, the exchange coupling constants evaluated based on the NEVPT2 relative spin-state energies are more negative than the ones obtained

<table><caption>**Table 3** Exchange coupling constants *J* (cm⁻¹) from CASSCF and CASSCF-NEVPT2 calculations for complexes **1** and **2**</caption><thead><tr><th rowspan="2">Active Space</th><th rowspan="2">Orbital Composition</th><th colspan="2">J(1)</th><th colspan="2">J(2)</th></tr><tr><th>CASSCF</th><th>NEVPT2</th><th>CASSCF</th><th>NEVPT2</th></tr></thead><tbody><tr><td>(2,2)</td><td>D</td><td>-4.1</td><td>-23.2</td><td>7.55</td><td>6.50</td></tr><tr><td>(14,8)</td><td>B, D</td><td>-8.0</td><td>-42.4</td><td>4.55</td><td>0.25</td></tr><tr><td>(18,10)</td><td>C, D</td><td>-4.3</td><td>-26.2</td><td>7.60</td><td>-1.20</td></tr></tbody></table>

from CASSCF or DMRG-SCF/CI energies. This reflects the critical role of dynamic electron correlation for the description of superexchange which leads to antiferromagnetic coupling in **1** and **2**. In the following, the results are discussed in view of the contribution from different orbital sets to the predicted exchange coupling constants.

As expected, the coupling constants obtained with a minimal active space that includes only the two singly occupied orbitals exhibit considerable deviations from the experimental value. For **1**, the CASSCF(2,2) energies result in $J(1) = -4.1$ cm⁻¹. While this value correctly indicates antiferromagnetic coupling in the electronic ground state, it is too small by a factor of about 44. Inclusion of dynamic effects by means of NEVPT2 leads to a significant numerical improvement ($J(1) = -23.2$ cm⁻¹) but is by far not sufficient to remedy the fundamental shortcomings of the underlying CASSCF(2,2) model. In case of complex **2**, the minimal active space even leads to a qualitatively wrong picture: Both, CASSCF and NEVPT2, predict a triplet ground state with small positive exchange coupling constants of 7.6 cm⁻¹ and 6.5 cm⁻¹, respectively.
<table><caption>**Table 4** Exchange coupling constants *J* (cm⁻¹) from DMRG-SCF/CI and DMRG-NEVPT2 calculations for complexes **1** and **2**</caption><thead><tr><th rowspan="2">Active space</th><th rowspan="2">Orbital composition</th><th colspan="2">J(1)</th><th colspan="2">J(2)</th></tr><tr><th>DMRG-SCF/CI</th><th>DMRG-NEVPT2</th><th>DMRG-SCF/CI</th><th>DMRG-NEVPT2</th></tr></thead><tbody><tr><td>(14,18)</td><td>B, D, E, F</td><td>-8.8</td><td>-38.2</td><td>5.3</td><td>-3.8</td></tr><tr><td>(14,24)</td><td>B, D, E, F, G</td><td>-7.5</td><td>-2.8</td><td>6.3</td><td>-1.3</td></tr><tr><td>(22,12)</td><td>A, B, D</td><td>-8.1</td><td>-47.9</td><td>4.5</td><td>-0.3</td></tr><tr><td>(22,14)</td><td>A, B, D, E</td><td>-8.8</td><td>-43.4</td><td>5.3</td><td>-1.1</td></tr><tr><td>(22,20)</td><td>A, B, D, E, G</td><td>-7.5</td><td>-36.8</td><td>6.4</td><td>0.9</td></tr><tr><td>(26,14)</td><td>A, C, D</td><td>-4.3</td><td>-37.9</td><td>7.6</td><td>-3.5</td></tr><tr><td>(26,16)</td><td>A, C, D, E</td><td>-4.5</td><td>-29.8</td><td>7.5</td><td>-3.9</td></tr><tr><td>(30,26)</td><td>B, C, D, E, F</td><td>-16.7</td><td>40.4</td><td>2.0</td><td>-0.1</td></tr><tr><td>(30,32)</td><td>B, C, D, E, F, G</td><td>-23.3</td><td>-</td><td>5.1</td><td>-</td></tr><tr><td>(38,30)</td><td>A, B, C, D, E, F</td><td>-20.7</td><td>-</td><td>-2.4</td><td>-</td></tr><tr><td>(38,36)</td><td>A, B, C, D, E, F, G</td><td>-30.9</td><td>-</td><td>4.9</td><td>-</td></tr></tbody></table>

Inclusion of the remaining Cu-3d orbitals (group C) in the active space has only a marginal effect on the relative CAS-SCF spin-state energies resulting in a shift below 1.0 cm⁻¹ of the corresponding exchange coupling constants. For NEVPT2, slightly larger effects are observed, shifting the exchange constants by -3.0 cm⁻¹ and -7.7 cm⁻¹, respectively. As a result, *J*(2) has a negative value, correctly corresponding to antiferromagnetic coupling. A similar shift of relative spin-state energies is observed when orbital group C is added to groups B, D, E and F, thereby enlarging the active space from (14,18) to (30,26). In the latter case, the NEVPT2 derived exchange coupling constants for **1** and **2** change by -2.2 cm⁻¹ and 3.7 cm⁻¹, respectively. In our view, these results emphasize the separation of the “doubly occupied” Cu 3d orbitals from the frontier orbitals and their concomitant subordinate role in the present context.

According to Anderson's model, superexchange is mediated by the bridging ligands between two metal-centered spins [87-89]. Therefore, the inclusion of ligand centered orbital group B in the active space is anticipated to considerably affect the calculated relative spin-state energies. In fact, addition of group B (O-2p orbitals) to the minimal active space results in a considerable lowering of all calculated exchange coupling constants. While the CASSCF derived exchange coupling of **1** undergoes only a slight change to a value of *J*(1) = -8.0 cm⁻¹, the corresponding NEVPT2 derived value is significantly lowered to *J*(1) = -42.4 cm⁻¹. This finding is in line with the results of a previous work by some of us on a bis-μ-oxo-bridged Mn dimer [6]. To our surprise, however, the predicted exchange coupling constants for **2** remain largely unaffected by inclusion of orbital group B to the active space. Interestingly, addition of orbital group G that incorporates O-3p dominated molecular orbitals results in a minor increase of the calculated exchange coupling constants. For example, upon an active space enlargement from (14,18) with groups B, D, E and F to (14,24)

with groups B, D, E, F and G the NEVPT2 derived *J*(1) increases by 5.4 cm⁻¹ to a value of -32.8 cm⁻¹. Similarly, *J*(2) increases by 2.5 cm⁻¹ to a value of -1.3 cm⁻¹.

In many instances, it has been found that radial correlation in the form of the so-called double-shell effect has a non-negligible impact on the energetics of late first-row transition metal complexes. Particularly when second-order perturbation theory is used to describe dynamic electron correlation on top of a CASSCF-Ansatz, the inclusion of a second set of transition metal d-orbitals to the active space has been found to be an important factor on the way to achieving accurate predictions. When orbital groups E and F that contain molecular orbitals with mostly Cu 4d character are appended to the (14,8) active space with orbital groups B and D, only minor changes between -3.8 cm⁻¹ and +3.8 cm⁻¹ are observed for the calculated exchange coupling constants. Further addition of the “doubly occupied” Cu-3d orbitals (group C) does not lead to significant changes. Thus, it appears that the double-shell effect plays only a subordinate role for antiferromagnetic coupling in **1** and **2**.

The leading configuration of both, the singlet and triplet state of **1** and **2**, features singly occupied orbitals of mostly Cu 3 $d_{x^2-y^2}$ character. Accordingly, the spin density is located along the Cu-O and the Cu-N bond axes. In view of this, spin distribution orbital group A consisting of Cu-N σ-bonding orbitals has been investigated for its role in mediating antiferromagnetic coupling in **1** and **2**. Notably, the active space consisting of orbital groups A, B and D yields the lowest *J*(1) value of -47.9 cm⁻¹ on the NEVPT2 level. In case of **2**, a value of *J*(2) = -3.5 cm⁻¹ is close to the lowest value of -3.9 cm⁻¹ that is obtained with a (26,16) active space consisting of orbital groups A, C, D and E.

The foregoing discussion allows us to identify orbital groups A, B and D as somewhat important for the mediation of antiferromagnetic coupling in **1** and **2**, while orbital
groups C, E, F and G contribute less strongly to the observed phenomenon. However, none of the presented exchange coupling constants in Table 4 is numerically close to the experimentally determined values of −180 cm⁻¹ and −45 cm⁻¹. Even with the largest affordable (30, 26) active space, the NEVPT2 results (−40.4 cm⁻¹ and −0.1 cm⁻¹) differ significantly from these values. For larger active spaces, only DMRG-SCF/CI calculations are feasible. Although for **1** a considerable negative shift of *J* for the (30, 32), (38, 30) and (38, 36) active spaces is observed, their relative size amounts to only 12–17% of the experimental value. These percentages are notably below the ~50% that have been observed for the previously studied bis-μ-oxo-bridged Mn dimer for which numerically accurate results could be obtained on the DMRG-SCF + NEVPT2 level of theory.

We thus conclude that the observed antiferromagnetic coupling in **1** and **2** is the result of the cumulative effect of a large number of configurations involving many orbitals that are not necessarily located around the bridging [Cu₂O₂] motif. In other words, dynamic electron correlation effects largely outweigh static correlation effects in the present cases. Finally, it should be noted that an enlargement of the active space does not necessarily induce a lowering of the calculated exchange coupling constants. In multiple cases, the addition of orbital groups results in a positive shift of the predicted values. This demonstrates the subtle connection between the contributions from different orbital sets in the present cases.

### 3.4 D. Difference-dedicated configuration interaction

The previous section documented that even at the current limits of computational feasibility, a brute-force approach that relies on active space expansion converges too slowly to be of practical utility for the copper dimers under study. An alternative strategy that has been employed successfully for systems like the present complexes is based on the difference-dedicated configuration interaction method. In this case, only the minimal active space of magnetic orbitals is used as reference model space, i.e., CAS(2,2) for Cu(II) dimers, and the correlation energy is recovered variationally by taking progressively into account different classes of excitation. When one-hole and one-particle excitations are considered, the method is denoted as DDCI1, equivalent to CAS plus single excitations from inactive to active (1 *h*), active to virtual (1 *p*) and inactive to virtual (1 *h*–1*p*) orbitals. DDCI2 additionally accounts for 2 *h* and 2*p* excitations. Finally, DDCI3 incorporates 2 *h*–1*p* and 1 *h*–2*p* excitations. It is at this level that quantitative predictions for magnetic coupling are expected, although the lower DDCI2 level has occasionally been deemed satisfactory.

The results from DDCI calculations on the two complexes considered in the present study are presented in Table 5. Here, we also demonstrate the dependence of the results on the parameter $T_{sel}$, which sets the threshold for inclusion of excited configuration state functions in the variational procedure based on the strength of their interaction with the 0th-order approximations to the target states. The two major conclusions are that DDCI2 is inadequate and that DDCI3 is necessary to obtain reliable results. In addition, the $T_{sel}$ value needs to be considerably tightened in order to avoid erratic behavior (e.g., notice the positive *J* value for complex **2** with DDCI3 at the default $T_{sel}$ setting used in ORCA). Provided these conditions are met, which can happen only at great computational cost, DDCI3 provides values for the exchange coupling constants that are indeed the best approximations to the experimental values compared to all other methods reported in the present work. This is in line with past studies of exchange-coupled Cu(II) systems [90, 91]. The superiority of DDCI3 compared to the DMRG-NEVPT2 approach even when the latter is based on a very large active space presumably suggests that excitations important for capturing a significant amount of dynamic electron correlation in these Cu(II) complexes may be of higher order than captured by NEVPT2. On the other hand, it should be kept in mind that DDCI at this level is hardly applicable, if at all, to higher-nuclearity complexes or to dimers with several unpaired electrons, and therefore, it cannot be viewed as a generally applicable approach in the same sense that DMRG-based methods are.

### 3.5 E. Broken-symmetry coupled cluster theory

In the previous sections, we showed that DFT is a computationally accessible approach yet unreliable as a predictive tool, and that DMRG-based CI and NEVPT2 fall short even with relatively large active spaces. Given that DDCI3 is not a generally applicable method, it is important to ask whether one can find an alternative approach that combines the simplicity of broken-symmetry DFT but brings in the possibility to include dynamic correlation at the level of a reliable wave function theory. This is in principle possible

<table><caption>**Table 5** Exchange coupling constants (in cm⁻¹) for complexes **1** and **2** derived from MR-DDCI(2,2) calculations, comparing the effect of the excitation level and the $T_{sel}$ parameter</caption><thead><tr><th rowspan="2">T<sub>sel</sub></th><th colspan="2">DDCI2</th><th colspan="2">DDCI3</th></tr><tr><th>J(1)</th><th>J(2)</th><th>J(1)</th><th>J(2)</th></tr></thead><tbody><tr><td>1.00 × 10<sup>-6</sup></td><td>6.9</td><td>80.8</td><td>-132.6</td><td>74.9</td></tr><tr><td>1.00 × 10<sup>-7</sup></td><td>-44.3</td><td>15.2</td><td>-181.6</td><td>-41.8</td></tr><tr><td>1.00 × 10<sup>-8</sup></td><td>-34.5</td><td>5.0</td><td>-176.3</td><td>-38.3</td></tr></tbody></table>
with broken-symmetry coupled cluster theory (BS-CC), which is only starting to be explored as a practical tool for the study of exchange-coupled systems.

Schurkus et al. have previously shown that the broken-symmetry approach in conjunction with the coupled cluster method produces useful results for bridged transition metal dimers [27]. This approach relied on using the highest spin and the broken-symmetry state to calculate *J* via the Yamaguchi approach (eq. (4)). However, the canonical coupled cluster procedure is expensive for realistically sized molecular systems. Several methods have been developed that can approximate the canonical solution within a reasonable accuracy at a reduced computational cost. One such approximation is based on the concept of local pair-natural orbitals (LPNOs) [45, 46]. Here, we employ and investigate, for the first time, the LPNO implementation of CCSD as a tool for the study of exchange coupling in a broken-symmetry approach.

The reasons that approximate methods like the LPNO are efficient include the localization of internal orbitals, due to which the number of electron pairs to be correlated is reduced, and the use of PNOs that span the virtual space and “compress” it. The accuracy of the LPNO results depends on two user-controlled cutoff parameters: *T*<sub>CutPNO</sub> and *T*<sub>CutPair</sub>. *T*<sub>CutPNO</sub> is the occupation number at which PNO expansion of a given pair is terminated. Thus, it controls the number of PNOs per pair of localized internal orbitals to be kept. It has a default value of 3.37 × 10<sup>-7</sup>. *T*<sub>CutPair</sub> is the parameter that controls which electron pairs to include based on their estimated pair-correlation energy. The pairs with the estimated pair-correlation energy lower than *T*<sub>CutPair</sub> are excluded. It has a default value of 10<sup>-4</sup> Eh. The pair-correlation energy of weakly correlated pairs is calculated by either semi-canonical or full iterative MP2. Collectively, the above settings are associated with the keywords NormalPNO and Tight-PNO, respectively. Both of these keywords are truncation thresholds and have their respective default values of the cutoff parameters as described in the work of Liakos et al. [92] Nevertheless, the values of these cutoff parameters can still be changed depending on the type of problem being studied. Unless otherwise indicated, the present calculations were performed with the TightPNO keyword as default starting point and with additional modifications of specific parameters.

In this section, we present the calculation of *J* with all of the three equations (eqs. (2) – (4)) mentioned in section B. In these calculations, *E*<sub>HS</sub> and *E*<sub>BS</sub> are the high-spin and broken-symmetry energies from the CCSD calculations, respectively. *S*<sub>max</sub> is the total spin of the high-spin state, and ⟨*S*<sup>2</sup>⟩<sub>HS</sub> and ⟨*S*<sup>2</sup>⟩<sub>BS</sub> are the spin expectation values for the high-spin and broken-symmetry HF states, respectively.

As an initial evaluation of the approach, we first report in Table 6 results of calculations done on one of the

molecules studied by Schurkus et al. that is comparable to molecule 1 and 2 in its molecular structure but is of smaller size (see Fig. 5). Compound 3 is tested for the behavior of *J* by changing the cutoff parameters *T*<sub>CutPNO</sub> and *T*<sub>CutPair</sub>. The purpose is to find out the limit of these thresholds at which the results of the LPNO-CCSD converge to the canonical CCSD result. The experimentally calculated *J* for 3 is −19 cm<sup>-1</sup>. Schurkus et al. calculated *J* with the canonical CCSD and CCSD(T) methods using the cc-pVDZ basis set, where they employed the Yamaguchi formulation and found the results to be equal to 4 cm<sup>-1</sup> and −12 cm<sup>-1</sup>, respectively. The result of our canonical

<figure><img src="image_6.png" /><figcaption>**Fig. 5** Compound 3 for which Schurkus et al. calculated the exchange coupling constant with broken-symmetry CCSD and CCSD(T) approaches [10].</figcaption></figure>

<table><caption>**Table 6** Effect of cutoff parameters *T*<sub>CutPNO</sub> and *T*<sub>CutPair</sub> on BS-CC computed exchange coupling constants (cm<sup>-1</sup>) for compound 3</caption><thead><tr><th><i>T</i><sub>CutPNO</sub></th><th><i>T</i><sub>CutPair</sub></th><th><i>J</i><sub>1</sub></th><th><i>J</i><sub>2</sub></th><th><i>J</i><sub>3</sub></th></tr></thead><tbody><tr><td>3.37 × 10<sup>-7</sup> (NormalPNO)</td><td>1.00 × 10<sup>-4</sup></td><td>-1104.5</td><td>-552.3</td><td>-1103.7</td></tr><tr><td>1.00 × 10<sup>-7</sup></td><td>1.00 × 10<sup>-5</sup></td><td>-612.4</td><td>-306.2</td><td>-611.9</td></tr><tr><td colspan="5"><i>Series 1</i></td></tr><tr><td>1.00 × 10<sup>-8</sup></td><td>1.00 × 10<sup>-5</sup></td><td>-508.5</td><td>-254.2</td><td>-508.1</td></tr><tr><td>1.00 × 10<sup>-8</sup></td><td>1.00 × 10<sup>-6</sup></td><td>-206.0</td><td>-103.0</td><td>-205.9</td></tr><tr><td>1.00 × 10<sup>-8</sup></td><td>1.00 × 10<sup>-7</sup></td><td>-149.8</td><td>-74.9</td><td>-149.7</td></tr><tr><td>1.00 × 10<sup>-8</sup></td><td>1.00 × 10<sup>-8</sup></td><td>-149.8</td><td>-74.9</td><td>-149.7</td></tr><tr><td colspan="5"><i>Series 2</i></td></tr><tr><td>1.00 × 10<sup>-9</sup></td><td>1.00 × 10<sup>-5</sup></td><td>-534.4</td><td>-267.2</td><td>-534.1</td></tr><tr><td>1.00 × 10<sup>-9</sup></td><td>1.00 × 10<sup>-6</sup></td><td>-183.1</td><td>-91.5</td><td>-182.9</td></tr><tr><td>1.00 × 10<sup>-9</sup></td><td>1.00 × 10<sup>-7</sup></td><td>-123.8</td><td>-61.9</td><td>-123.7</td></tr><tr><td>1.00 × 10<sup>-9</sup></td><td>1.00 × 10<sup>-8</sup></td><td>-123.6</td><td>-61.8</td><td>-123.5</td></tr><tr><td colspan="5"><i>Series 3</i></td></tr><tr><td>1.00 × 10<sup>-10</sup></td><td>1.00 × 10<sup>-5</sup></td><td>-520.8</td><td>-260.4</td><td>-520.4</td></tr><tr><td>1.00 × 10<sup>-10</sup></td><td>1.00 × 10<sup>-6</sup></td><td>-159.2</td><td>-79.6</td><td>-159.1</td></tr><tr><td>1.00 × 10<sup>-10</sup></td><td>1.00 × 10<sup>-7</sup></td><td>-29.1</td><td>-14.6</td><td>-29.1</td></tr><tr><td>1.00 × 10<sup>-10</sup></td><td>1.00 × 10<sup>-8</sup></td><td>-29.0</td><td>-14.5</td><td>-28.9</td></tr><tr><td colspan="5"><i>Series 4</i></td></tr><tr><td>1.00 × 10<sup>-11</sup></td><td>1.00 × 10<sup>-5</sup></td><td>-530.0</td><td>-265.0</td><td>-529.6</td></tr><tr><td>1.00 × 10<sup>-11</sup></td><td>1.00 × 10<sup>-6</sup></td><td>86.7</td><td>43.4</td><td>86.6</td></tr><tr><td>1.00 × 10<sup>-11</sup></td><td>1.00 × 10<sup>-7</sup></td><td>-29.9</td><td>-14.9</td><td>-29.8</td></tr><tr><td>1.00 × 10<sup>-11</sup></td><td>1.00 × 10<sup>-8</sup></td><td>-29.7</td><td>-14.8</td><td>-29.7</td></tr><tr><td>Canonical CCSD</td><td></td><td>21.0</td><td>10.5</td><td>21.0</td></tr><tr><td>Experimental</td><td></td><td></td><td>-19</td><td></td></tr></tbody></table>
CCSD calculation using the same formulation but with the def2-TZVP basis set was $20 \text{ cm}^{-1}$.

In the following discussion, all values for *J* were obtained using the Yamaguchi formulation since it is the generalized formula to calculate the magnetic exchange coupling constant. A first set of calculations was done invoking the NormalPNO and TightPNO keywords and keeping the default values of the cutoff parameters associated with those keywords. The *J*<sub>3</sub> obtained with NormalPNO cutoff parameter values was $-1103.7 \text{ cm}^{-1}$ which considerably deviates from the experimental value. Upon using the TightPNO cutoff parameter values, a significant shift in the value of *J*<sub>3</sub> is observed, but it still falls very far from the canonical value. We find that the default set of parameters associated with both NormalPNO and TightPNO are entirely unrealistic for attaining the accuracy required to describe the magnetic exchange coupling constant.

These initial calculations are followed by a systematic series of calculations, where $T_{\text{CutPNO}}$ is held constant at values between $1.0 \times 10^{-8}$ (series 1 in Table 6) and $1.0 \times 10^{-11}$ (series 4), while $T_{\text{CutPair}}$ is varied to inspect its effect on the predicted *J*. In series 1, as the $T_{\text{CutPair}}$ is tightened, the calculated *J*<sub>3</sub> value converges to $-149.7 \text{ cm}^{-1}$. A similar pattern is observed in series 2 as well where the *J*<sub>3</sub> value converges to $-123.5 \text{ cm}^{-1}$. The final results for both these series point to the fact that with the cutoff parameter values selected for these series, *J*<sub>3</sub> does not ultimately converge to the canonical or the experimental value. For series 3 and 4, *J*<sub>3</sub> converges to $-28.9 \text{ cm}^{-1}$ and $-29.7 \text{ cm}^{-1}$, respectively. These results suggest that by tightening the cutoff parameters, it is possible to converge *J* reasonably close to the canonical values, although the latter are not exactly reproduced. Of course, the computational cost increases considerably when cutoff parameters are tightened. In case of **1** and **2**, this increase turns out to be problematic.

Table 7 contains the results of a similar set of broken-symmetry LPNO-CCSD calculations for compound **1**. In this case, however, no canonical CCSD calculation could be performed owing to its enormous computational cost. Generally, similar trends are observed as above. Vastly erroneous results are observed for calculations that use the default NormalPNO or TightPNO truncation thresholds, but results are shifted considerably toward the experimental value after tightening of these thresholds. Yet, only a subset of the above presented calculations was tractable for **1** as the computational cost increases notably with tighter thresholds. Therefore, even the most accurate broken-symmetry LPNO-CCSD calculations reported here yield exchange coupling constants that are too small by a factor of $\approx 2$. Similar results were obtained when B3LYP high-spin and broken-symmetry reference functions are utilized (see Supporting Information).

<table><caption>**Table 7** Effect of cutoff parameters $T_{\text{CutPNO}}$ and $T_{\text{CutPair}}$ on BS-CC computed exchange coupling constants ($\text{cm}^{-1}$) for complex **1**</caption><thead><tr><th>$T_{\text{CutPNO}}$</th><th>$T_{\text{CutPair}}$</th><th>$J_1$</th><th>$J_2$</th><th>$J_3$</th></tr></thead><tbody><tr><td>$3.37 \times 10^{-7}$<br>(NormalPNO)</td><td>$1.00 \times 10^{-4}$</td><td>$-2759.5$</td><td>$-1379.7$</td><td>$-2757.3$</td></tr><tr><td>$1.00 \times 10^{-7}$</td><td>$1.00 \times 10^{-5}$</td><td>$-1378.3$</td><td>$-689.1$</td><td>$-1377.2$</td></tr><tr><td colspan="5"><em>Series 1</em></td></tr><tr><td>$1.00 \times 10^{-8}$</td><td>$1.00 \times 10^{-5}$</td><td>$-1147.6$</td><td>$-573.8$</td><td>$-1146.6$</td></tr><tr><td>$1.00 \times 10^{-8}$</td><td>$1.00 \times 10^{-6}$</td><td>$-1006.7$</td><td>$-503.3$</td><td>$-1005.9$</td></tr><tr><td>$1.00 \times 10^{-8}$</td><td>$1.00 \times 10^{-7}$</td><td>$-604.1$</td><td>$-302.0$</td><td>$-603.6$</td></tr><tr><td colspan="5"><em>Series 2</em></td></tr><tr><td>$1.00 \times 10^{-9}$</td><td>$1.00 \times 10^{-5}$</td><td>$-1086.1$</td><td>$-543.1$</td><td>$-1085.3$</td></tr><tr><td>$1.00 \times 10^{-9}$</td><td>$1.00 \times 10^{-6}$</td><td>$-911.2$</td><td>$-455.6$</td><td>$-910.5$</td></tr></tbody></table>

We conclude that for realistically sized molecules, broken-symmetry LPNO-CCSD has the theoretical possibility of approximating the canonical results, but the cost becomes astronomical as the thresholds are tightened. In addition, the currently available implementation does not include triples corrections, which inherently limits the accuracy of results. The approach remains promising, and it is suggested that it should be further explored in future work. However, in its current form it cannot be recommended as an approach of practical utility for systems such as those studied in the present work.

## 4 Conclusions

In this work, all available approaches to the calculation of exchange coupling constants are tested for a pair of antiferromagnetically coupled bis-μ-hydroxo Cu(II) dimers. The well-established broken-symmetry DFT Ansatz yields moderately accurate results low computational cost when the TPSSh functional is used together with Bencini's formula to extract the exchange coupling constant from the “raw” DFT energies. However, this combination does not perform universally well for all first-row transition metal systems. Instead, on account of the strong dependence of the obtained results on the functional form and the chosen formalism to evaluate exchange coupling constants, any DFT approach will have to be carefully benchmarked on a number of archetypical systems before it can be used in a predictive way.

In contrast to previous studies of early transition metal complexes, DMRG-based multireference calculations failed to provide accurate results for compounds **1** and **2** even when large active spaces were used and/or dynamic electron correlation was taken into account by second-order perturbation theory. Hence, the observed antiferromagnetism appears to be dominated by dynamic electron correlation effects that require an accurate description
beyond second-order perturbation theory of excited con-
figurations that are not covered by the large active spaces
described in this work. Our DDCI calculations on top of a
minimal active space confirm this assumption as they yield
exchange coupling constants close to the experimentally
determined values provided that tight selection thresholds
are employed and excitations with 3 degrees of freedom
are taken into account (DDCI3). Unfortunately, however,
this computational approach is not practical for larger sys-
tems and/or systems with many unpaired electrons due to
the intrinsically high computational costs.

In principle, the BS-LPNO-CCSD tested in this
work allows for a high-level treatment of dynamic elec-
tron correlation while also accounting for static elec-
tron correlation effects. However, test calculations on
$[Cu_2Cl_6]^{4-}$ revealed that only when extremely tight thresh-
olds for pair selection and pair-natural orbital selection in
the PNO framework are invoked, accurate results for the
exchange coupling constants can be expected. As a result,
the associated computational cost renders this approach
unfeasible even for the moderately sized compounds **1** and
**2**.

This work demonstrates that none of the tested set of
theoretical methods provides a fully satisfying approach to
the description of exchange coupling in systems with two
or more Cu(II) ions. Yet, it highlights the elevated relative
importance of dynamic electron correlation effects in this
kind of systems as compared to previously studied early
transition metal systems. Therefore, methodological devel-
opments aiming at efficient formulations of approximate
MR-CI (or even MR-CC) will likely be more useful in
the current context rather than those aiming at combina-
tions of large active spaces and second-order perturbation
theory.

**Supplementary Information** The online version contains supplementary material available at https://doi.org/10.1007/s00214-021-02830-0.

**Acknowledgements** This work was funded by the Agence nationale de la recherche (ANR) and the Deutsche Forschungsgemeinschaft (DFG) through the joint project CuBISM (ANR grant CUBISM, grant no. ANR-18 CE092 0040 01 and DFG grant no. RO 5688/3-1). MR would also like to thank the DFG for funding through project RO 5688/1-1. The authors also acknowledge financial support of this work by the France-Germany Hubert Curien Program—German Academic Exchange Service (DAAD) (Procope 2019-2020 project 42525PB/DAAD project 57445526).

**Funding** Open Access funding enabled and organized by Projekt DEAL. This work was funded by the Agence nationale de la recherche (ANR) and the Deutsche Forschungsgemeinschaft (DFG) through the joint project CuBISM (ANR grant CUBISM, grant no. ANR-18 CE092 0040 01 and DFG grant no. RO 5688/3-1). MR would also like to thank the DFG for funding through project RO 5688/1-1. The authors also acknowledge financial support of this work by the France-Germany Hubert Curien Program—German Academic Exchange Service

(DAAD) (Procope 2019–2020 project 42525 PB/DAAD project 57445526).

**Declarations**

**Conflicts of interest** The authors declare that they have no conflict of interest.

**Data Availability** Cartesian coordinates, geometrical parameters, additional BS-DFT results and additional BS-CC results are provided in Supporting Information.

**References**

1. Heisenberg W (1926) Mehrkorperproblem und Resonanz in der Quantenmechanik. Z Phys 38:411–426

2. Heisenberg WZ (1928) Zur theorie des ferromagnetismus. Z Phys 49:619–636

3. Dirac PAM (1929) Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 123:714–733

4. Vleck JHV (1932) The theory of electric and magnetic susceptibilities. Oxford University Press, Oxford

5. Malrieu JP, Caballol R, Calzado CJ, de Graaf C, Guihéry N (2014) Magnetic interactions in molecules and highly correlated materials: physical content, analytical derivation, and rigorous extraction of magnetic Hamiltonians. Chem Rev 114:429–492

6. Roemelt M, Krewald V, Pantazis DA (2018) Exchange coupling interactions from the density matrix renormalization group and N-electron valence perturbation theory: application to a biomimetic mixed-valence manganese complex. J Chem Theory Comput 14:166–179

7. Sharma P, Truhlar DG, Gagliardi L (2020) Magnetic coupling in a tris-hydroxo-bridged chromium dimer occurs through ligand mediated superexchange in conjunction with through-space coupling. J Am Chem Soc 142:16644–16650

8. Jezierski B (2010) Multireference coupled-cluster Ansatz. Mol Phys 108:3043–3054

9. Szalay PG, Müller T, Gidofalvi G, Lischka H, Shepard R (2012) Multiconfiguration self-consistent field and multireference configuration interaction methods and applications. Chem Rev 112:108–181

10. Lischka H, Nachtigallova D, Aquino AJ, Szalay PG, Plasser F, Machado FB, Barbatti M (2018) Multireference approaches for excited states of molecules. Chem Rev 118:7293–7361

11. Park JW, Al-Saadon R, MacLeod MK, Shiozaki T, Vlaisavljevich B (2020) Multireference electron correlation methods: journeys along potential energy surfaces. Chem Rev 120:5878–5909

12. Castell O, Caballol R (1999) Ab Initio configuration interaction calculation of the exchange coupling constant in hydroxo doubly bridged Cr(III) dimers. Inorg Chem 38:668–673

13. Calzado CJ, Cabrero J, Malrieu JP, Caballol R (2002) Analysis of the magnetic coupling in binuclear complexes. I. Physics of the coupling. J Chem Phys 116:2728–2747

14. Calzado CJ, Cabrero J, Malrieu JP, Caballol R (2002) Analysis of the magnetic coupling in binuclear complexes. II. Derivation of valence effective Hamiltonians from ab initio CI and DFT calculations. J Chem Phys 116:3985–4000

15. Miralles J, Daudey J-P, Caballol R (1992) Variational calculation of small energy differences. the singlet-triplet gap in $[Cu_2Cl_6]_2$. Chem Phys Lett 198:555–562
16. Miralles J, Castell O, Caballol R, Malrieu J-P (1993) Specific CI calculation of energy differences: transition energies and bond energies. Chem Phys 172:33–43

17. Chan GK-L, Head-Gordon M (2002) Highly correlated calculations with a polynomial cost algorithm: a study of the density matrix renormalization group. J Chem Phys 116:4462

18. Legeza Ö, Röder J, Hess BA (2003) Controlling the accuracy of the density-matrix renormalization-group method: the dynamical block state selection approach. Phys Rev B 67:125114

19. Chan GK-L, Sharma S (2011) The density matrix renormalization group in quantum chemistry. Annu Rev Phys Chem 62:465–481

20. Harris TV, Kurashige Y, Yanai T, Morokuma K (2014) Ab initio density matrix renormalization group study of magnetic coupling in dinuclear iron and chromium complexes. J Chem Phys 140:054303

21. Sharma S, Sivalingam K, Neese F, Chan GK-L (2014) Low-energy spectrum of iron-sulfur clusters directly from many-particle quantum mechanics. Nat Chem 6:927–933

22. Roemelt M, Pantazis DA (2019) Multireference approaches to spin-state energetics of transition metal complexes utilizing the density matrix renormalization group. Advanced Theory and Simulations 2:1800201

23. V. Krewald, D. A. Pantazis, (2019) In Transition Metals in Coordination Environments, Springer, pp. 91–120.

24. Noodleman L (1981) Valence bond description of antiferromagnetic coupling in transition metal dimers. J Chem Phys 74:5737

25. Caballol R, Castell O, Illas F, I. de P. R. Moreira, J. P. Malrieu, (1997) Remarks on the proper use of the broken symmetry approach to magnetic coupling. J Phys Chem A 101:7860–7866

26. Neese F (2009) Prediction of molecular properties and molecular spectroscopy with density functional theory: from fundamental theory to exchange-coupling. Coord Chem Rev 253:526–563

27. Schurkus H, Chen D-T, Cheng H-P, Chan G, Stanton J (2020) Theoretical prediction of magnetic exchange coupling constants from broken-symmetry coupled cluster calculations. J Chem Phys 152:234115

28. Olivares-Amaya R, Hu W, Nakatani N, Sharma S, Yang J, Chan GK-L (2015) The ab-initio density matrix renormalization group in practice. J Chem Phys 142:034102

29. Zgid D, Nooijen M (2008) The density matrix renormalization group self-consistent field method: orbital optimization with the density matrix renormalization group method in the active space. J Chem Phys 128:144116

30. Ghosh D, Hachmann J, Yanai T, Chan GK-L (2008) Orbital optimization in the density matrix renormalization group, with applications to polyenes and (β-carotene. J Chem Phys 128:144117

31. Angeli C, Cimiraglia R, Malrieu JP (2002) n-electron valence state perturbation theory: a spinless formulation and an efficient implementation of the strongly contracted and of the partially contracted variants. J Chem Phys 117:9138

32. Guo S, Watson MA, Hu W, Sun Q, Chan GK-L (2016) N-electron valence state perturbation theory based on a density matrix renormalization group reference function, with applications to the chromium dimer and a trimer model of poly(p-Phenylenevinylene). J Chem Theory Comput 12:1583–1591

33. Freitag L, Knecht S, Angeli C, Reiher M (2017) Multireference perturbation theory with cholesky decomposition for the density matrix renormalization group. J Chem Theory Comput 13:451

34. Khedkar A, Roemelt M (2019) Active space selection based on natural orbital occupation numbers from n-electron valence perturbation theory. J Chem Theory Comput 15:3522–3536

35. Neese F (2018) Software update: the ORCA program system, version 4.0. WIREs Comput. Mol. Sci. 8:e1327

36. Arcus C, Fivizzani K, Pavkovic S (1977) Preparation of Cu(tmen)OH(Cl) and molecular structure of

di-μ-hydroxo-bis(N,N,N',N'-tetramethylethylenediamine)-dicopper(II) perchlorate, [Cu(tmen)OH]2(ClO4)2. J Inorg Nucl Chem 39:285–287

37. Chaudhuri P, Ventur D, Wieghardt K, Peters EM, Peters K, Simon A (1985) Preparation, magnetism, and crystal structures of the tautomers [LCu(?2-OH)2CuL](ClO4)2 (Blue) and [LCu(?2-OH2) (?2-O)CuL](ClO4?2)2 (Green): ?-Aqua-?-oxo vs. Di-?-hydroxo linkage. Angew Chem Int Ed 24:57–59

38. Tao J, Perdew JP, Staroverov VN, Scuseria GE (2003) Climbing the density functional ladder: nonempirical meta-generalized gradient approximation designed for molecules and solids. Phys Rev Lett 91:146401

39. Weigend F, Ahlrichs R (2005) Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy. Phys Chem Chem Phys 7:3297–3305

40. Weigend F (2006) Accurate coulomb-fitting basis sets for H to Rn. Phys Chem Chem Phys 8:1057–1065

41. Dunlap BI, Connolly JWD, Sabin JR (1979) On some approximations in applications of X α theory. J Chem Phys 71:3396–3402

42. Vahtras O, Almlöf J, Feyereisen MW (1993) Integral approximations for LCAO-SCF calculations. Chem Phys Lett 213:514–518

43. Meyer W (1973) PNO–CI Studies of electron correlation effects. I. Configuration expansion by means of nonorthogonal orbitals, and application to the ground state and ionized states of methane. J Chem Phys 58:1017–1035

44. Ahlrichs R, Driessler F, Lischka H, Staemmler V, Kutzelnigg W (1975) PNO–CI (pair natural orbital configuration interaction) and CEPA–PNO (coupled electron pair approximation with pair natural orbitals) calculations of molecular systems. II. The molecules BeH2, BH, BH 3, CH4, CH–3, NH3 (planar and pyramidal), H2O, OH+3, HF and the Ne atom. J Chem Phys 62:1235–1247

45. Neese F, Wennmohs F, Hansen A (2009) Efficient and accurate local approximations to coupled-electron pair approaches: an attempt to revive the pair natural orbital method. J Chem Phys 130:114108

46. Neese F, Hansen A, Liakos DG (2009) Efficient and accurate approximations to the local coupled cluster singles doubles method using a truncated pair natural orbital basis. J Chem Phys 131:064103

47. Lee CT, Yang WT, Parr RG (1988) Development of the Colle-Salvetti correlation-energy formula into a functional of the electron density. Phys Rev B 37:785–789

48. Becke AD (1993) Density-functional thermochemistry. III. The role of exact exchange. J Chem Phys 98:5648–5652

49. Sharma S, Chan GK-L (2012) Spin-adapted density matrix renormalization group algorithms for quantum chemistry. J Chem Phys 136:124121

50. Neese F, Wennmohs F, Hansen A, Becker U (2009) Efficient, approximate and parallel Hartree–Fock and hybrid DFT calculations. a ‘chain-of-spheres’ algorithm for the Hartree–Fock exchange. Chem Phys 356:98–109

51. Izsák R, Neese F (2011) An overlap fitted chain of spheres exchange method. J Chem Phys 135:144105

52. Izsák R, Neese F, Klopper W (2013) Robust fitting techniques in the chain of spheres approximation to the fock exchange: the role of the complementary space. J Chem Phys 139:094111

53. Hättig C (2005) Optimization of auxiliary basis sets for RI-MP2 and RI-CC2 calculations: Core-valence and quintuple-ζ basis sets for H to Ar and QZVPP basis sets for Li to Kr. Phys Chem Chem Phys 7:59–66

54. Bencini A, Totti F, Daul CA, Doclo K, Fantucci P, Barone V (1997) Density functional calculations of magnetic exchange interactions in polynuclear transition metal complexes. Inorg Chem 36:5022–5030
55. Comba P, Hausberg S, Martin B (2009) Calculation of exchange coupling constants of transition metal complexes with DFT. J Phys Chem A 113:6751–6755

56. Ruiz E (2011) Exchange coupling constants using density func-
tional theory: Long-range corrected functionals. J Comput Chem
32:1998–2004

57. Luo S, Nie D, Li Z, Sun X, Hu L, Liu X (2020) Effects of car-
boxylic acid auxiliary ligands on the magnetic properties of azido-
Cu (II) complexes: A density functional theory study. Polyhedron
182:114506

58. Jana A, Konar S, Das K, Ray S, Golen JA, Rheingold AL, Carrella LM, Rentschler E, Mondal TK, Kar SK (2012) Azide bridged dicopper(II), dicobalt(II) complexes and a rare double µ-chloride bridged ferromagnetic dicobalt(II) complex of a pyrazolyl-pyrimidine ligand: Synthesis, crystal structures, magnetic and DFT studies. Polyhedron 38:258–266

59. Zaltariov MF, Vieru V, Zalibera M, Cazacu M, Martins NM, Martins LM, Rapta P, Novitchi G, Shova S, Pombeiro AJ (2017) A Bis (µ-chlorido)-bridged cobalt (II) complex with silyl-containing schiff base as a catalyst precursor in the solvent-free oxidation of cyclohexane. Eur J Inorg Chem 2017:4324–4332

60. Pantazis DA (2019) Meeting the challenge of magnetic coupling in a triply-bridged chromium dimer: complementary broken-symmetry density functional theory and multireference density matrix renormalization group perspectives. J Chem Theory Comput 15:938–948

61. Orio M, Pantazis DA, Petrenko T, Neese F (2009) Magnetic and spectroscopic properties of mixed valence manganese(III,IV) dimers: a systematic study using broken symmetry density functional theory. Inorg Chem 48:7251–7260

62. Joshi RP, Phillips JJ, Mitchell KJ, Christou G, Jackson KA, J. E. (2020) Accuracy of density functional theory methods for the calculation of magnetic exchange couplings in binuclear iron (III) complexes. Peralta, Polyhedron 176:114194

63. Grimme S (2005) Accurate calculation of the heats of formation for large main group compounds with spin-component scaled MP2 methods. J Phys Chem A 109:3067–3077

64. Quintal MM, Karton A, Iron MA, Boese AD, Martin JM (2006) Benchmark Study of DFT functionals for late-transition-metal reactions †. J Phys Chem A 110:709–716

65. Staroverov VN, Scuseria GE, Tao J, Perdew JP (2003) assessment of a new nonempirical density functional: molecules and hydrogen-bonded complexes. J Chem Phys 119:12129–12137

66. Zhao Y, Truhlar DG (2008) The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals. Theor Chem Acc 120:215–241

67. Ruiz E, Alvarez S, Cano J, Polo V (2005) About the calculation of exchange coupling constants using density-functional theory: The role of the self-interaction error. J Chem Phys 123:164110

68. Adamo C, Barone V, Bencini A, Broer R, Filatov M, Harrison N, Illas F, Malrieu J, I. de PR Moreira, (2006) About the calculation of exchange coupling constants using density-functional theory: The role of the self-interaction error. J Chem Phys 124:107101

69. Ruiz E, Cano J, Santiago Alvarez V, Polo, (2006) About the calcula-
tion of exchange coupling constants using density-functional theory:
The role of the self-interaction error. J Chem Phys 124:107102

70. Bencini A, Totti F (2009) A few comments on the application of density functional theory to the calculation of the magnetic structure of oligo-nuclear transition metal clusters. J Chem Theory Comput 5:144–154

71. Ruiz E, Cano J, Alvarez S, Alemany P (1999) Broken symmetry approach to calculation of exchange coupling constants for homobinuclear and heterobinuclear transition metal complexes. J Comput Chem 20:1391–1400

72. Yamaguchi K, Fukui H, Fueno T (1986) Molecular orbital (MO) theory for magnetically interacting organic compounds. ab-initio mo calculations of the effective exchange integrals for cyclophane-type carbene dimers. Chem Lett 15:625–628

73. Santiago PH, Tiago FS, Castro MS, Souza PE, Martins JB, Gatto CC (2020) DFT analysis, spectroscopic study and biological activity of a newly synthesized benzoylhydrazone binuclear Cu (II) complex. J Inorg Biochem 204:110949

74. Luo S, Mei H, Sun X, Zheng P (2020) Effect of 3d heterometallic ions on the magnetic properties of azido-Cu (II) with isonicotinic acid coligands: A theoretical perspective. J. Mol. Graphics Modell 97:107562

75. Ruiz E, Alemany P, Alvarez S, Cano J (1997) Toward the prediction of magnetic coupling in molecular systems: hydroxo- and alkoxo-bridged Cu(II) binuclear complexes. J Am Chem Soc 119:1297–1303

76. Ruiz E, Alemany P, Alvarez S, Cano J (1997) Structural modeling and magneto-structural correlations for hydroxo-bridged copper(II) binuclear complexes. Inorg Chem 36:3683–3688

77. Reis NV, Barros WP, Oliveira WX, Pereira CL, Rocha WR, Pinheiro CB, Lloret F, Julve M, Stumpf HO (2018) Crystal structure and magnetic properties of an oxamato-bridged heterobimetallic tetranuclear [NiIICuII] 2 complex of the rack type. Eur J Inorg Chem 2018:477–484

78. Simões TR, Walace D, Metz KC, Ribeiro MA, Valente DC, Cardozo TM, Pinheiro CB, Pedroso EF, Horta BA, Pereira CL (2020) Dinuclear copper(ii) complexes containing oxamate and blocking ligands: crystal structure, magnetic properties, and DFT calculations. New J Chem 44:2597–2608

79. Becke AD (1988) Density-functional exchange-energy approximation with correct asymptotic behavior. Phys Rev A 38:3098

80. Perdew JP, Wang Y (1992) Accurate and simple analytic representation of the electron-gas correlation energy. Phys Rev B 45:13244

81. Perdew JP, Burke K, Ernzerhof M (1996) Generalized gradient approximation made simple. Phys Rev Lett 77:3865–3868

82. Perdew JP (1986) Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys Rev B 33:8822

83. Adamo C, Barone V (1999) Toward reliable density functional methods without adjustable parameters: The PBE0 model. J Chem Phys 110:6158–6170

84. Ruiz E (2004) Principles and Applications of Density Functional Theory in Inorganic. Chemistry II:71–102

85. Stein CJ, Pantazis DA, Krewald V (2019) Orbital entanglement analysis of exchange-coupled systems. J Phys Chem Lett 10:6762–6770

86. Chalupský J, Srnec M, Yanai T (2021) Interpretation of exchange interaction through orbital entanglement. J Phys Chem Lett 12:1268–1274

87. Anderson P (1950) Antiferromagnetism. theory of superexchange interaction. Phys Rev 79:350

88. Anderson PW (1959) New approach to the theory of superexchange interactions. Phys Rev 115:2

89. P. W. Anderson (1963), In *Solid State Phys., Vol. 14*, Elsevier, pp. 99–214.

90. I. de PR Moreira, F. Illas, (2006) A unified view of the theoretical description of magnetic coupling in molecular chemistry and solid state physics. Phys Chem Chem Phys 8:1645–1659

91. Maurice R, Sivalingam K, Ganyushin D, Guihéry N, de Graaf C, Neese F (2011) Theoretical determination of the zero-field splitting in copper acetate monohydrate. Inorg Chem 50:6229–6236

92. Liakos DG, Sparta M, Kesharwani MK, Martin JM, Neese F (2015) Exploring the accuracy limits of local pair natural orbital coupled-cluster theory. J Chem Theory Comput 11:1525–1539

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.