Singlet–triplet energy gaps and the degree of diradical character in binuclear copper molecular magnets characterized by spin-flip density functional theory

Natalie Orms and Anna I. Krylov *
Department of Chemistry, University of Southern California, Los Angeles, California 90089-0482, USA. E-mail:

Received 31st October 2017 , Accepted 22nd January 2018

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.

I. Introduction

Molecular magnets (MMs) are organic polyradicals that have a high-spin ground state in the absence of an applied magnetic field.1,2 They can be used as building blocks to create novel, light-weight, and tunable magnetic materials, an alternative to conventional dense metallic or metal oxide magnets. MMs are also of interest in the context of quantum computing and quantum information storage.3–5 Since the 1990's much attention has been given to the design and study of organometallic complexes containing d- and f-block elements—including terbium, manganese, chromium, nickel, and copper—whose electronic structure (and consequent magnetic properties) make them suitable candidates for use in magnetic materials.6 Suitable candidates for magnetic applications should have one or more unpaired, weakly interacting electrons. Here, we limit ourselves to systems containing only two unpaired electrons (diradicals).7–11 Electronic configurations that can be generated for two electrons in two orbitals are shown in Fig. 1.
image file: c7cp07356a-f1.tif
Fig. 1 Wave functions of diradicals that are eigenfunctions of S2 (only configurations with positive spin projections are shown). Wave function (i) corresponds to the high-spin Ms = 1 triplet state. Wave functions (ii)–(v) correspond to the low-spin states: Ms = 0 singlets and triplets. Note that all Ms = 0 configurations can be formally generated by a spin-flipping excitation of one electron from the high-spin Ms = 1 configuration. The character of the wave functions (i.e., covalent versus ionic) depends on the nature of molecular orbitals ϕ1 and ϕ2 and the value of λ.

The triplet state of a diradical can be represented by either a single Slater determinant (high-spin Ms = 1 configuration (i)) or by a linear configuration of two Ms = 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)
with little-to-no contributions from ionic configurations
[ϕ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 Ms = 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 Ms = 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) methods20–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) model14–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 model14–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.

II. Theoretical methods

A. Diradicals and spin-flip approach

As described above, the electronic states of diradicals,7,8,10,59 such as bicopper MMs, arise from two electrons distributed in two nearly degenerate MOs. For same-spin electrons, there is only one possible arrangement: Ms ± 1. For the Ms = 0 states, more configurations can be generated, as illustrated in Fig. 1. All configurations depicted in Fig. 1 can contribute to diradical spin-states of interest, which makes choosing an electronic structure method that treats all states on an equal footing, as SF does, of the utmost importance. Because in the SF approach there is no assumptions of the relative importance of configurations from Fig. 1, these calculations allow one to unambiguously determine the state character, e.g., to quantify the extent of the interaction between the unpaired electrons and the weights of ionic configurations in the respective wave functions. Thus, SF methods are capable of distinguishing whether a given MM candidate falls within the regime of spatially separated, weakly-interacting spin moments, as assumed within the HDvV model.

