What Are the Physical Contents of Hubbard and Heisenberg

Article pubs.acs.org/JCTC

What Are the Physical Contents of Hubbard and Heisenberg Hamiltonian Interactions Extracted from Broken Symmetry DFT Calculations in Magnetic Compounds? Grégoire David,† Nathalie Guihéry,*,‡ and Nicolas Ferré† †

Aix Marseille Univ, CNRS, ICR, 13397 Marseille, France LCPQ, IRSAMC, Université de Toulouse 3, Paul Sabatier, 31400 Toulouse, France

‡

S Supporting Information *

ABSTRACT: Analytical expressions of the interactions present in the Heisenberg−Dirac van Vleck and Hubbard Hamiltonians have been derived as functions of both the energy of several broken symmetry DFT solutions and their expectation value of the S2 spin operator. Then, following a strategy of decomposition of the magnetic exchange coupling into its main contributions (direct exchange, kinetic exchange, and spin polarization) and using a recently proposed method of spin decontamination, values of these interactions have been extracted. As already observed, they weakly depend on the correlation functional but strongly depend on the exchange one. In order to distinguish between the effect of the delocalization of the magnetic orbitals and that of the amount of Hartree−Fock exchange (HFX) when hybrid exchange-correlation functionals are used, we have disentangled these two contributions by either freezing the magnetic orbitals and varying the amount of HFX or varying the magnetic orbitals while keeping the same amount of HFX. As expected, increasing the amount of HFX induces a slight relocalization of the magnetic orbitals on the magnetic center which results in a weak increase of the repulsion energy U parameter and a weak decrease of both the direct exchange Kab and hopping |t| parameters. Conversely, the amount of HFX has a huge effect on all the parameters, even when some of the parameters should be exchange-independent, like U. Indeed, it is analytically demonstrated that the physical content of the U parameter extracted from several broken-symmetry solutions depends on the amount of HFX and that this pathological behavior has the same origin as the self-interaction error. This result is interesting not only to theoretical chemists working in the field of magnetic systems but also to DFT methodologists interested in using this theory for studying either excited states or strongly correlated systems. Finally, the performance of the range-separated ωB97XD functional for both ferromagnetic and antiferromagnetic transition-metal compounds and organic systems must be noted.

■

INTRODUCTION Remarkable properties, such as magnetism, superconductivity, magnetoresistive effects, and magnetic anisotropy, find their origin in the quantum behavior of matter. Some of them are still not fully understood and arouse great interest in fundamental research. The systems that house these properties are usually difficult to treat theoretically as they are highly correlated (i.e. with a strongly multireference character) and often of large size. Consequently, both experimentalists and theoreticians use model Hamiltonians, such as Heisenberg− Dirac-van Vleck (HDvV), Hubbard, t-J, etc., in order to characterize them.1 The advantage of these models over the allelectron exact Hamiltonian is 2-fold: (i) the number of explicitly considered electronic configurations is dramatically reduced as only a few orbitals and electrons are explicitly considered (the magnetic electrons in their magnetic orbitals for instance). The subsequent use of these models on larger systems is therefore possible enabling one to introduce © 2017 American Chemical Society

collective effects at work in large systems or highly correlated materials. (ii) The model interactions are in limited number and provide a good understanding of the physics. Indeed, they may capture complex physics within rather simple mechanisms, such as charge hopping or spin exchange, which is helpful to streamline the properties. One must stress that, as they are not observable, their experimental determination is indirect and may be subject to controversy.2 The methods of choice, in order to extract them from theoretical calculations,3 are based on the wave function (WF) theory, and one may in particular quote the Difference Dedicated Configuration Interaction4 which is currently the best available method for this purpose. Indeed, these methods provide the energies and wave functions of the multireference lowest states. A direct identification between the energies of the model states and those of the Received: September 19, 2017 Published: October 17, 2017 6253

DOI: 10.1021/acs.jctc.7b00976 J. Chem. Theory Comput. 2017, 13, 6253−6265

Article

Journal of Chemical Theory and Computation

