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

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.


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][2][3][4] 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 Ŝ I and Ŝ J is associated with a coupling parameter J IJ . Negative values for J IJ 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 (1) H HDvV = −2 ∑ I<J J IJŜIŜJ 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 selfconsistent 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][6][7]. Well-established approaches to accomplish this crucial task are based on perturbation theory, configuration interaction or coupled cluster theory [8][9][10][11].
A common concern for all multireference methods is their high computational cost, arising from the underlying CAS-SCF 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][13][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][18][19]. Relying on a socalled 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][21][22][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 ⟨ S 2 ⟩ . 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 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.
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][39][40][41][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][44][45][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][51][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.

A. Description of the Cu dimers
The structures of 1 and 2 display the same coordination core [Cu 2 (μ-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).

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][55][56][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][68][69] regarding the correct way of calculating the exchange coupling constant from the high-spin and BS energies, since three different schemes can be employed within the BS-DFT framework (eqs. (2)- (5)). [70] These formalisms were described in works from Noodleman's (Eqs. (2) and (5) where E HS and E BS are the high-spin and broken-symmetry energies, respectively, S ab is the orbital overlap, S max is the total spin of the high-spin state, and ⟨ S 2 ⟩ HS and ⟨ S 2 ⟩ BS are the spin expectation values for the high-spin and brokensymmetry 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 orbitals [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-tostrong 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 equations 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 contrary, Luo et al. [74] and Ruiz et al. [75,76] described dinuclear 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 underestimated 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.
Prior to a detailed discussion of the numerical performance of the various tested functionals in the current context, it is instructive to examine the nature of antiferromagnetic 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 analyzed. For both complexes, the magnetic orbitals are the 3 d x 2 −y 2 Cu orbitals in σ* combination with ligand p orbitals. Exchange is mediated mostly by overlap over the bridging hydroxyl groups (2p x and 2p y 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 integral is correlated with the nature of the functional, specifically with the percentage of exact exchange. This in turn is correlated with the energetic separation between the BS and 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 −1 ) and absolute percentage devia-tion (APD) for the 4 different spin projection formalisms. The experimental J value is − 180 cm −1 Functionals % of HF exchange    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 1 and J 3 , 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 2 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.

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. Finally, groups E, F and G (Fig. 4) correspond to "secondshell" 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 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 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 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 −1 . 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 −1 ) 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 −1 and 6.5 cm −1 , respectively.  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 −1 of the corresponding exchange coupling constants. For NEVPT2, slightly larger effects are observed, shifting the exchange constants by − 3.0 cm −1 and − 7.7 cm − 1 , 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 −1 and 3.7 cm −1 , 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][88][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 −1 , the corresponding NEVPT2 derived value is significantly lowered to J(1) = − 42.4 cm −1 . 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 −1 to a value of − 32.8 cm −1 . Similarly, J(2) increases by 2.5 cm −1 to a value of − 1.3 cm −1 .
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 −1 and + 3.8 cm −1 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 −1 on the NEVPT2 level. In case of 2, a value of J(2) = − -3.5 cm −1 is close to the lowest value of − 3.9 cm −1 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  (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 2 O 2 ] 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.

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 (1p) and inactive to virtual (1 h-1p) orbitals. DDCI2 additionally accounts for 2 h and 2p excitations. Finally, DDCI3 incorporates 2 h-1p and 1 h-2p 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.

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 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 brokensymmetry 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 brokensymmetry 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 CutPNO and T CutPair . T CutPNO 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 -7 . T CutPair 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 CutPair are excluded. It has a default value of 10 -4 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 HS and E BS are the high-spin and broken-symmetry energies from the CCSD calculations, respectively. S max is the total spin of the high-spin state, and ⟨ S 2 ⟩ HS and ⟨ S 2 ⟩ BS 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 CutPNO and T CutPair . 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 −1 . 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 −1 and − 12 cm −1 , respectively. The result of our canonical CCSD calculation using the same formulation but with the def2-TZVP basis set was 20 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 Nor-malPNO and TightPNO keywords and keeping the default values of the cutoff parameters associated with those keywords. The J 3 obtained with NormalPNO cutoff parameter values was − 1103.7 cm −1 which considerably deviates from the experimental value. Upon using the TightPNO cutoff parameter values, a significant shift in the value of J 3 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 CutPNO is held constant at values between 1.0 × 10 -8 (series 1 in Table 6) and 1.0 × 10 -11 (series 4), while T CutPair is varied to inspect its effect on the predicted J. In series 1, as the T CutPair is tightened, the calculated J 3 value converges to − 149.7 cm −1 . A similar pattern is observed in series 2 as well where the J 3 value converges to − 123.5 cm −1 . The final results for both these series point to the fact that with the cutoff parameter values selected for these series, J 3 does not ultimately converge to the canonical or the experimental value. For series 3 and 4, J 3 converges to − 28.9 cm −1 and − 29.7 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 brokensymmetry 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 ≈ 2. Similar results were obtained when B3LYP high-spin and broken-symmetry reference functions are utilized (see Supporting Information).
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.

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 secondorder 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 configurations 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 systems 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 electron correlation while also accounting for static electron correlation effects. However, test calculations on [Cu 2 Cl 6 ] 4− revealed that only when extremely tight thresholds 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 developments 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 combinations of large active spaces and second-order perturbation theory.