To describe diradicals within the SF framework, one begins with the high-spin Ms = 1 triplet state, which is well represented by a single Slater determinant. The low-spin states (singlets and the Ms = 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

image file: c7cp07356a-t1.tif(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.

B. The Heisenberg–Dirac–van-Vleck Hamiltonian for localized, weakly interacting electrons

The key assumption used to construct phenomenological spin Hamiltonians is that the unpaired electrons are spatially separated and weakly interacting. The resulting magnetic Hamiltonian is:
image file: c7cp07356a-t2.tif(4)
where the ŜA are local electronic spin operators on nuclei A, the ÎA are nuclear spin moments, and D, A, and J represent the zero-field splitting, hyperfine, and exchange-coupling interaction, respectively.2,12,13 Magnetic parameters D, A, and J are extracted from experiment, typically by fitting electron paramagnetic resonance, magnetic susceptibility, or other macroscopic properties to a model equation.16,65–68 Solving this Hamiltonian allows one to predict magnetic behavior of the system in the presence or absence of the magnetic field.

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

image file: c7cp07356a-t3.tif(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) = −JAB·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., ΔEST in the case of diradicals).

C. Quantifying radical character and bonding patterns by density-based analysis

Due to the complexity of many-body wave functions and arbitrariness in orbital choices, the interpretation of realistic wave functions (and often even of Kohn–Sham states) in terms simple two-electron-in-two orbitals picture (as in Fig. 1) is not straightforward. In particular, in large open-shell systems canonical Hartree–Fock of Kohn–Sham orbitals often provide a poor representation of the frontier MOs even for relatively simple high-spin states.11Fig. 2 illustrates the difference between canonical Kohn–Sham orbitals and spin-density difference for copper diradicals: while the excess spin density exhibits an expected pattern (unpaired electrons are localized on d-orbitals of copper), the two formally singly occupied canonical MOs are delocalized over the ligands. In this case, even assigning the state characters and verifying the correctness of the solution of the self-consistent field (SCF) procedure are problematic.
image file: c7cp07356a-f2.tif
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.

image file: c7cp07356a-f3.tif
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 H2 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, nu and nu,nl, proposed by Head-Gordon:43

image file: c7cp07356a-t4.tif(7)
image file: c7cp07356a-t5.tif(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), nu,nl delivers more physically meaningful and consistent results for more complex wave functions. By comparing occupations of the frontier natural orbitals and nu,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 approach11 to a set of 8 binuclear copper diradicals. We quantify the degree of radical character using Head-Gordon's indices,43eqn (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.

III. Computational details

Fig. 3 and Fig. S1 (ESI) show the experimental structures used in calculations of singlet–triplet gaps of benchmark systems. Table 1 provides their associated experimental exchange coupling constants. Counterions were removed from all structures. The ferrocene group in the experimental PATFIA complex was also removed (as in previous studies46,58), unless otherwise specified.
Table 1 Experimental exchange-coupling constants for eight binuclear copper diradicals shown in Fig. 3
Complex J (cm−1) Ref.
BISDOW −382 50
CUAQAC02 −286 51
CAVXUS −19 52 and 53
PATFIA −11 54
Cu2Cl62− 0, −40 55
YAFZOU 111 57
CITLAT 113 49

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-D78 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, ΔEST, defined as:

within the HDvV model, ΔEST = J, by virtue of eqn (6).

We compared the following DFT functionals in the SF-TDDFT calculations of ΔEST: 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 Ms = 0 triplet states (ΔEST).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 module44 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, Cu2Cl6, 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 kernel35,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 IQmol106 and natural orbitals were rendered using Jmol.

IV. Results and discussion

A. The comparison of singlet–triplet gaps computed by different methods

We begin by comparing computed singlet–triplet gaps against experimentally derived exchange-coupling constants. Mean absolute errors (MAEs) for each DFT functional and EOM-SF-CCSD are presented in Fig. 4. Tabulated mean errors (ME), mean absolute errors (ΔMAE), and standard deviations of the error (ΔSTD) are provided in the ESI for the non-collinear kernel, the PBE0, PBE50, and B5050LYP functionals, and various all-electron basis sets and ECPs.
image file: c7cp07356a-f4.tif
Fig. 4 Mean absolute error (MAE) in the singlet–triplet gap for the eight copper benchmark systems. Error relative to experimental values of the exchange-coupling constant, J, are presented for 18 density functionals and EOM-SF-CCSD. The all-electron cc-pVDZ basis set was used for all atoms.

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 Ms = 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.

image file: c7cp07356a-f5.tif
Fig. 5 Theoretical singlet–triplet gaps computed with the PBE0, ω-PBEh, B3LYP, and B5050LYP functionals (y-axis) versus experimentally derived values of the exchange-coupling parameter (x-axis). The inset in the bottom right zooms into the region spanning the range of −50 to 10 cm−1. Dashed diagonal line marks perfect agreement between theoretical and experimental values. The all-electron cc-pVDZ basis was applied to all atoms. The collinear SF-TDDFT kernel was used with the B5050LYP functional, and non-collinear SF-TDDFT kernel with all other functionals.
Table 2 EOM-SF-CCSD/cc-pVDZ singlet–triplet gaps (computed using the X-ray structures) and experimental exchange-coupling constants, in cm−1
Complex EOM-CCSD ΔEST Ŝ2refa Ŝ2ta Ŝ2sa Exp. J
a S2ref denotes 〈S2〉 of the Ms = 1 reference determinant. 〈S2t, and 〈S2s denote the Ms = 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 −208b 2.115 2.013 0.016 −382
CUAQAC02 −180c 2.006 2.000 0.005 −286
Cu2Cl62− −16c 2.009 1.997 0.009 0, −40
XAMBUI 0.0 2.006 2.003 0.001 2
YAFZOU 86c 2.022 111
CITLAT 59d 2.008 1.980 0.028 113

B. SF-TDDFT: the effect of ECP and basis sets

Tables 3–6 provide singlet–triplet gaps calculated with four hybrid functionals, and seven ECP/basis set combinations. For the PBE0 and ωPBEh functionals, results obtained with the collinear kernel are tabulated in the ESI. The use of smaller basis sets and ECPs affords computational savings, which is attractive for large systems, however, careful benchmarking is required to determine the effect of additional approximations on the computed quantities. We find that regardless of functional, there is little difference in quality of results when a triple-zeta, rather than double-zeta, basis set is used, justifying the choice of a smaller DZ-quality basis for calculations of exchange-coupling in binuclear copper diradicals. Use of an ECP (regardless of ECP choice) appears to yield slightly better agreement with experiment for AF complexes, while little-to-no effect is observed for F complexes.
Table 3 Singlet–triplet gaps (cm−1) for complexes 1–8 calculated using SF-TDDFT with the LRC-ωPBEh method and the non-collinear TDDFT kernel
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

Table 4 Singlet–triplet gaps (cm−1) for complexes 1–8 calculated using SF-TDDFT with the PBE0 method and the non-collinear TDDFT kernel
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

Table 5 Singlet–triplet gaps (cm−1) for complexes 1–8 calculated using SF-TDDFT with the PBE50 method
Basis Complex
1 2 3 4 5 6 7 8
NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb
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

Table 6 Singlet–triplet gaps (cm−1) for complexes 1–8 calculated using SF-TDDFT with the B5050LYP method
Basis Complex
1 2 3 4 5 6 7 8
NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb NCa Cb
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

C. SF-TDDFT: collinear versus non-collinear kernel

Fig. 4, Tables 3–6 present error analysis and values of the singlet–triplet gap computed with the collinear SF-TDDFT kernel and using the geometries from crystal structures. Fig. 6 compares the singlet–triplet gap calculated for the BISDOW (AF), PATFIA (weakly AF), and CITLAT (F) complexes using the collinear and non-collinear SF-TDDFT kernels. The ECP10MDF/cc-pVDZ basis is used, where the ECP is applied only to copper atoms. Errors are non-uniform with respect to both functional and kernel choice, a result that is consistent with findings for other complexes not included in Fig. 6.
image file: c7cp07356a-f6.tif
Fig. 6 Singlet–triplet gaps for complexes 1–8 calculated using the PBE0 (A), LRC-ωPBEh (B), PBE50 (C), and B5050LYP (D) density functionals. The performance of the non-collinear and collinear TDDFT kernels is compared in (C) and (D). The ECP10MDF pseudopotential and matching cc-pVDZ-PP basis set were applied to Cu atoms. The cc-pVDZ basis was used for all other atoms.

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.

D. BS-DFT versus SF-DFT

Two studies compared the performance of BS-DFT and SF-DFT for the binuclear copper MMs.25,46 Ziegler and co-workers focused on the performance of SF constricted variational DFT, SF-CV-DFT; they reported excellent performance of SF-CV-DFT with the B5050LYP functional and extensively discussed limitations of the BS-DFT and REKS/ROKS approaches. They reported BS-DFT results obtained with several functionals, illustrating strong dependence of the exchange-coupling values on the functional. The errors of BS-DFT with B5050LYP functional were comparable to those of SF-CV-DFT.

Truhlar and co-workers46 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.

E. The effect of geometry on the computed gaps

To investigate the possible effect of geometry on the computed singlet–triplet gaps, we carried out geometry optimizations for three of the eight benchmark complexes: BISDOW, PATFIA, and CITLAT. The resulting structures are shown in Fig. S2 in ESI. The most important structural parameter is the distance between the radical sites (Cu atoms). We found that the internuclear Cu distance of the optimized triplet geometries of BISDOW and CITLAT changed very little, with an increase of 0.03 Å and a decrease of 0.02 Å, respectively. We observed a larger change in the optimized PATFIA complex (0.19 Å), possibly due to the absence of the ferrocene group in the simplified structure, a topic that is examined more thoroughly in the next section.

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.

Table 7 Singlet–triplet gaps (cm−1) for optimized geometries of complexes 1, 4 and 8, compared with gaps obtained using the experimentally determined molecular geometries
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 DRSa
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

F. Wavefunction analysis

Energies alone are not sufficient to assess the performance of electronic structure methods.11 To gain further insight into performance of SF-TDDFT, we now investigate the characters of the underlying electronic states and compare relevant quantities (such as the number of effectively unpaired electrons) computed by SF-TDDFT and EOM-SF-CCSD. Table 8 presents density-based analysis of the lowest singlet and triplet SF-TDDFT states of all eight complexes at the PBE50/cc-pVDZ level of theory. Fig. 7 shows the associated frontier natural orbitals.
Table 8 NC-PBE50/cc-pVDZ energy splittings (cm−1) and wave function properties of lowest singlet and Ms = 0 triplet states of all eight benchmark complexes
Molecule ΔEST 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
Cu2Cl62− 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

image file: c7cp07356a-f7.tif
Fig. 7 Frontier natural orbitals of lowest singlet and Ms = 0 triplet states of copper diradicals at the PBE50/cc-pVDZ level of theory. With the exception of YAFZOU, for other systems only α orbitals of the triplet states are shown, since spatial extent of relevant paired α and β natural orbitals—and the appearance of frontier natural orbitals associated with the lowest singlet and triplet states—does not differ. [n with combining macron] = |nα + nβ|, with |nαnβ| provided in parentheses. [n with combining macron]s and [n with combining macron]t correspond to [n with combining macron] values obtained from the occupancies of the singlet and triplet natural orbitals, respectively.

As expected, for all triplet states the number of effectively unpaired electrons as given by nu,nl is exactly 2 (nu gives slightly larger values, similarly to our previous study or organic diradicals11). 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 nu and nu,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 Cu2Cl62− 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 nu,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, Cu2Cl62−, XAMBUI, and CITLAT.

The frontier natural orbitals depicted in Fig. 7 are all of the dxy or dyz 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 Ms = 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 dyz 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).