rationalize the performances of the DFT functionals and therefore help in guiding the choice of the appropriated ones is to provide tools to analyze the physical content of model effective interactions. In contrast with multireference wave functions, the analysis of BS-DFT solutions is not straightforward as they may simultaneously incorporate all or part of several physical effects, thus no longer allowing their distinction or differentiation. The present paper aims at providing an appropriate procedure of analysis. In order to control at best the physical effects introduced in the calculated interactions, we have developed a strategy that combines (i) analytical derivations of the interactions as functions of energy differences between various BS-DFT solutions and their expectation values of the S2 electron spin operator and (ii) the use of selective freezing of orbital sets that has shown to provide good decompositions of magnetic exchange couplings into their main contributions45,46 and the use of a recently proposed method of spin decontamination47 that enables one to accurately account for spin polarization effects. Finally, the roles of the spin delocalization and of %HFX have been analyzed separately in order to appreciate how each contribution affects the interactions. It is worth noting that this procedure of analysis could be used to determine the physical content in BS-DFT of more subtle interactions such as three-body or four-body ones and the interactions of other model Hamiltonians (double exchange, anisotropic spin models, etc.). In each case new analytical derivations would be required. The paper is organized as follows. The next section derives the analytical expressions of both the Heisenberg and Hubbard interactions and recalls the method of freezing of orbitals from which the decomposition of the interactions into their main contributions becomes possible. The results are discussed in section III. Summary and perspectives are given in section IV.

computed states can therefore be performed, from which all interactions can usually be extracted. Another advantage of WFbased methods is the possibility to decompose the interactions in their various contributions and to analyze the physical effects from the weight of the determinants introduced in the configuration interaction expansion. Although WF methods are quite expensive computationally, extractions of HdvV,5,6 anisotropic spin,7,8 t-J,9 double exchange10,11 models have been successfully performed in both molecular systems and highly correlated materials. On the other hand, the broken-symmetry (BS) Kohn−Sham density functional theory (DFT) approach combined with a strategy first proposed by Noodleman12 appears as a powerful tool for the evaluation of model Hamiltonian parameters. While its applicability to large systems is probably its main advantage over WF methods, it is also particularly suited for the calculation of subtle interactions such as the biquadratic exchange in systems of spin S = 113 or the double exchange interactions14 and for geometry optimizations.15,16 Its performances are well documented in several review papers.17−21 Nevertheless, the values of the parameters depend on the choice of the exchange-correlation functional. Moreover, the spin-contamination problem of the broken-symmetry solutions severely hampers its routine use.22−24 Concerning magnetic compounds, previous works have shown that the computed magnetic couplings strongly depend on the amount of Hartree−Fock exchange (hereafter simply denoted %HFX) when hybrid functionals are used.24−35 The reason invoked in several papers is the exceeding spin delocalization obtained with most of the functionals which results in overestimations of both antiferromagnetic and ferromagnetic couplings.36−39 The article of Calzado et al.5 provides an explanation for the overestimation of antiferromagnetic couplings, showing that the exceeding delocalization resulting from the use of a too small % HFX leads to too small values of the magnetic electron repulsion (the U parameter of the Hubbard Hamiltonian) and therefore to an overestimation of the kinetic exchange contribution. On the other hand, as exceeding delocalization may also play in favor of a larger ferromagnetic direct exchange contribution between the magnetic electrons, it could also rationalize the overprediction of ferromagnetic couplings. It was also noted that the %HFX needed to reproduce the experimental values is close to 33%.5,23,40−43 Phillips et al.44 have proposed that spin delocalization is probably not the single source of errors as the nonself-consistent calculation of magnetic couplings with semilocal functionals does not give a balanced improvement of both ferro- and antiferromagnetic couplings. The physical content of magnetic couplings and, more generally, of effective interactions is rather complex as it results from subtle effects of electron correlation. For instance, while the direct exchange always provides a ferromagnetic contribution, the spin polarization which is also captured in the effective K parameter of the generalized Hubbard Hamiltonian might be either ferromagnetic (parallel alignment of the magnetic spins) or antiferromagnetic (antiparallel alignment of the magnetic spins) and of much larger magnitude than the direct exchange. Likewise U, that reflects the impact of the electron correlation, proves to be an extremely sensitive parameter to the choice of the functional. This may seem paradoxical as it is a purely Coulombic term in Hartree−Fock theory that would thus be expected to be exactly calculated in DFT. An important issue if one wants to understand the various physical effects that

