Natalie
Orms
and
Anna I.
Krylov
*
Department of Chemistry, University of Southern California, Los Angeles, California 90089-0482, USA. E-mail: krylov@usc.edu
First published on 22nd January 2018
Molecular magnets, defined here as organic polyradicals, can be used as building blocks in the fabrication of novel and structurally diverse magnetic light-weight materials. We present a theoretical investigation of the lowest spin states of several binuclear copper diradicals. In contrast to previous studies, we consider not only the energetics of the low-lying states (which are related to the exchange-coupling parameter within the Heisenberg–Dirac–van-Vleck model), but also the character of the diradical states themselves. We use natural orbitals, their occupations, and the number of effectively unpaired electrons to quantify bonding patterns in these systems. We compare the performance of spin-flip time-dependent density functional theory (SF-TDDFT) using various functionals and effective core potentials against the wave function based approach, equation-of-motion spin-flip coupled-cluster method with single and double substitutions (EOM-SF-CCSD). We find that SF-TDDFT paired with the PBE50 and B5050LYP functionals performs comparably to EOM-SF-CCSD, with respect to both singlet–triplet gaps and states' characters. Visualization of frontier natural orbitals shows that the unpaired electrons are localized on copper centers, in some cases exhibiting slight through-bond interaction via copper d-orbitals and p-orbitals of neighboring ligand atoms. The analysis reveals considerable interactions between the formally unpaired electrons in the antiferromagnetic diradicaloids, meaning that they are poorly described by the Heisenberg–Dirac–van-Vleck model. Thus, for these systems the experimentally derived exchange-coupling parameters are not directly comparable with the singlet–triplet gaps. This explains systematic discrepancies between the computed singlet–triplet energy gaps and the exchange-coupling parameters extracted from experiment.
The triplet state of a diradical can be represented by either a single Slater determinant (high-spin M_{s} = 1 configuration (i)) or by a linear configuration of two M_{s} = 0 determinants with equal weights (configuration (ii)). Configurations (iii)–(v) give rise to three singlet states. The energy gaps, relative state ordering, and relative weights of the Slater determinants (i.e., coefficients λ) depend on the nature and energy separation of the respective frontier molecular orbitals (MOs) ϕ_{1} and ϕ_{2}. The character of the MOs also determines the character of the wave function, e.g., whether it is predominantly covalent (i.e., two electrons residing on different parts of a molecule) or ionic. Detailed analysis of different types of diradical electronic structure can be found in classic papers.^{7,8} In the context of MMs, the two states of interest are the lowest singlet and triplet states. In systems like binuclear copper complexes, one expects these two states to have covalent wave functions in which the unpaired electrons are localized on the two metal centers:
Ψ^{s,t}(1,2) ∼ [ϕ_{Cu1}(1)ϕ_{Cu2}(2) ± ϕ_{Cu2}(1)ϕ_{Cu1}(2)] × [α(1)β(2) ∓ β(1)α(2)] | (1) |
[ϕ_{Cu1}(1)ϕ_{Cu1}(2) + ϕ_{Cu2}(1)ϕ_{Cu2}(2)] × [α(1)β(2) − β(1)α(2)] | (2) |
In these expressions, ϕ_{Cu1} and ϕ_{Cu2} denote orbitals localized on the two copper centers, such as copper d-orbitals (perhaps including small contributions from the nearest ligand atoms). If the actual MOs hosting the unpaired electrons are delocalized and can be described as (nearly degenerate) bonding and antibonding combinations of ϕ_{Cu1} and ϕ_{Cu2} (case 1 in Fig. 1), as is the case in the MMs studied here, the triplet states are described by configurations (i) and (ii) and have pure covalent character, as Ψ^{t} in eqn (1). The character of the lowest singlet state can vary, depending on the exact weights of configurations (iii) through (v). A purely covalent singlet wave function, Ψ^{s} of eqn (1), corresponds to configuration (iv) with λ = 1. A smaller value of λ gives rise to the ionic configurations mixed into the wave function. This happens when the interaction between the two centers stabilizes the bonding MO relative to the antibonding one, either due to through-space or through-bond interactions. The ionic configurations can also appear in the singlet wave functions due to mixing with configuration (iii), which, in contrast to (ii), has pure ionic character. Since all configurations, (iii)–(v), can contribute into the singlet state, the correct description of this state requires an electronic structure method that treats (iii)–(v) on an equal footing. In order to describe relative energies of singlet and triplet states, the method should provide a balanced and unbiased description of all four M_{s} = 0 configurations from Fig. 1.
From a theoretical perspective, the search for promising MMs begins with first-principle calculations of the relevant terms in the phenomenological spin Hamiltonian.^{2,12,13} Of all terms in the spin Hamiltonian, the most energetically significant one is the exchange-coupling interaction between unpaired, spatially separated electrons.^{14–16} The sign and magnitude of electronic exchange-coupling between atomic spin centers determines whether a molecule will behave ferromagnetically or antiferromagnetically upon exposure to an external magnetic field, and thus, whether the molecule is suitable for application in a magnetic material. In bimetallic diradical complexes, like those considered here, the exchange-coupling constant equals the energy difference between the lowest singlet and triplet states.^{17} Thus, exchange-coupling terms of the spin Hamiltonian can be computed by ab initio methods as singlet–triplet energy gaps. This approach can be generalized to systems with more unpaired electrons and multiple polyradical centers.^{18,19}
The main challenge in applying this simple strategy is in the multiconfigurational character of the low-spin wave functions, which becomes evident by inspecting Fig. 1. The high-spin states, such as M_{s} = 1 triplet (i), are well represented by a single Slater determinant and, therefore, their energies can be reliably computed by standard single-reference methods, i.e., coupled-cluster theory or DFT. In contrast, wave functions of low-spin states (ii)–(v) require at least two Slater determinants. Consequently, they are poorly described by single-reference methods. This imbalance in the description of high-spin and low-spin states results in large errors in the computed singlet–triplet gaps.
Several strategies have been employed for describing the open-shell states and exchange-coupling in MMs. Historically, the most popular are the broken symmetry (BS) methods^{20–22} and spin-restricted Kohn–Sham (REKS/ROKS) methods.^{22–24} Both approaches suffer from imbalance in their treatment—and sometimes, outright exclusion—of important configurations depicted in Fig. 1. For example, in BS approach all singlet and triplet configurations are scrambled. While spin projection allows one to formally separate singlet and triplet manifolds, it does not distinguish between open- and closed-shell singlet states (iii)–(v). Despite belonging to entirely different states (which may even have different spatial symmetry), these configurations remain scrambled in spin-projected BS solutions. BS solutions are often contaminated by mixing with other electronic states, i.e., higher multiplets.^{25} Even more disturbing is the fact that spin-projectors are not uniquely defined and the results are affected by specific choices of a projector.^{20,26–29} Furthermore, BS-DFT does not scale well with the number of radical sites, as the number of BS solutions grows rapidly. The most important concern is the quality of wave functions and spin-dependent properties extracted from BS-DFT. As pointed out by Neese,^{12} although broken-symmetry solutions represent charge density reasonably well, they result in erroneous spin density. Consequently, the quality of broken-symmetry description of properties that are determined by charge density and by spin density can be vastly different. Thus, despite being reasonably successful in terms of predicting energy splittings between the different components of multiplet spectra, BS-DFT approach does not offer a fully satisfactory solution to modeling MMs. The drawback of the REKS/ROKS scheme is that closed- and open-shell states are described by different sets of equations and that the effect of different dynamic correlation in the two manifolds is not fully accounted for.^{30}
A more balanced and robust alternative is the spin-flip (SF) family of methods.^{31–36} Within SF approach, a single-determinant high-spin triplet state is used as a reference from which one computes the manifold of low-spin states arising from configurations (ii)–(iv) in Fig. 1 by using spin-flipping excitations. The performance of various variants of SF method for organic polyradicals has been extensively benchmarked, focusing on the energy gaps between the electronic states and their structures.^{9,34–42} Recently, we investigated the performance of SF methods in terms of characters of the underlying wave functions.^{11} Specifically, we compared how different functionals and wave function methods describe frontier orbitals, their occupations, and the number of effectively unpaired electrons. Towards this end, we employed density-based wave function analysis tools.^{43–45} By using natural orbitals and their occupations, this analysis allows one to map correlated multiconfigurational wave functions into a simple two-electrons-in-two-orbitals picture from Fig. 1.
The goal of this study is to investigate the performance of SF methods for MM. In contrast to previously studied organic di- and tri-radicals, typical singlet–triplet gaps in MM are much smaller, within several hundreds of wave numbers, which calls for sub-chemical accuracy (1 kcal mol^{−1} = 350 cm^{−1}). Additional challenges arise due to the presence of heavy atoms, which might require using special basis sets and/or effective core potentials (ECPs). Although encouraging results of SF calculations for several MMs with transition metals have been reported before,^{18,19,25,46–48} no systematic studies of the effect of ECP (and/or basis sets) on multiplet splittings within the SF-TDDFT framework has been carried out. In the present work, we focus on a set of eight binuclear copper complexes whose magnetic properties have been extensively studied.^{22,46,49–58} Given the above considerations and their large sizes, binuclear copper diradicals represent a stringent test for the theory. We report the results of the EOM-SF-CCSD and SF-TDDFT calculations of singlet–triplet gaps (which are equal to exchange-coupling within the assumptions of Heisenberg–Dirac–van-Vleck (HDvV) model^{14–16}) and the bonding pattern (the shapes of frontier natural orbitals and the extent of interaction between the formally unpaired electrons). We analyze the effects of basis sets, ECPs, TDDFT kernel, the functional, and molecular geometry.
The structure of the article is as follows: the next section provides a brief theoretical overview of SF methods, the HDvV model^{14–16} for exchange-coupling, and density-based analysis. The following section outlines computational details. We then illustrate the effects of the above mentioned variables on computed singlet–triplet gaps and wave function properties (frontier natural orbitals, their occupation, and the number of effectively unpaired electrons) of the eight binuclear copper compounds. In conclusion, we comment on the applicability of the HDvV model for this family of diradicals and the importance of density-based analysis in the description of open-shell systems.
To describe diradicals within the SF framework, one begins with the high-spin M_{s} = 1 triplet state, which is well represented by a single Slater determinant. The low-spin states (singlets and the M_{s} = 0 triplets) appear as single excitations from the high-spin reference state. Thus, in the SF approach, a high-spin state is used as the reference state and the target manifold of states is generated by applying a linear spin-flipping excitation operator to the:^{60}
(3) |
Using different approaches to describe correlation in the reference state gives rise to different SF methods.^{31–36,61} Our primary focus here is SF-TDDFT approach,^{34,35,62–64} which is very attractive owing to its low computational cost. We compare the performance of several functionals chosen on the basis of the results of the previous studies.^{34,35} For the benchmarking purposes, we also report the results of EOM-SF-CCSD,^{33} a highly accurate wave-function method including advanced treatment of dynamical correlation. A density-based wave function analysis allows quantitative comparison of the character of underlying wave functions computed by different methods.
(4) |
The parameters of eqn (4)—such as hyperfine couplings, zero-field splittings, exchange-couplings—can be extracted from experimental measurements or computed from correlated many-electron wave functions or DFT.^{2,12,13} All phenomenological terms in eqn (4) can be computed using existing relativistic and non-relativistic electronic structure methods.^{21,22,46,47,58,69–75} In this work, we focus on the third, non-relativistic term, which gives rise to the HDvV Hamiltonian:^{14–16}
(5) |
The HDvV Hamiltonian describes weakly coupled localized electrons neglecting hyperfine couplings and zero-field splittings. This phenomenological Hamiltonian can be mapped to the exact many-body Hamiltonian using effective Hamiltonian theory.^{13} By solving eqn (5), one can relate the energies of different spin states (singlets and triplets, in the case of diradicals) to the experimental magnetic parameter, J, known as the exchange-coupling constant.^{17} Exchange-couplings represent dressed, or effective interactions between the spins.^{13} The sign of J determines whether a molecular magnet behaves ferromagnetically or antiferromagnetically in an external magnetic field.
Analytic solutions of eqn (5) are used to construct a formula connecting macroscopic magnetic susceptibility with exchange-couplings.^{76} This formula is used to extract the latter from the temperature-dependent measurements of magnetic susceptibility. Thus, all experimental values of J derived from these measurements equal to the true energy gaps only if the key assumption of the HDvV model, that the spins are weakly interacting and the contributions of the ionic configurations into the total wave function are small, hold true.
Eqn (5) assumes two (or more) weakly interacting, spatially separated electronic spins localized on two (or more) atomic centers, A and B. Obviously, the configurations arising from paired (ionic) electronic configurations or locally excited (i.e., non-Hund) configurations are not admitted.^{77} The justification for this choice of the model Hamiltonian is that for the weakly interacting unpaired electrons, the ionic configurations are much higher in energy than covalent triplet and singlet configurations of eqn (1). Under these assumptions, the HDvV Hamiltonian can be derived from the Hubbard model by using quasidegenerate perturbation theory (see, for example, ref. 13).
What happens when ionic configurations are mixed with large weights into the singlet states in antiferromagnetic species? In this case, the energy of the singlet diradical state is stabilized via configuration interaction, leading to an increased singlet–triplet gap. While the macroscopic properties, such as temperature dependence of magnetic susceptibility will reflect this phenomenon, the application of the HDvV-based model to fit the data would produce J values that are too large, because larger values of J are necessary to describe the increased singlet–triplet gap in the absence of explicit ionic configurations in the model. Thus, one should expect that for strongly AF species the experimentally derived J values are systematically overestimated.
In the two-center systems, J is directly related to the energies of spin states via the Landé interval rule.^{17} With an electronic structure method that offers a balanced treatment of open and closed-shell states, one can calculate J as an energy difference between the spin states:
E(S) − E(S − 1) = −J_{AB}·S | (6) |
In diradicals S and S − 1 are triplet and singlet states: E(S) is the energy of the lowest-energy triplet state and E(S − 1) is the energy of the lowest singlet state. A positive J gives rise to a ferromagnetic ground state, while negative J gives rise to an antiferromagnetic ground state. Note that a different form of eqn (5) might be used,^{18,19} with a factor of 2 in front of the sum, giving rise to a different relationship between energy gaps and exchange-coupling. Thus, for meaningful comparisons between different studies, it is important to carefully check which form was used (and of course it is always safer to report energy gaps between physical states, e.g., ΔE_{ST} in the case of diradicals).
Fig. 2 Spin difference densities, unrestricted singly occupied molecular and natural frontier orbitals (SOMOs and SONOs, respectively) of triplet CUAQAC02 (left) and CITLAT (right) at the B5050LYP/cc-pVTZ level of theory. CUAQAC02 has 202 electrons while CITLAT has 278, making their low-lying states and associated orbital surfaces numerous and complex. Reproduced with permission from ref. 11. |
Fig. 3 Eight binuclear copper complexes included in this study. Structures are denoted with their Cambridge Structural Database names.^{107,108} Complexes 1, 3, 4, 6, and 7 have a charge of 2+; complex 2 is neutral; complex 5 has a charge of 2−; and complex 8 has a charge of 1+. The respective Lewis structures are shown in Fig. S1 in ESI.† |
This problem is effectively addressed by using density-based analysis:^{43–45} natural orbitals and their occupations provide a clear and unambiguous picture of the essential electronic structure of open-shell systems. Natural orbitals, which are eigen-functions of the one-particle density matrix, afford most compact representation of the wave function. The respective eigen-values can be interpreted as occupations. Thus, by computing natural orbitals for a given one-particle density matrix one can obtain frontier MOs representing, for example, the unpaired electrons in diradicals. The shapes of these orbitals allow one to asses the extent of localization of the unpaired electrons and their through-space and through bond interactions with each other and with other moieties. We note that when two frontier orbitals are exactly singly occupied (such as in perfect diradicals), their choice is not unique and any linear combination provides a legitimate representation (for example, in H_{2} with a completely broken bond, one can consider either a pair of localized atomic-like orbitals or a delocalized bonding–antibonding pair).
By using natural occupations, one can define and compute the number of effectively unpaired electrons. In this work, we use two indices, n_{u} and n_{u,nl}, proposed by Head-Gordon:^{43}
(7) |
(8) |
In both equations, the sum runs over all natural orbitals and the contributions of the doubly occupied and unoccupied orbitals are exactly zero. While both formulas give identical answers for the limiting cases (such as a two-electron triplet state), n_{u,nl} delivers more physically meaningful and consistent results for more complex wave functions. By comparing occupations of the frontier natural orbitals and n_{u,nl} one can assess how well a two-electrons-in-two-orbitals picture represents the real wave function.
We recently applied this technique to prototypical di- and triradicals.^{11} Here, we extend this approach^{11} to a set of 8 binuclear copper diradicals. We quantify the degree of radical character using Head-Gordon's indices,^{43}eqn (7) and (8), and visualize natural frontier orbitals to characterize the localization and the interactions between the unpaired electrons. Visual inspection of natural orbitals, paired with the analysis of the respective natural occupations and the computed number of effectively unpaired electrons, obtained from a variety of SF methods, provides a robust way of validating the applicability of the HDvV model (and the Landé interval rule) for describing this family of MM candidates.
Experimental geometries are used in all EOM-SF-CCSD and SF-TDDFT calculations, unless otherwise specified. To test possible effects of uncertainty of the structure, we also considered optimized structures of the high-spin triplet states of BISDOW, PATFIA (without the ferrocene group), and CITLAT. These optimizations were performed using the ωB97X-D^{78} functional and all-electron cc-pVTZ basis. The results for optimized structures are presented in ESI.†
In the tables below, we report singlet–triplet gaps, ΔE_{ST}, defined as:
ΔE_{ST} = E_{S} − E_{T}. | (9) |
We compared the following DFT functionals in the SF-TDDFT calculations of ΔE_{ST}: LDA (with Slater exchange and VWN correlation), several members of the Becke-exchange/LYP correlation family: BLYP,^{79,80} B3LYP,^{81} and B5050LYP (50% Hartree–Fock + 8% Slater + 42% Becke for exchange and 19% VWN + 81% LYP for correlation).^{35} Of the P86 correlation (with Becke exchange) and PW91 families, we chose the BP86, B3P86,^{82,83} PW91, and B3PW91 functionals.^{84–87} From the PBE family, we selected the PBE, PBE0 (75% PBE and 25% Hartree–Fock exchange, 100% PBE correlation), PBE50 (50% PBE and 50% Hartree–Fock exchange and 100% PBE correlation), and ωPBEh (80% PBE, 20% short-range Hartree–Fock exchange and 100% long-range Hartree–Fock exchange, PBE correlation) functionals.^{88,89} Of the Minnesota family of functionals, we chose M06 (hybrid with 27% Hartree–Fock exchange),^{90} M06-L (meta-GGA),^{90} GAM (GGA),^{91} MN15-L (meta-GGA),^{92} and MN15 (hybrid with 44% Hartree–Fock exchange and MN15 correlation).^{93} All of these functionals (with the exception of GAM, MN15-L and MN15, which are recent additions to the Minnesota family) have been extensively benchmarked with SF-TDDFT for organic polyradicals, wherein hybrid functionals such as PBE50, B5050LYP, and PBE0 were shown to produce relative state energies approaching chemical accuracy.^{34,35} In this paper, we present a similar error analysis of collinear and non-collinear SF-TDDFT with the above functionals in calculating energy gaps of binuclear copper complexes.
We also present EOM-SF-CCSD/cc-pVDZ energy splittings for selected complexes. The core electrons were frozen in all EOM calculations. For all SF-TDDFT and EOM-SF-CCSD calculations, we report only electronic energy separations between the singlet and M_{s} = 0 triplet states (ΔE_{ST}).^{60} Because of the similarity of the electronic structure of the singlet and triplet states in the case of very weakly interacting electrons, we expect that both the structures and vibrational frequencies of the two states are very close.
We performed density and wave function analysis using the libwfa module^{44} contained in the Q-Chem electronic structure package.^{94,95} We analyzed singlet and triplet states of all eight benchmark complexes at the PBE50/cc-pVDZ level of theory, and six of the eight complexes (BISDOW, CUAQAC02, Cu_{2}Cl_{6}^{−}, XAMBUI, YAFZOU, and CITLAT) at the EOM-SF-CCSD/cc-pVDZ level of theory. PBE50, B97, and LDA functionals were used with the non-collinear SF-TDDFT kernel^{35,96,97} in the analysis of CUAQAC02. Collinear SF-TDDFT with the B5050LYP functional, which was recommended in the original SF-TDDFT paper,^{34} was also used in the density analysis of CUAQAC02. We also present frontier natural orbitals of the full experimental structure of the PATFIA complex, which includes a ferrocene group, at the collinear SF-TDDFT/cc-pVDZ level with the B5050LYP functional. We compare these results to those obtained from the simplified structure without the ferrocene group.
We compared the performance of Dunning's cc-pVDZ and cc-pVTZ all-electron basis sets for all eight complexes and several functionals. Table S1 in the ESI† provides an atom-by-atom description of the effective core potentials (ECPs) and matching basis functions used in our study. Included are LANL2DZ (a non-relativistic ECP calibrated with Hartree–Fock total energies),^{98} SRSC (a relativistic ECP fitted with all electron-eigenvalues and charges),^{99} and CRENBL (a relativistic ECP fitted with Hartree–Fock valence orbital energies).^{100,101} Also included is the ECP10MDF pseudopotential,^{102,103} a relativistic ECP designed to reproduce valence energy spectra. ECP10MDF has been shown to perform well in the EOM-CCSD analysis of small copper compounds.^{104,105}
We performed all calculations with the Q-Chem electronic structure package.^{94,95} Molecular orbitals were rendered using IQmol^{106} and natural orbitals were rendered using Jmol.
Hybrid functionals—in particular, LRC-ωPBEh, PBE0, PBE50, and B5050LYP—outperform LDA and GGA functionals and approach the accuracy of EOM-SF-CCSD. As will be shown in subsequent sections, errors against the experiment are non-uniform among antiferromagnetic (AF) and ferromagnetic (F) complexes and the functionals and kernels that yield the closest agreement with experiment and EOM-SF-CCSD vary depending on the sign of J (i.e., whether the complex exhibits a singlet or triplet ground state). Although the lowest MAE is observed for the MN15 functional, energy differences between the high-spin triplet reference and the M_{s} = 0 component of the triplet state are often large (on the order of 2 eV, in contrast with typical values of 0.1–0.2 observed for other functionals, see, for example, Table S2 in ESI†) and spin-contamination of target states high, ranging between 0.01 and 0.33, whereas for other functionals such as B5050LYP and PBE50, spin-contamination is typically within 0.01–0.02. M06 shows larger errors relative to MN15. Compared to B5050LYP, M06 MAE is almost three times larger, likely because of a too-small fraction of the exact exchange (27%).
The non-collinear kernel outperforms the collinear SF-TDDFT kernel for all functionals, with the exception of hybrids PBE50 and B5050LYP, where the collinear kernel shows modest improvement. For GGA functionals, use of the collinear kernel is tantamount to Koopmans' estimate of the gaps using Kohn–Sham orbital energies, as the spin-flipped configurations are not coupled in the absence of the Hartree–Fock exchange.^{34} Consequently, energy errors are high for these functionals, relative to results obtained with the non-collinear kernel, and spin contamination of target states significant (see ESI†).
The scatter plot in Fig. 5 illustrates the singlet–triplet gaps of all eight complexes computed by several hybrid functionals and EOM-SF-CCSD. An idealized line showing perfect agreement between experimental values of J and theoretical singlet–triplet gaps is also shown. All five selected functionals capture the qualitative trend and reproduce the sign of the exchange-coupling. The errors are relatively small for complexes with near-zero experimental exchange-coupling and for complexes with a triplet ground state. We observe better agreement with experiment when hybrid functionals with 50% HF exchange are applied to molecules with a triplet ground state. For most functionals, the discrepancy between theoretical and experimental values are quite large for the two complexes with singlet ground states, BISDOW and CUAQAC02. Table 2, in which EOM-SF-CCSD/cc-pVDZ results for exchange-coupling are presented, shows the same trend. Again, we note good agreement with experiment for complexes that exhibit near-zero or ferromagnetic coupling, while substantial discrepancy is observed for anti-ferromagnetic complexes. It will be shown in Section 4.6 that this discrepancy arises not from electronic structure method failure, but the limitations of the HDvV model in the case of the singlet ground states of more strongly AF complexes.
Complex | EOM-CCSD ΔE_{ST} | 〈Ŝ^{2}〉_{ref}^{a} | 〈Ŝ^{2}〉_{t}^{a} | 〈Ŝ^{2}〉_{s}^{a} | Exp. J |
---|---|---|---|---|---|
a 〈S^{2}〉_{ref} denotes 〈S^{2}〉 of the M_{s} = 1 reference determinant. 〈S^{2}〉_{t}, and 〈S^{2}〉_{s} denote the M_{s} = 0 target triplet and singlet states, respectively. b Cholesky decomposition with a tolerance of 1 × 10^{−2}. Virtual orbitals frozen above 4.5 hartree (total of 14). c Cholesky decomposition with a tolerance of 1 × 10^{−3}. d EOM-SF-MP2 approximation. Virtual orbitals above 4 hartree (total of 20). | |||||
BISDOW | −208^{b} | 2.115 | 2.013 | 0.016 | −382 |
CUAQAC02 | −180^{c} | 2.006 | 2.000 | 0.005 | −286 |
Cu_{2}Cl_{6}^{2−} | −16^{c} | 2.009 | 1.997 | 0.009 | 0, −40 |
XAMBUI | 0.0 | 2.006 | 2.003 | 0.001 | 2 |
YAFZOU | 86^{c} | 2.022 | — | — | 111 |
CITLAT | 59^{d} | 2.008 | 1.980 | 0.028 | 113 |
Basis | Complex | |||||||
---|---|---|---|---|---|---|---|---|
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |
cc-pVDZ | −444 | −319 | −24 | −35 | 39 | 0 | 198 | 231 |
ECP10MDF/cc-pVDZ | −416 | −296 | −24 | −44 | 40 | 0 | 184 | 213 |
LANL2DZ | −510 | −532 | −20 | 38 | −2 | −1 | 244 | 263 |
Exp. J | −382 | −286 | −19 | −11 | (0, −40) | −2 | 111 | 113 |
Basis | Complex | |||||||
---|---|---|---|---|---|---|---|---|
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |
cc-pVDZ | −523 | −354 | −27 | −35 | −23 | 2 | 201 | 232 |
cc-pVTZ | −506 | −362 | −26 | −33 | 35 | 1 | 197 | 227 |
ECP10MDF/cc-pVDZ | −486 | −326 | −27 | −40 | 23 | −1 | 186 | 213 |
ECP10MDF/cc-pVTZ | −456 | −326 | −25 | −40 | 40 | 1 | 84 | 204 |
LANL2DZ | −615 | −603 | −21 | 39 | 0 | 0 | 250 | 267 |
SRSC | −471 | −449 | −19 | 21 | −10 | 1 | 191 | 208 |
CRENBL | −474 | −459 | −19 | 31 | −75 | 0 | 202 | 223 |
Exp. J | −382 | −286 | −19 | −11 | (0, −40) | −2 | 111 | 113 |
Basis | Complex | |||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |||||||||
NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | |
a Non-collinear TDDFT kernel. b Collinear TDDFT kernel. | ||||||||||||||||
cc-pVDZ | −173 | −219 | −121 | −154 | −9 | −11 | 0 | 60 | 55 | −3 | 0 | 0 | 93 | 77 | 120 | 78 |
cc-pVTZ | −158 | −211 | −119 | −153 | −7 | −11 | 1 | −57 | 63 | 7 | 1 | −1 | 90 | 60 | 118 | 77 |
ECP10MDF/cc-pVDZ | −158 | −213 | −115 | −148 | 8 | 11 | 15 | −60 | −52 | −3 | −1 | 1 | 88 | 57 | 112 | 75 |
ECP10MDF/cc-pVTZ | −148 | −199 | −111 | −143 | −6 | −11 | −8 | −56 | 61 | 15 | 0 | 0 | 84 | 56 | 110 | 74 |
LANL2DZ | −203 | −262 | −208 | −269 | −9 | −10 | 15 | −60 | −18 | −85 | 0 | 0 | 141 | 80 | 114 | 86 |
SRSC | −160 | −208 | −156 | −198 | −7 | −9 | −15 | −46 | 7 | — | 0 | 0 | 79 | 60 | 99 | 73 |
CRENBL | −161 | −211 | −158 | −201 | −7 | −15 | −15 | −50 | −35 | −90 | 0 | 0 | 92 | 60 | 99 | 74 |
Exp. J | −382 | −286 | −19 | −11 | 0, −40 | −2 | 111 | 113 |
Basis | Complex | |||||||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | |||||||||
NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | NC^{a} | C^{b} | |
a Non-collinear TDDFT kernel. b Collinear TDDFT kernel. | ||||||||||||||||
cc-pVDZ | −155 | −212 | −114 | −149 | −8 | −11 | 2 | −57 | 45 | −7 | −1 | 0 | 84 | 61 | 109 | 77 |
cc-pVTZ | −150 | −205 | −115 | −149 | −6 | −10 | 2 | −56 | 53 | 2 | −1 | 0 | 82 | 60 | 108 | 77 |
ECP10MDF/cc-pVDZ | −142 | −205 | −103 | −140 | 7 | 11 | 9 | −57 | 44 | 0 | 1 | 1 | 75 | 56 | 98 | 74 |
ECP10MDF/cc-pVTZ | −133 | −191 | −99 | −136 | −2 | −10 | −8 | −56 | 52 | 12 | 1 | 0 | 73 | 53 | 92 | 73 |
LANL2DZ | −204 | −255 | −211 | −261 | −7 | −10 | 7 | −31 | −27 | −83 | 0 | −1 | 97 | 79 | 110 | 85 |
SRSC | −162 | −199 | −158 | −192 | −7 | −9 | 22 | −38 | −7 | — | 0 | 0 | 72 | 58 | 90 | 72 |
CRENBL | −165 | −205 | −162 | −199 | −9 | −9 | 6 | −42 | −48 | −95 | −2 | 0 | 79 | 60 | 94 | 73 |
Exp. J | −382 | −286 | −19 | −11 | 0, −40 | −2 | 111 | 113 |
For BISDOW, the most anti-ferromagnetic complex in our study, the modest Hartree–Fock exchange of the PBE0 and LRC-ωPBEh functionals paired with the non-collinear kernel appears to yield results in closest agreement with the experimental values of J (but in disagreement with EOM-SF-CCSD). However, even with the PBE0 and LRC-ωPBEh functionals, the discrepancy is always larger than for ferromagnetic complexes. Exchange-coupling between copper centers in CITLAT, the most ferromagnetic complex in our study, is best captured by functionals with a higher percentage of Hartree–Fock exchange, such as in PBE50 and B5050LYP. The sign of the weakly AF complex, PATFIA, is captured by both the PBE0 and LRC-ωPBEh functionals, regardless of kernel choice, while the magnitude is better represented by the PBE50 or B5050LYP functionals and the non-collinear kernel.
Truhlar and co-workers^{46} carried out a detailed benchmark study in which they compared collinear SF-TDDFT with BS-DFT (and with REKS) using several functionals with varying amount of Hartree–Fock exchange (X): B3LYP (X = 20), M06 (X = 27), B3LYP40 (X = 40), B1LYP40 (X = 40), B1PW40 (X = 40), BMK (X = 42), MPW1K (X = 42.8), B3LYP54 (X = 54), M06-2X (X = 54), M06-HF (X = 100). Of this set, B3LYP54 is very similar to B5050LYP. Their test set comprised 12 binuclear copper MMs, including all 8 MMs investigated here. This study reported that SF-TDDFT approach performs systematically better than the spin-projected weighted-average BS strategy, although the optimal percentage of Hartree–Fock exchange was found to be slightly different for the two approaches (about 40% for SF-TDDFT and 35% for spin-projected weighted-average BS-DFT). In agreement with our findings, they pointed out that the SF-TDDFT exchange-couplings show a much stronger dependence on the percentage of Hartree–Fock exchange than on the other details of the exchange functionals or the nature of the correlation functionals.
Table 7 summarizes exchange-coupling constants computed for the two types of geometries, four hybrid functionals and three ECP/basis combinations. Overall, we find that singlet–triplet gaps computed at a given level of theory changes very little with subtle changes in nuclear geometry, indicating that possible uncertainties in the experimentally derived geometries cannot account for the discrepancies between the computed and experimental exchange-couplings. Only for PATFIA, we observe a slightly better agreement between theory and experiment when using the optimized structure. Overall, these results justify using crystal structures for exchange-couplings calculations.
Method | Basis | Complex 1 | Complex 4 | Complex 8 | |||
---|---|---|---|---|---|---|---|
X-Ray | Optimized | X-Ray | Optimized | X-Ray | Optimized | ||
a SCF calculation converged to a wrong reference state. | |||||||
PBE0 | cc-pVDZ | −523 | −524 | −35 | −57 | 109 | DRS^{a} |
ECP10MDF/ccpvdz | −486 | −459 | −40 | −61 | 213 | 233 | |
LANL2DZ | −615 | −571 | 39 | −12 | 267 | 313 | |
LRC-ωPBEh | cc-pVDZ | −444 | −418 | −35 | −59 | 231 | 259 |
ECP10MDF/ccpvdz | −416 | −394 | −44 | −64 | 213 | 233 | |
LANL2DZ | −510 | −473 | 38 | −15 | 263 | 306 | |
PBE50 | cc-pVDZ | −173 | −161 | 0 | −24 | 120 | 113 |
ECP10MDF/ccpvdz | −158 | −148 | 15 | −31 | 112 | 104 | |
LANL2DZ | −203 | −180 | 15 | −16 | 114 | 134 | |
B5050LYP | cc-pVDZ | −155 | −152 | 2 | −25 | 109 | 102 |
ECP10MDF/ccpvdz | −142 | −133 | 9 | −32 | 98 | 90 | |
LANL2DZ | −204 | −165 | 7 | −19 | 110 | 115 | |
Exp. J | −382 | −11 | 113 |
Molecule | ΔE_{ST} | n _{u} | n _{u,nl} | 〈Ŝ^{2}〉 | |||
---|---|---|---|---|---|---|---|
Singlet | Triplet | Singlet | Triplet | Singlet | Triplet | ||
BISDOW | −173 | 1.81 | 2.01 | 1.96 | 2.00 | 0.01 | 2.02 |
CUAQAC02 | −121 | 1.84 | 2.01 | 1.97 | 2.00 | 0.01 | 2.01 |
CAVXUS | −9 | 1.96 | 2.01 | 2.00 | 2.00 | 0.07 | 1.97 |
PATFIA | 0 | 1.86 | 1.98 | 1.98 | 2.00 | 0.36 | 1.68 |
Cu_{2}Cl_{6}^{2−} | 55 | 1.85 | 2.02 | 1.97 | 2.00 | 0.02 | 2.03 |
XAMBUI | 0 | 2.00 | 2.01 | 2.00 | 2.00 | 0.01 | 2.02 |
YAFZOU | 93 | 1.97 | 2.01 | 2.00 | 2.00 | 0.01 | 2.02 |
CITLAT | 120 | 1.96 | 2.01 | 2.00 | 2.00 | 0.01 | 2.02 |
As expected, for all triplet states the number of effectively unpaired electrons as given by n_{u,nl} is exactly 2 (n_{u} gives slightly larger values, similarly to our previous study or organic diradicals^{11}). We observe that for the BISDOW and PATFIA complexes—for which the recommended SF-TDDFT methods and EOM-SF-CCSD consistently disagree with experiment—the number of unpaired electrons (i.e., values of n_{u} and n_{u,nl}) in the ground-state singlets are less than 2, indicative of weaker diradical character and mixing in of the ionic configurations. These effects are outside of the domain of the HDvV Hamiltonian. The weakly AF complexes Cu_{2}Cl_{6}^{2−} and PATFIA, for which top-performing SF-TDDFT methods disagree with experiment in both magnitude and sign, also exhibit weaker diradical character. CAVXUS (a weakly AF complex), XAMBUI (a complex that exhibits zero exchange-coupling between radical centers), and the ferromagnetic YAFZOU and CITLAT complexes all have singlet states with the n_{u,nl} values at the ideal diradical value of 2. These are also the four complexes for which top-performing DFT functionals and EOM-SF-CCSD agree well with the experimental exchange-coupling values.
Spin-contamination appears to be minor: for most complexes the deviation from the exact 〈Ŝ^{2}〉 values are 0.01–0.03. In CAVCUS, the deviation is slightly larger (0.07 for the singlet state). The largest spin-contamination is observed in PATFIA where the 〈Ŝ^{2}〉 of the singlet state is 0.36. Spin-balance is also evidenced by the natural orbitals occupations: |n_{α} − n_{β}| values of exactly zero for BISDOW, CUAQAC02, Cu_{2}Cl_{6}^{2−}, XAMBUI, and CITLAT.
The frontier natural orbitals depicted in Fig. 7 are all of the d_{xy} or d_{yz} type, exhibiting weak interaction with p-orbitals of nearest ligand atoms. Mixing of open- and closed-shell configurations is symmetry-allowed for CAVXUS, PATFIA, and YAFZOU. When spin-traced occupation numbers are exactly equal, both localized and delocalized frontier orbitals provide a valid representation of the density. Thus, the apparent localization of the frontier orbitals in CAVXUS, PATFIA, and the M_{s} = 0 triplet state of CITLAT does not result from the spatial or spin symmetry breaking. We note that populations associated with frontier natural orbitals are consistent with the computed number of the unpaired electrons, meaning that electronic structure of these compounds maps well into two-electrons-in-two-orbitals problem.
The d_{yz} bonding/anti-bonding pair of frontier orbitals observed in CUAQAC02 is unique among the other copper diradicals considered in our study. The δδ* natural orbital orientation suggests a level of electron–electron interaction that might preclude the use of the HDvV Hamiltonian for modeling exchange-coupling, and lead to the observed imperfect diradical character of the singlet ground state.
Using the interesting case of CUAQAC02, Table 9 and Fig. 8 illustrate the agreement between the SF-TTDFT/PBE50, SF-TDDFT/B5050LYP, and EOM-SF-CCSD, not just for the relative energetics of the states of interest, but their associated densities. While LDA and B97 fail to accurately describe state energetics and the degree of diradical character associated with the singlet (and, in the case of B97, the triplet) states, the shapes of the frontier natural orbitals predicted by LDA are consistent with EOM-SF-CCSD, PBE50 and B5050LYP. The main difference between LDA and other methods is in natural occupations (and, consequently, the number of effectively unpaired electrons): LDA's singlet state has much stronger closed-shell character. Shown in the ESI† (Fig. S3), B97 singly occupied natural orbitals exhibit a σσ-type interaction with equatorial nearest-neighbor oxygen atoms, meaning that the lowest energy triplet state returned by this functional is altogether inconsistent with the other four methods. Using standard troubleshooting tools, we were not able to find an SCF solution with this functional that corresponds to the same bonding pattern as predicted by other functionals. Using a wrong reference state leads to large errors in the singlet triplet gaps and large spin-contamination (Table 9). This case illustrates the utility of using natural orbitals for detecting problematic situations with SCF convergence (as one can see from Fig. 2, the inspection of canonical orbitals is not very instructive for these highly non-Koopmans systems).
Method | ΔE_{ST} | n _{u} | n _{u,nl} | 〈Ŝ^{2}〉 | |||
---|---|---|---|---|---|---|---|
Singlet | Triplet | Singlet | Triplet | Singlet | Triplet | ||
NC-PBE50 | −123 | 1.84 | 2.01 | 1.97 | 2.00 | 0.01 | 2.01 |
COL-B5050LYP | −148 | 1.81 | 2.01 | 1.96 | 2.00 | 0.01 | 2.01 |
NC-LDA | −1420 | 0.57 | 2.00 | 0.47 | 2.00 | 0.01 | 2.01 |
NC-B97 | 71 | 2.75 | 2.82 | 2.59 | 2.69 | 1.07 | 1.77 |
EOM-SF-CCSD | −171 | 1.79 | 1.91 | 1.96 | 2.00 | 0.00 | 2.00 |
Finally, we consider the case of PATFIA when the full experimental structure that includes the ferrocene group is used in the SF-TDDFT calculation. Results at the B5050LYP/cc-pVDZ level of theory (with the collinear TDDFT kernel) are summarized for simplified PATFIA and PATFIA + Fe(C_{5}H_{5})_{2} in Table 10 and Fig. 9. The observed localization of unpaired electron density on the copper atoms in PATFIA + (C_{5}H_{5})_{2}, and the consistency in the computed values of the singlet–triplet gap, n_{u} and n_{u,nl}, suggest that the ferrocene group indeed does not interact with unpaired electrons on the copper centers, and that the ubiquitous but thus-far unsubstantiated use of the simplified structure in theoretical calculations of exchange-coupling is valid.
Molecule | ΔE_{ST} | n _{u} | n _{u,nl} | 〈Ŝ^{2}〉 | |||
---|---|---|---|---|---|---|---|
Singlet | Triplet | Singlet | Triplet | Singlet | Triplet | ||
PATFIA | −57 | 1.80 | 2.01 | 1.95 | 2.00 | 0.02 | 2.02 |
PATFIA + Fe(C_{5}H_{5})_{2} | −60 | 1.80 | 2.01 | 1.96 | 2.00 | 0.05 | 1.98 |
Small variations in nuclear geometries—specifically, a slight increase or decrease in the internuclear separation between the Cu atoms that host the unpaired electrons—does not impact the computed value of ΔE_{ST}, validating the use of experimentally determined geometries (or slightly modified symmetrized geometries) in calculations of exchange-coupling. Furthermore, the presence of the ferrocene group that is often excluded in ab initio studies of exchange-coupling in the PATFIA complex did not alter the value of ΔE_{ST} or the underlying density, justifying the use of the simplified structure.
We employ density-based analysis to gain insight into the interaction between the unpaired electrons in copper diradicals and, in so doing, discover that only five of the eight binuclear copper complexes exhibit perfect or near-perfect diradical character in their lowest singlet and triplet states. This finding suggests that while the HDvV Hamiltonian might be an appropriate model for exchange-coupling in CAVXUS, PATFIA, Cu_{2}Cl_{6}^{2−}, XAMBUI, and CITLAT, quality of results should be expected to suffer for complexes like BISDOW and CUAQAC02, wherein ionic configurations are mixed into the singlet ground states. For these species, the experimentally derived values of J are expected to be systematically overestimated. Thus, the observed discrepancies between theoretical and experimental values for AF molecules arise not due to the limitations of theoretical methods, but due to the limitations of the HDvV model. Spin-contamination appears to be minor in this class of compounds.
We have illustrated that the PBE50 or B5050LYP functionals (with a modest choice of basis set and ECP) can be paired with the SF-TDDFT approach to paint a robust picture of the energetics of MM candidates, with a level of accuracy comparable to EOM-SF-CCSD. When combined with the density-based analysis of the relevant spin states, one can determine whether or not the unpaired electrons in a given MM obey the underlying physics of the HDvV Hamiltonian. Screening of MM candidates based solely on degree of radical character observed in the low-lying spin states is made possible. Recent work by Mayhall and Head-Gordon extends the formalism of SF to computation of exchange-coupling constants in complexes with greater than two radical sites and/or unpaired electrons.^{18,19}
Footnote |
† Electronic supplementary information (ESI) available: Relevant Cartesian geometries, details of ECP, additional results for collinear kernels. See DOI: 10.1039/c7cp07356a |
This journal is © the Owner Societies 2018 |