Table 9 SF-TDDFT/cc-pVTZ and EOM-SF-CCSD/cc-pVDZ energy splittings (cm−1) and wave function properties of lowest singlet and Ms = 0 triplet spin states of CUAQAC02
Method ΔEST 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

image file: c7cp07356a-f8.tif
Fig. 8 Frontier α natural orbitals of lowest singlet and Ms = 0 triplet states of CUAQAC02 at the SF-TDDFT/cc-pVTZ and EOM-SF-CCSD/cc-pVDZ levels of theory. Natural orbital surfaces obtained with three different density functionals are compared. The collinear TDDFT kernel was used with the B5050LYP functional. [n with combining macron] = |nα + nβ|, with |nαnβ| provided in parentheses. [n with combining macron]s and [n with combining macron]t correspond to [n with combining macron] values obtained from the occupancies of the singlet and triplet natural orbitals, respectively.

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(C5H5)2 in Table 10 and Fig. 9. The observed localization of unpaired electron density on the copper atoms in PATFIA + (C5H5)2, and the consistency in the computed values of the singlet–triplet gap, nu and nu,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.

Table 10 B5050LYP/cc-pVDZ SF-TDDFT energy splittings (cm−1) and wave function properties of lowest singlet and Ms = 0 triplet spin states of simplified PATFIA and the full experimental structure, PATFIA + Fe(C5H5)2. The collinear TDDFT kernel was applied
Molecule ΔEST 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(C5H5)2 −60 1.80 2.01 1.96 2.00 0.05 1.98