■

DERIVATION AND EXTRACTION OF THE HDVV AND HUBBARD MODELS FROM BROKEN-SYMMETRY DFT SOLUTIONS General Presentation of the Models. For a system containing two magnetic centers A and B, the Heisenberg− Dirac van Vleck Hamiltonian writes HDvV Ĥ = −JSA⃗̂ . SB⃗̂ (1) ̂ ̂ where S⃗A (respectively S⃗B) is the spin momentum operator of center A (respectively B), and J is the amplitude of the magnetic coupling between the two centers. In this expression J is negative for an antiferromagnetic coupling. The three main contributions to the magnetic coupling are (i) the ferromagnetic direct exchange Kab between the magnetic orbitals a and b localized on centers A and B, (ii) the antiferromagnetic kinetic exchange (or superexchange) provided by the intersite delocalization of the magnetic orbitals, and (iii) the spin and charge polarizations effects brought by the nonmagnetic orbitals. This Hamiltonian can be derived from the generalized Hubbard Hamiltonian written as48 Hubbard Ĥ = K (aA+, ↑aB+, ↓aA , ↓aB , ↑ + H.c.) + t(aA+, MsaB , Ms + H.c.)

+U

∑ i=A ,B

ni , ↑ni , ↓ (2)

in the case of two centers bearing a spin S = 1/2. H.c. stands for Hermitian conjugate. The effective exchange integral K 6254


Article

Journal of Chemical Theory and Computation accounts for both the direct exchange and the spin polarization effects, t is the hopping integral between the two centers, â†A,Ms and âB,Ms are respectively the creation and annihilation operators of an electron at sites A and B with Ms spin component, and U is the energy difference between the repulsion of the two electrons occupying the same local magnetic orbital and that exerted when singly occupying the two local magnetic orbitals. The charge polarization effects result in a screening of the value of U. n̂i,Ms is the occupation number operator. With the limit t/U being sufficiently small, i.e. for magnetic systems, the HDvV Hamiltonian becomes relevant for the description of the lowest states of the Hubbard Hamiltonian that are essentially spread over neutral configurations, while the ionic configurations effects are accounted for by the antiferromagnetic contribution to the magnetic coupling. Extraction of the Various Interactions and Contributions to the HDvV and Hubbard Models. The method of extraction presented in refs 45 and 46 proceeds through consecutive freezings and relaxations of different classes of orbitals. In the first place, a restricted open shell (RO) DFT solution is computed for the Ms = 1 spin component (with ⟨S2⟩T,RO = 2) ΦT,RO =

∏ i i ̅gu

=

i

where UFC means unrestricted formalism using frozen core orbitals, and a′ and b′ are partially delocalized orbitals that contain delocalization tails of each magnetic orbital on the other magnetic center:49 a′ = a cos φ + b sin φ b′ = a sin φ + b cos φ

The energy of this solution as a function of the Hubbard parameters is given by E BS,UFC =

g+u 2

and

b=

a′ = a cos φ + b sin φ b′ = b cos φ + a sin φ

ΔJKE = =

(4)

|t | =

i

2 (2U cos 2 φ sin 2 φ + 4t cos φ sin φ) 1 + 4 cos2 φ sin 2 φ (12)

(E T,RO − E BS,UFC) + Kab

U=2

1 − ⟨S2⟩BS,UFC

(13)

(E T,RO − E BS,UFC) + Kab 1 − ⟨S2⟩BS,UFC

− 2Kab

(14)

Finally, the spin polarization is introduced by relaxing the core orbitals in the field of the magnetic a and b orbitals for the Ms = 1 solution and in the field of the partially delocalized magnetic orbitals a′ and b′ for the Ms = 0 one

(6)

=

2 (E BS,UFC − E T,RO) − 2Kab ⟨S2⟩T,RO − ⟨S2⟩BS,UFC

Using eqs 10 and 11 in eqs 9 and 12 one gets