image file: c7cp07356a-f9.tif
Fig. 9 Frontier α natural orbitals of lowest singlet and Ms = 0 triplet states of simplified PATFIA and the full experimental structure, PATFIA + Fe(C5H5)2. Natural orbitals were obtained at the B5050LYP/cc-pVDZ level of theory with the collinear TDDFT kernel applied. [n with combining macron] = |nα + nβ|, with |nαnβ| provided in parentheses. [n with combining macron]s and [n with combining macron]t correspond to [n with combining macron] values obtained from the occupancies of the singlet and triplet natural orbitals, respectively.

V. Conclusions

We report a benchmark study of eight binuclear copper diradicals, focusing on the singlet–triplet energy gaps (which are related to exchange-coupling constants within the HDvV model) and the associated wave functions. We carefully compare SF-TDDFT methods with selected functionals against a high-level wave function based method, EOM-SF-CCSD. In agreement with previous calibration studies of organic di- and triradicals,11,34,35 we find that hybrid functionals outperform LDA and GGA-type functionals. The non-collinear kernel vastly improves the quality of SF-TDDFT results for GGA and LDA functionals, with modest improvement for most hybrids. Performance of DZ and TZ-quality Dunning's basis sets are comparable, while the use of various ECPs has little effect on the computed value of ΔEST, suggesting that computational savings can be achieved without sacrificing quality of results when extending the approaches outlined in this work to larger binuclear (or trinuclear) copper complexes.

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 ΔEST, 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 ΔEST 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, Cu2Cl62−, 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

Conflicts of interest

The authors declare the following competing financial interest(s): A. I. K. is a member of the Board of Directors and a part-owner of Q-Chem, Inc.