In order to calculate the kinetic exchange energy responsible for the antiferromagnetic contribution and coming from the ionic valence bond forms |aa|̅ and |bb̅|, the magnetic orbitals are relaxed for the Ms = 0 solution while keeping the core orbitals frozen

∏ i i ̅ a′ b ̅ ′

(11)

The kinetic exchange energy that is responsible for the antiferromagnetic contribution to the magnetic coupling is

with ⟨S2⟩BS,RO = 1. From the energy difference between these two previous solutions, the direct exchange integral Kab can be extracted:

ΦBS,UFC =

(10)

The expectation value of S can also be expressed as a function of the rotation angle φ:

(5)

Kab = E BS,RO − E T,RO

∂E ∂φ

2

∏ i i ̅ab ̅ i

i

= 0 leads to the equation −t cos φ sin φ = U + 2Kab

and its minimization

Then a spin flip of one of the magnetic electrons leads to the Ms = 0 spin solution computed in the set of RO-DFT frozen orbitals ΦBS,RO =

∏ i i ̅ a′ b ̅ ′ (9)

(3)

g−u 2

Ĥ

= E BS,RO + (2U + 4Kab)cos 2 φ sin 2 φ + 4t cos φ sin φ

where i denotes all the doubly occupied nonmagnetic (hereafter called core) orbitals, g and u are the delocalized (canonical) magnetic orbitals, and a and b are the localized magnetic orbitals obtained through the rotations: a=

∏ i i ̅ a′ b ̅ ′ i

∏ i i ̅ab i

(8)

ΦT,UFM =

∏ jα jβ ab j

ΦBS,UFM =

∏ kα kβa′b ̅ ′ k

∏ i i ̅(ab ̅ cos2 φ + ba ̅ sin 2 φ

(15)

(16)

where UFM means unrestricted formalism using frozen magnetic orbitals, while jα, jβ, kα, and kβ denote the relaxed core orbitals of the Ms = 1 and Ms = 0 solutions. Using a recently proposed method of spin decontamination,47 the spin polarization contribution is

i

+ (aa ̅ + bb ̅ )sin φ cos φ) (7)

6255


Article

Journal of Chemical Theory and Computation

Figure 1. Schematic representation of the studied molecular complexes. An experimental (when available) or theoretical magnetic exchange coupling is reported for 1,51 2,52−56 3,57 4,58 5,6 6,59 and 7.1

ΔJSP =

2(E BS,UFM − E T,UFM) 2 − (⟨S

2

⟩ + ⟨S

BS,UFC

2 BS,UFM

⟩)/2 + ⟨S2 BS,UFC⟩(⟨S2 BS,UFM⟩ − ⟨S2 BS,UFC⟩)/2

ΔJSP 2

(18)

According to the HDvV Hamiltonian (eq 1), the magnetic exchange can be calculated from all these contributions: JCALC = 2Kab + ΔJKE + ΔJSP

(19)

Alternatively, the spin decontaminated magnetic coupling can also be determined directly from the Ms = 0 and Ms = 1 BS solutions47 JCOMP = [2(E BS,U − E T,U)]/[2 − (⟨S2⟩BS,UFC + ⟨S2⟩BS,U ) /2 + ⟨S2⟩BS,UFC (⟨S2⟩BS,U − ⟨S2⟩BS,UFC )/2] (20)

where U means unrestricted formalism for all the orbitals (none are frozen). Finally, in order to estimate the consistency of the extraction, the magnetic coupling can be calculated from the Hubbard parameters: U 2 + 16t 2 (21) 2 One should note that for the UFC solution, the Yamaguchi procedure for spin decontamination50 has been safely used. Indeed, in the absence of spin polarization, there is no need for the most sophisticated method.47 Studied Systems and Computational Information. As the magnetic couplings and, more generally, the interactions of highly correlated molecular models strongly depend on the geometrical structure, either available experimental X-ray diffraction geometries or optimized geometries when the corresponding X-ray structure was not available have been used for the following studied complexes. Ph(HNO)251 and C6H415,16,52−56 labeled 1 and 2, respectively, are strongly antiferromagnetic organic compounds. [{Cu(H2 O)}2 (μAcO)4],57 labeled 3, is the well-known copper acetate molecule in which two copper ions are bridged by four acetato groups in a paddle-wheel core. The [{Cu(dpt)}]2{μ-O2C-(η5-C5H4)Fe(η5-C5H5)}2](ClO4)2 system,58 labeled 4, consists of two JIC = 2K +