This work is supported by the Department of Energy through the DE-FG02-05ER15685 grant (A. I. K.). A. I. K. thanks Prof. Joachim Sauer for stimulating discussions and his suggestions regarding the effects of ECPs and basis sets.


  1. P. M. Lahti, Magnetic properties of organic materials, ed. P. M. Lahti, Marcel Dekker, Inc, 1999 Search PubMed.
  2. M. R. Pederson and T. Baruah, in Handbook of Magnetism and Magnetic Materials, ed. S. Parkin, J. Wiley and Sons, London, 2007, p. 1 Search PubMed.
  3. M. N. Leuenberger and D. Loss, Quantum computing in molecular magnets, Nature, 2001, 410, 789 CrossRef PubMed.
  4. C. Schlegel, J. vanSlageren, M. Manoli, E. K. Brechin and M. Dressel, Direct observation of quantum coherence in single-molecule magnets, Phys. Rev. Lett., 2008, 101, 147203 CrossRef CAS PubMed.
  5. M. D. Jenkins, D. Zueco, O. Roubeau, G. Aromi, J. Majere and F. Luis, A scalable architecture for quantum computation with molecular nanomagnets, Dalton Trans., 2016, 45, 16682 RSC.
  6. R. A. Layfield, Organometallic single-molecule magnets, Organometallics, 2014, 33, 1084–1099 CrossRef CAS.
  7. L. Salem and C. Rowland, The electronic properties of diradicals, Angew. Chem., Int. Ed., 1972, 11, 92–111 CrossRef CAS.
  8. V. Bonacčić-Koutecký, J. Koutecký and J. Michl, Neutral and charged biradicals, zwitterions, funnels in S1, and proton translocation: Their role in photochemistry, photophysics, and vision, Angew. Chem., Int. Ed., 1987, 26, 170–189 CrossRef.
  9. L. V. Slipchenko and A. I. Krylov, Singlet-triplet gaps in diradicals by the spin-flip approach: A benchmark study, J. Chem. Phys., 2002, 117, 4694–4708 CrossRef CAS.
  10. A. I. Krylov, Quantum Chemistry of Open-Shell Species, in Reviews in Comp. Chem., ed. A. L. Parrill and K. B. Lipkowitz, J. Wiley & Sons, 2017, vol. 30 Search PubMed.
  11. N. Orms, D. Rehn, A. Dreuw and A. I. Krylov, Characterizing bonding patterns in diradicals and triradicals by density-based wave function analysis: A uniform approach, J. Chem. Theory Comput., 2018 DOI:10.1021/acs.jctc.7b01012.
  12. F. Neese, in High Resolution EPR: Applications to Metalloenzymes and Metals in Medicine, ed. L. Berliner and G. Hanson, Springer New York, New York, NY, 2009, p. 175 Search PubMed.
  13. J. P. Malrieu, R. Caballol, C. J. Calzado, C. de Graaf and N. Guihéry, Magnetic interactions in molecules and highly correlated materials: Physical content, analytical derivation, and rigorous extraction of magnetic Hamiltonians, Chem. Rev., 2013, 114, 429–492 CrossRef PubMed.
  14. W. Heisenberg, Zur Theorie des Ferromagnetismus, Z. Phys., 1928, 49, 619 CrossRef CAS.
  15. P. A. M. Dirac, On the theory of quantum mechanics, Proc. R. Soc. London, Ser. A, 1926, 112, 661 CrossRef CAS.
  16. J. H. vanVleck, The Theory of Electric and Magnetic Susceptibilities, Clarendon Press, Oxford, 1932 Search PubMed.
  17. A. Landé, Termstruktur und Zeemaneffekt der Multipletts, Z. Phys., 1923, 15, 189–205 CrossRef.
  18. N. J. Mayhall and M. Head-Gordon, Computational quantum chemistry for single Heisenberg spin couplings made simple: Just one spin flip required, J. Chem. Phys., 2014, 141, 134111 CrossRef PubMed.
  19. N. J. Mayhall and M. Head-Gordon, Computational quantum chemistry for multiple-site Heisenberg spin couplings made simple: Still only one spin-flip required, J. Phys. Chem. Lett., 2015, 6, 1982–1988 CrossRef CAS PubMed.
  20. L. Noodleman, Valence bond description of antiferromagnetic coupling in transition metal dimers, J. Chem. Phys., 1981, 74, 5737–5743 CrossRef CAS.
  21. S. Sinnecker, F. Neese, L. Noodleman and W. Lubitz, Calculating the electron paramagnetic resonance parameters of exchange coupled transition metal complexes using broken symmetry density functional theory: application to a Mn(III)/Mn(IV) model compound, J. Am. Chem. Soc., 2004, 126, 2613–2622 CrossRef CAS PubMed.
  22. I. P. R. Illas, R. Moreira, M. Costa and F. Filatov, Restricted ensemble-referenced Kohn–Sham versus broken symmetry approaches in density functional theory: Magnetic coupling in Cu binuclear complexes, J. Chem. Theory Comput., 2007, 3, 764–774 CrossRef PubMed.
  23. M. Filatov and S. Shaik, Chem. Phys. Lett., 1998, 288, 689 CrossRef CAS.
  24. M. Filatov and S. Shaik, J. Chem. Phys., 1999, 110, 116 CrossRef CAS.
  25. H. R. Zhekova, M. Seth and T. Ziegler, Calculation of the exchange coupling constants of copper binuclear systems based on spin-flip constricted variational density functional theory, J. Chem. Phys., 2011, 135, 184105 CrossRef PubMed.
  26. J.-M. Mouesca, J. L. Chen, L. Noodleman, D. Bashford and D. A. Case, Density functional/Poisson-Boltzmann calculations of redox potentials for iron-sulfur clusters, J. Am. Chem. Soc., 1994, 116, 11898 CrossRef CAS.
  27. M. Nishino, S. Yamanaka, Y. Yoshioka and K. Yamaguchi, Theoretical approaches to direct exchange couplings between divalent chromium ions in naked dimers, tetramers, and clusters, J. Phys. Chem. A, 1997, 101, 705–712 CrossRef CAS.
  28. E. Ruiz, J. Cano, S. Alvarez and P. Alemany, Broken symmetry approach to calculation of exchange coupling constants for homobinuclear and heterobinuclear transition metal complexes, J. Comput. Chem., 1999, 20, 1391–1400 CrossRef CAS.
  29. T. Soda, Y. Kitagawa, T. Onishi, Y. Takano, Y. Shigeta, H. Nagao, Y. Yoshioka and K. Yamaguchi, Ab initio computations of effective exchange integrals for HH, HHeH and Mn2O2 complex: comparison of broken-symmetry approaches, Chem. Phys. Lett., 2000, 319, 223–230 CrossRef CAS.
  30. A. Kazaryan, J. Heuver and M. Filatov, Excitation energies from spin-restricted ensemble-referenced Kohn-Sham method: A state-average approach, J. Phys. Chem. A, 2008, 112, 12980–12988 CrossRef CAS PubMed.
  31. A. I. Krylov, The spin-flip equation-of-motion coupled-cluster electronic structure method for a description of excited states, bond-breaking, diradicals, and triradicals, Acc. Chem. Res., 2006, 39, 83–91 CrossRef CAS PubMed.
  32. A. I. Krylov, Size-consistent wave functions for bond-breaking: The equation-of-motion spin-flip model, Chem. Phys. Lett., 2001, 338, 375–384 CrossRef CAS.
  33. S. V. Levchenko and A. I. Krylov, Equation-of-motion spin-flip coupled-cluster model with single and double substitutions: Theory and application to cyclobutadiene, J. Chem. Phys., 2004, 120, 175–185 CrossRef CAS PubMed.
  34. Y. Shao, M. Head-Gordon and A. I. Krylov, The spin-flip approach within time-dependent density functional theory: Theory and applications to diradicals, J. Chem. Phys., 2003, 118, 4807–4818 CrossRef CAS.
  35. Y. A. Bernard, Y. Shao and A. I. Krylov, General formulation of spin-flip time-dependent density functional theory using non-collinear kernels: Theory, implementation, and benchmarks, J. Chem. Phys., 2012, 136, 204103 CrossRef PubMed.
  36. D. Lefrancois, M. Wormit and A. Dreuw, Adapting algebraic diagrammatic construction schemes for the polarization propagator to problems with multi-reference electronic ground states exploiting the spin-flip ansatz, J. Chem. Phys., 2015, 143, 124107 CrossRef PubMed.
  37. L. V. Slipchenko and A. I. Krylov, Electronic structure of the trimethylenemethane diradical in its ground and electronically excited states: Bonding, equilibrium structures and vibrational frequencies, J. Chem. Phys., 2003, 118, 6874–6883 CrossRef CAS.
  38. L. V. Slipchenko and A. I. Krylov, Electronic structure of the 1,3,5,-tridehydrobenzene triradical in its ground and excited states, J. Chem. Phys., 2003, 118, 9614–9622 CrossRef CAS.
  39. L. V. Slipchenko, T. E. Munsch, P. G. Wenthold and A. I. Krylov, 5-dehydro-1,3-quinodimethane: A hydrocarbon with an open-shell doublet ground state, Angew. Chem., Int. Ed., 2004, 43, 742–745 CrossRef CAS PubMed.
  40. V. Vanovschi, A. I. Krylov and P. G. Wenthold, Structure, vibrational frequencies, ionization energies, and photoelectron spectrum of the para-benzyne radical anion, Theor. Chim. Acta, 2008, 120, 45–58 CrossRef CAS.
  41. E. Hossain, S. M. Deng, S. Gozem, A. I. Krylov, X.-B. Wang and P. G. Wenthold, Photoelectron spectroscopy study of quinonimides, J. Am. Chem. Soc., 2017, 139, 11138–11148 CrossRef CAS PubMed.
  42. D. Lefrancois, D. Rehn and A. Dreuw, Accurate adiabatic singlet-triplet gaps in atoms and molecules employing the third- order spin-flip algebraic diagrammatic construction scheme for the polarization propagator, J. Chem. Phys., 2016, 145, 084102 CrossRef PubMed.
  43. M. Head-Gordon, Characterizing unpaired electrons from the one-particle density matrix, Chem. Phys. Lett., 2003, 372, 508–511 CrossRef CAS.
  44. F. Plasser, M. Wormit and A. Dreuw, New tools for the systematic analysis and visualization of electronic excitations. I. Formalism, J. Chem. Phys., 2014, 141, 024106 CrossRef PubMed.
  45. A. V. Luzanov, in Practical Aspects of Computational Chemistry IV, J. Leszczynski and M. Shukla, ed. Springer-Verlag, New York, 2016 Search PubMed.
  46. R. Valero, F. Illas and D. G. Truhlar, Magnetic coupling in transition-metal binuclear complexes by spin-flip time-dependent density functional theory, J. Chem. Theory Comput., 2011, 7, 3523–3531 CrossRef CAS PubMed.
  47. I. Seidu, H. R. Zhekova, M. Seth and T. Ziegler, Calculation of exchange coupling constants in triply-bridged dinuclear Cu(II) compounds based on spin-flip constricted variational density functional theory, J. Phys. Chem. A, 2012, 116, 22682277 CrossRef PubMed.
  48. H. Zhao, C. Fang, J. Gao and C. Liu, Spin-state energies of heme-related models from spin-flip TDDFT calculations, Phys. Chem. Chem. Phys., 2016, 18, 29486–29494 RSC.
  49. S. Youngme, J. Phatchimkun, N. Wannarit, N. Chaichit and S. Meejoo, vanAlbada, G.A.; Reedijk, J. Refined crystal structure of tetra-p-acetato-bisaquodicopper(l1), Polyhedron, 2008, 27, 304–318 CrossRef CAS.
  50. O. Castillo, I. Muga, A. Luque, J. M. Gutierrez-Zorrilla, J. Sertucha, P. Vitoria and P. Roman, Synthesis, chemical characterization, X-ray crystal structure and magnetic properties of oxalato-bridged copper(II) binuclear complexes with 2,2′-bipyridine and diethylenetriamine as peripheral ligands, Polyhedron, 1999, 18, 1235 CrossRef CAS.
  51. P. deMeester, S. R. Fletcher and A. C. Skapski, Refined crystal structure of tetra-p-acetato-bisaquodicopper(l1), J. Chem. Soc., Dalton Trans., 1973, 23, 2575–2578 RSC.
  52. J. Sletten, The structures of two oxalato-bridged Cu dimers; [Cu2(Me4en)2(C2O4)(H2O)2](PF6)2·2H2O and [Cu2(Et5dien)2(C2O4)](PF6)2, Acta Chem. Scand., Ser. A, 1983, 37, 569 CrossRef.
  53. M. Julve, M. Verdaguer, A. Gleizes, M. Pohiloche-Levisalle and O. Kahn, Design of μ-oxalato copper(II) binuclear complexes exhibiting expected magnetic properties, Inorg. Chem., 1984, 23, 3808 CrossRef CAS.
  54. C. López, R. Costa, F. Illas, C. de Graaf, M. M. Turnbull, C. P. Landee, E. Espinosa, I. Mata and E. Molins, Magneto-structural correlations in binuclear copper(II) compounds bridged by a ferrocenecarboxylato(1) and an hydroxo- or methoxo-ligands, Dalton Trans., 2005, 2322 RSC.
  55. O. Castell, J. Miralles and R. Caballol, Structural dependence of the singlet-triplet gap in doubly bridged copper dimers: a variational CI calculation, Chem. Phys., 1994, 179, 377 CrossRef CAS.
  56. C. López, R. Costa, F. Illas, E. Molins and E. Espinosa, Ferromagnetic copper(II) complex containing ferrocenecarboxylato bridging ligands, Inorg. Chem., 2000, 39, 4560 CrossRef.
  57. T. Tokii, N. Hamamura, M. Nakashima and Y. Muto, Crystal structures and magnetic properties of novel μ-carboxylato-μ-hydroxo-bridged binuclear copper(II) complexes with 1,10-phenanthroline, Bull. Chem. Soc. Jpn., 1992, 65, 1214 CrossRef CAS.
  58. H. Zhekova, M. Seth and T. Ziegler, Introduction of a new theory for the calculation of magnetic coupling based on spinflip constricted variational density functional theory. Application to trinuclear copper complexes which model the native intermediate in multicopper oxidases, J. Chem. Theory Comput., 2011, 7, 18581866 CrossRef PubMed.
  59. A. I. Krylov, Triradicals, J. Phys. Chem. A, 2005, 109, 10638–10645 CrossRef CAS PubMed.
  60. In the SF approach, the high-spin triplet state is used as a “stepping stone”, to generate a balanced set of configurations describing the target-state manifold. The energy of the reference triplet is different from the Ms = 0 triplet due to their different treatment by the SF ansatz. Only in the limit of the full account of electron correlation, e.g., EOM-SF-CCSD for two-electron system, the reference and the Ms = 0 triplets become degenerate. In SF calculations, the energy gaps between singlet and triplets are computed using the energies of the target Ms = 0 states. The energy gap between the reference and the SF triplet state can be used to quantify the imperfections in the treatment of electron correlation in the reference state. Table S2 in ESI shows energy splittings for CUAQACO2 (complex 2).
  61. J. S. Sears, C. D. Sherrill and A. I. Krylov, A spin-complete version of the spin-flip approach to bond breaking: What is the impact of obtaining spin eigenfunctions?, J. Chem. Phys., 2003, 118, 9084–9094 CrossRef CAS.
  62. F. Wang and T. Ziegler, The performance of time-dependent density functional theory based on a noncollinear exchange-correlation potential in the calculations of excitation energies, J. Chem. Phys., 2005, 122, 074109 CrossRef PubMed.
  63. Z. Rinkevicius and H. Årgen, Spin-flip time dependent density functional theory for single–triplet splittings in σ,σ-biradicals, Chem. Phys. Lett., 2010, 491, 132–135 CrossRef CAS.
  64. Z. Rinkevicius, O. Vahtras and H. Årgen, Spin-flip time dependent density functional theory applied to excited states with single, double, or mixed electron excitaion character, J. Chem. Phys., 2010, 133, 114104 CrossRef PubMed.
  65. S. A. Kozyrev and B. M. Al'tshuler, Electron Paramagnetic Resonance, Academic Press, New York, 1964 Search PubMed.
  66. A. Caneschi, I. Mirabeau, M. Hennion, H. Casalta, H. Andres, H. U. Güdel and A. V. Irodova, Low-energy magnetic excitations of the Mn12-acetate spin cluster observed by neutron scattering, Phys. Rev. Lett., 1999, 83, 628 CrossRef.
  67. G. Amoretti, R. Caciuffo, J. Combet, A. Murani and A. Caneschi, Inelastic neutron scattering below 85 eV and zero-field splitting parameters in the Fe8 magnetic cluster, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 62, 3022 CrossRef CAS.
  68. J. A. Volkert, C. M. Larrabee, E. T. Alessi, J. O. Asiedu, K. R. Cook, L. J. Hoerning, G. S. Klinger, S. G. Okin and T. L. Santee, Magnetic circular dichroism spectroscopy as a probe of geometric and electronic structure of cobalt(II)-substituted proteins: ground-state zero-field splitting as a coordination number indicator, J. Am. Chem. Soc., 1997, 119, 4182 CrossRef.
  69. M. Atanasov, D. Aravena, E. Suturina, E. Bill, D. Maganas and F. Neese, First prinicples approach to the electronic structure, magnetic anisotropy and spin relaxation in mononuclear 3d-transition metal single molecular magnets, Coord. Chem. Rev., 2015, 289–290, 177–214 CrossRef CAS.
  70. F. Neese and E. I. Solomon, Calculation of zero-field splittings, g-values, and the relativistic nephelauxetic effect in transition metal complexes. Application to high-spin ferric complexes, Inorg. Chem., 1998, 37, 6568 CrossRef CAS PubMed.
  71. S. Sinnecker and F. Neese, Spin–spin contributions to the zero-field splitting tensor in organic triplets, carbenes and biradicals. A density functional and ab initio study, J. Phys. Chem. A, 2006, 110, 12267–12275 CrossRef CAS PubMed.
  72. D. Ganyushin and F. Nesse, First-principles calculations of zero-field splitting parameters, J. Chem. Phys., 2006, 125, 024103 CrossRef PubMed.
  73. H. Zhekova, M. Seth and T. Ziegler, Calculation of the exchange coupling constants of copper binuclear systems based on spin-flip constricted variational density functional theory, J. Chem. Phys., 2012, 135, 184105 CrossRef PubMed.
  74. M. Muzarová and M. Kaupp, A critical validation of density functional and coupled-cluster approaches for the calculation of EPR hyperfine coupling constants in transition metal complexes, J. Phys. Chem. A, 1999, 103, 9966–9983 CrossRef.
  75. T. N. Yanai, Y. Lan and T. Kurashige, Scalar relativistic calculations of hyperfine coupling constants using ab initio density matrix renormalization group method in combination with third-order DouglasKrollHess transformation: Case studies on 4d transition metals, J. Chem. Theory Comput., 2014, 11, 73–81 Search PubMed.
  76. B. Bleaney and K. D. Bowers, Anomalous paramagnetism of copper acetate, Proc. R. Soc. London, Ser. A, 1952, 214, 451–465 CrossRef CAS.
  77. The effect of ionic configurations is included in the HDvV exchange-couplings in an effective fashion, via perturbation theory. The perturbative treatment assumes that their contributions are small.
  78. J.-D. Chai and M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom-atom dispersion interactions, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620 RSC.
  79. A. D. Becke, A multicenter numerical integration scheme for polyatomic molecules, J. Chem. Phys., 1988, 88, 2547 CrossRef CAS.
  80. C. Lee, W. Yang and R. G. Parr, Development of the Colle-Salvetti correlation-energy formula into a functional of the electron-density, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789 CrossRef CAS.
  81. A. D. Becke, Density-functional thermochemistry. iii. The role of exact exchange, J. Chem. Phys., 1993, 98, 5648 CrossRef CAS.
  82. A. D. Becke, Density-functional exchange-energy approximation with correct asymptotic-behavior, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100 CrossRef CAS.
  83. J. P. Perdew, Phys. Rev. B: Condens. Matter Mater. Phys., 1986, 33, 8822 CrossRef.
  84. J. P. Perdew In Electronic Structure of Solids, ed. P. Ziesche and H. Eschrig, Akadernie Verlag, Berlin, 1991 Search PubMed.
  85. J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh and C. Fiolhais, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 46, 6671 CrossRef CAS.
  86. J. P. Perdew and Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244–13249 CrossRef.
  87. A. D. Becke, Density functional thermochemistry. II. The effect of the Perdew Wang generalized gradient correlation correction, J. Chem. Phys., 1992, 97, 9173–9177 CrossRef CAS.
  88. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868 CrossRef CAS PubMed.
  89. C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys., 1999, 110, 6158–6170 CrossRef CAS.
  90. Y. Zhao and D. G. Truhlar, 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., 2008, 120, 215–241 CrossRef CAS.
  91. H. S. Yu, W. Zhang, P. Verma, X. He and D. G. Truhlar, Nonseparable exchange–correlation functional for molecules, including homogeneous catalysis involving transition metals, Phys. Chem. Chem. Phys., 2015, 17, 12146–12160 RSC.
  92. H. S. Yu, X. He and D. G. Truhlar, MN15-L: A new local exchange-correlation functional for Kohn–Sham density functional theory with broad accuracy for atoms, molecules, and solids, J. Chem. Theory Comput., 2016, 12, 1280–1293 CrossRef CAS PubMed.
  93. H. S. Yu, S. Haoyu, X. He, S. L. Li and D. G. Truhlar, MN15: A Kohn–Sham global-hybrid exchange-correlation density functional with broad accuracy for multi-reference and single-reference systems and noncovalent interactions, Chem. Sci., 2016, 7, 5032–5051 RSC.
  94. Y. Shao, Z. Gan, E. Epifanovsky, A. T. B. Gilbert, M. Wormit, J. Kussmann, A. W. Lange, A. Behn, J. Deng and X. Feng, et al., Advances in molecular quantum chemistry contained in the Q-Chem 4 program package, Mol. Phys., 2015, 113, 184–215 CrossRef CAS.
  95. A. I. Krylov and P. M. W. Gill, Q-Chem: An engine for innovation, WIREs Comput. Mol. Sci., 2013, 3, 317–326 CrossRef CAS.
  96. F. Wang and T. Ziegler, Time-dependent density functional theory based on a noncollinear formulation of the exchange-correlation poential, J. Chem. Phys., 2004, 121, 12191 CrossRef CAS PubMed.
  97. O. Vahtras and Z. Rinkevicius, General excitations in time-dependent density functional theory, J. Chem. Phys., 2007, 126, 114101 CrossRef PubMed.
  98. P. J. Hay and W. R. Wadt, Ab initio effective core potentials for molecular calculations. potentials for K to Au including the outermost core orbitals, J. Chem. Phys., 1985, 82, 299 CrossRef CAS.
  99. C. Hartwigsen, S. Goedecker and J. Hutter, Relativistic separable dual-space Gaussian pseudopotentials from H to Rn, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 58, 3641–3662 CrossRef CAS.
  100. M. M. Hurley, L. F. Pacios, P. A. Christiansen, R. B. Ross and W. C. Ermler, Ab initio relativistic effective potentials with spinorbit operators. II. K through Kr, J. Chem. Phys., 1986, 84, 6840 CrossRef CAS.
  101. L. F. Pacios and P. A. Christiansen, Ab initio relativistic effective potentials with spinorbit operators. I. Li through Ar, J. Chem. Phys., 1985, 82, 2664 CrossRef.
  102. D. Figgen, G. Rauhut, M. Dolg and H. Stoll, Energy-consistent pseudopotentials for group 11 and 12 atoms: adjustment to multi-configuration Dirac-Hartree-Fock data, Chem. Phys., 2005, 311, 227–244 CrossRef CAS.
  103. K. A. Peterson and C. Puzzarini, Systematically convergent basis sets for transition metals. II. Pseudopotential-based correlation consistent basis sets for the group 11 (Cu, Ag, Au) and 12 (Zn, Cd, Hg) elements, Theor. Chem. Acc., 2005, 114, 283–296 CrossRef CAS.
  104. T.-C. Jagau, D. B. Dao, N. S. Holtgrewe, A. I. Krylov and R. Mabbs, Same but different: Dipole-stabilized shape resonances in CuF and AgF, J. Phys. Chem. Lett., 2015, 6, 2786–2793 CrossRef CAS PubMed.
  105. N. Orms and A. I. Krylov, Modeling photoelectron spectra of CuO, Cu2O, and CuO2 anions with equation-of-motion coupled-cluster methods: An adventure in Fock space, J. Phys. Chem. A, 2018 DOI:10.1021/acs.jpca.7b10620.
  106. IQmol molecular viewer, A. T. B. Gilbert,
  107. F. H. Allen, The CSD system: The Cambridge structural database: A quarter of a million crystal structures and rising, Acta Crystallogr., Sect. B: Struct. Sci., 2002, B58, 380–388 CrossRef CAS.
  108. I. J. Bruno, J. C. Cole, P. R. Edgington, M. Kessler, C. F. Macrae, P. McCabe, J. Pearson and R. Taylor, ConQuest: New software for searching the Cambridge structural database and visualizing crystal structures, Acta Crystallogr., Sect. B: Struct. Sci., 2002, B58, 389 CrossRef CAS.


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