(17)

Cu(dpt) moieties (dpt = dimethylpropilenetriamine) bridged by two ferrocenecarboxylato ligands. The square-pyramidal coordination, with parallel basal planes, of the copper ion results in parallel magnetic orbitals unfavorable to the superexchange pathways, thus resulting in a weak ferromagnetic coupling. In the strongly antiferromagnetic Cu2(N3)2(NH3)6,6 labeled 5, the two copper ions are bridged by two end-to-end azido groups. [(Et5dien)2Cu2(μ-C2O4)]2+ion (Et5dien = pentaethyldiethylenetriamine),59 labeled 6, is made of two copper ions bridged through an oxalato bis-chelating ligand. Finally, the Cu2Cl62− ion,1,6 labeled 7, planar units exhibit a weak antiferromagnetic coupling. A schematic representation of these complexes can be found in Figure 1. Calculations have been performed with the Gaussian package60 and using selective freezing of orbitals thanks to the Local Self-Consistent-Field (LSCF) method.61 We have used several GGA functionals namely PBE-PBE,62 PBEVWN5,63 TPSS-TPSS,64 and TPSS-VWN5. In this FX-FC notation, the first acronym FX denotes the exchange functional, while the second one FC denotes the correlation functional. To appreciate the effect of %HFX, the hybrid GGA B3LYP65 functional is used varying the α parameter in the calculation of the exchange−correlation energy according to the equation

The exchange parameter of the Hubbard Hamiltonian, accounting for the spin polarization, is therefore K = Kab +

− 2Kab − ΔJKE

U−

E XC = αE XHF + (1 − α)E XLocal + βE XNonlocal + ECDFT

(22)

Local where EHF X is calculated using Hartree−Fock exchange, EX 66 Nonlocal and EX (β = 0.72) are computed with DFT exchange, and EDFT is the correlation energy computed at the DFT level. C

We have also considered the BH&HLYP functional which features 50%HFX, as it is implemented in the Gaussian09 package.60 Finally we have tested the long-range-corrected hybrid functional ωB97XD.67 In short, the long-range exchange is treated with Hartree−Fock, while the short-range one is treated with a hybrid density functional approximation. The basis set is 6-31G(d,p) for all atoms. All figures of molecular structures and orbitals have been obtained by using the JMol software.68

■

RESULTS AND DISCUSSION Analysis of the Self-Consistent Results. For compounds 1−7, values of JCALC (eq 19), JCOMP (eq 20), and JIC (eq 21) are reported in Table 1. While the extracted values for the Hubbard 6256


Article

Journal of Chemical Theory and Computation Table 1. Magnetic Couplings in cm−1a 1

2

3

4

5

6

7

JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF JCALC JIC JCOMP JREF

PBE PBE

PBE VWN5

TPSS TPSS

TPSS VWN5

B3LYP_20%

B3LYP_30%

B3LYP_50%

BH&H LYP

ωB97XD

−4192 −4791 −4342

−4115 −4733 −4252

−3946 −4287 −4156

−3879 −4228 −4079

−2737 −2713 −3338

−2057 −2059 −3279

−2088 −2090 −3217

−2529 −2502 −3309

−3272 −3346 −3332

−3304 −3390 −3357

−2917 −2951 −2975

−2787 −2823 −2847

−2200 −2206 −2264

−1879 −1880 −1951

−1956 −1958 −2028

−2383 −2389 −2464

−904 −909 −940

−891 −897 −949

−773 −775 −825

−777 −780 −828

−234 −234 −251

−106 −106 −119

−114 −114 −128

−256 −255 −271

18 17 12

16 16 10

18 18 12

16 15 11

−3186 −3201 −3608 ≃−200051 −2530 −2544 −2592 −242052−56 −354 −353 −379 −28657 3 3 1 258 −1794 −1798 −2023

What Are the Physical Contents of Hubbard and Heisenberg

Recommend Documents