Erik Donovan Hedegård

Julien Toulouse

Hans Jørgen

Aagaard Jensen

Multiconfigurational Short-Range Density-Functional Theory for Open-Shell Systems

Keywords: range-separated DFT, MCSCF, open shells

Many chemical systems cannot be described by quantum chemistry methods based on a single-reference wave function. Accurate predictions of energetic and spectroscopic properties require a delicate balance between describing the most important configurations (static correlation) and obtaining dynamical correlation efficiently. The former is most naturally done through a multiconfigurational wave function, whereas the latter can be done by e.g. perturbation theory. We have employed a different strategy, namely a hybrid between multiconfigurational wave functions and densityfunctional theory (DFT) based on range separation. The method is denoted multiconfigurational short-range DFT (MC-srDFT) and is more efficient than perturbative approaches as it capitalizes on the efficient treatment of the (short-range) dynamical correlation by DFT approximations. In turn, the method also improves DFT with standard approximations through the ability of multiconfigurational wave functions to recover large parts of the static correlation. Until now, our implementation was restricted to closed-shell systems and to lift this restriction, we present here the generalization of MC-srDFT to open-shell cases. The additional terms required to treat open-shell systems are derived and implemented in the DALTON program. This new method for open-shell systems is illustrated on dioxygen and [Fe(H 2 O) 6 ] 3+ .

I. INTRODUCTION

Quantum chemistry (QC) methods have become paramount to gain insight in experimental investigations and to support findings from experimental studies. Moreover, QC methods can be employed to predict novel molecular systems with a given desired property. Naturally, the above uses of theoretical methods require that these methods are both efficient and accurate, also in cases where reference data might not exist or be inconclusive. Density-functional theory (DFT) has to a large extent fulfilled these requirements [START_REF] Burke | Perspective on Density Functional Theory[END_REF] , and can for many systems be employed routinely. Still, DFT with standard approximations is known to fail for systems that cannot be described by a single-reference wave function. [START_REF] Cohen | Insights into Current Limitations of Density Functional Theory[END_REF][START_REF] Cohen | Challenges for density functional theory[END_REF] Such systems are often characterized by dense orbital manifolds and several (near)-degenerate states. In wave-function theory (WFT), such systems are best described with a multiconfigurational wave function. [START_REF] Szalay | Multiconfiguration Self-Consistent Field and Multireference Configuration Interaction methods and applications[END_REF] These exist in a number of different forms of which many have in common the definition of a complete active space CAS(m, n) of m electrons in n orbitals. In this active space, all configurations (fulfilling additional spin and symmetry constraints) are included. This is usually combined with optimization of orbital parameters in what is denoted a complete-activespace self-consistent field (CASSCF) procedure. Unfortunately, the computational demands increase dramatically with the size of the active space, and many systems require active spaces beyond the current limitations to give physically meaningful results. In recent years, several groups have focused on lifting the limitations for the size of active space with methods such as the density-matrix renormalization group (DMRG) [START_REF] Marti | The Density Matrix Renormalization Group Algorithm in Quantum Chemistry[END_REF][START_REF] Chan | The Density Matrix Renormalization Group in Quantum Chemistry[END_REF][START_REF] Kurashige | Multireference electron correlation methods with density matrix renormalisation group reference functions[END_REF][START_REF] Knecht | New Approaches for ab initio Calculations of Molecules with Strong Electron Correlation[END_REF] , quantum Monte Carlo (QMC) [START_REF] Filippi | Simple formalism for efficient derivatives and multideterminant expansions in Simple formalism for efficient derivatives and multi-determinant expansions in quantum Monte Carlo[END_REF][START_REF] Li Manni | Combining the Complete Active Space Self-Consistent Field Method and the Full Configuration Interaction Quantum Monte Carlo within a Super-CI Framework , with Application to Challenging Metal-Porphyrins[END_REF][START_REF] Thomas | Stochastic Multiconfigurational Self-Consistent Field Theory[END_REF] , restricted-active-space (RAS) [START_REF] Olsen | Determinant based Configuration-Interaction algorithms for complete and restricted Configuration-Interaction spaces[END_REF] or generalized-active-space (GAS) methods. [START_REF] Fleig | The generalized active space concept for the relativistic treatment of electron correlation. I. Kramers-restricted two-component configuration interaction[END_REF][START_REF] Thyssen | A direct relativistic four-component multiconfiguration self-consistent-field method for molecules[END_REF][START_REF] Ma | The Generalized Active space concept in Multiconfigurational Self-Consistent Field methods[END_REF] Yet, even with extended active orbital spaces, essential parts of the remaining dynamical electron correlation cannot be obtained, except for the smallest systems where the number of occupied valence orbitals is small. Typically, dynamical correlation for such multireference systems is obtained after initially obtaining a correct representation of a zeroth-order Hamiltonian (including static correlation). The exact nature of the subsequent steps responsible for recovering dynamical correlation depends on the chosen method, but well-known examples are multireference perturbation theories [START_REF] Andersson | Second-order perturbation theory with a CASSCF reference function[END_REF][START_REF] Andersson | Second-order perturbation theory with a Complete Active Space Self-Consistent Field reference function[END_REF][START_REF] Malmqvist | The restricted active space followed by second-order perturbation theory method: Theory and application to the study of CuO 2 and Cu 2 O 2 systems[END_REF][START_REF] Ma | Second-order perturbation theory for generalized active space self-consistent-field wave functions[END_REF][START_REF] Angeli | Introduction of n-electron valence states for multireference perturbation theory[END_REF] such as complete-active-space secondorder perturbation theory (CASPT2) [START_REF] Andersson | Second-order perturbation theory with a CASSCF reference function[END_REF][START_REF] Andersson | Second-order perturbation theory with a Complete Active Space Self-Consistent Field reference function[END_REF] and n-electron valence state perturbation theory (NEVPT2) [START_REF] Angeli | Introduction of n-electron valence states for multireference perturbation theory[END_REF] .

We have in a number of recent papers employed a different strategy, namely a hybrid between MCSCF and DFT that capitalizes on the fairly efficient treatment of the (shortrange) dynamical correlation by semi-local DFT approximations and the ability of CASSCF or more general MCSCF models to recover large parts of the static correlation. A number of different CAS-DFT hybrid methods have been suggested [START_REF] Grimme | A combination of Kohn-Sham Density Functional Theory and multi-reference Configuration-Interaction methods[END_REF][START_REF] Marian | Performance of the Density Functional Theory/Multireference Configuration-Interaction method on electronic excitation of extended π-Systems[END_REF][START_REF] Li Manni | Multiconfiguration Pair-Density Functional Theory[END_REF][START_REF] Gagliardi | Multiconfiguration Pair-Density Functional Theory: A New Way To Treat Strongly Correlated Systems[END_REF][START_REF] Savin | Density functionals for the Yukawa electron-electron interaction[END_REF][START_REF] Leininger | Combining long-range Configuration Interaction with short-range density functionals[END_REF][START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF][START_REF] Sharkas | A multiconfigurational hybrid density-functional theory[END_REF][START_REF] Fromager | Multi-Configuration Time-Dependent Density-Functional Theory based on range separation[END_REF][START_REF] Hedegård | Assessment of charge-transfer excitations with time-dependent, range-separated density functional theory based on long-range MP2 and multiconfigurational self-consistent field wave functions[END_REF] , many of which are under active development. There is one important difference between these methods: some add DFT after initial optimization of the multiconfigurational wave function [START_REF] Li Manni | Multiconfiguration Pair-Density Functional Theory[END_REF][START_REF] Gagliardi | Multiconfiguration Pair-Density Functional Theory: A New Way To Treat Strongly Correlated Systems[END_REF] , whereas others optimize the DFT and wave function parts simultaneously [START_REF] Grimme | A combination of Kohn-Sham Density Functional Theory and multi-reference Configuration-Interaction methods[END_REF][START_REF] Marian | Performance of the Density Functional Theory/Multireference Configuration-Interaction method on electronic excitation of extended π-Systems[END_REF][START_REF] Savin | Density functionals for the Yukawa electron-electron interaction[END_REF][START_REF] Leininger | Combining long-range Configuration Interaction with short-range density functionals[END_REF][START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF][START_REF] Sharkas | A multiconfigurational hybrid density-functional theory[END_REF][START_REF] Fromager | Multi-Configuration Time-Dependent Density-Functional Theory based on range separation[END_REF][START_REF] Hedegård | Assessment of charge-transfer excitations with time-dependent, range-separated density functional theory based on long-range MP2 and multiconfigurational self-consistent field wave functions[END_REF] . The method reported here belongs to the latter kind. A challenge for all CAS-DFT hybrids is to avoid double counting of electron correlation, as the CAS wave functions will invariably include some of the dynamical correlation. The method that we discuss in this paper, and that we have been developing since several years [START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF][START_REF] Fromager | Multi-Configuration Time-Dependent Density-Functional Theory based on range separation[END_REF][START_REF] Hedegård | Assessment of charge-transfer excitations with time-dependent, range-separated density functional theory based on long-range MP2 and multiconfigurational self-consistent field wave functions[END_REF][START_REF] Hedegård | Polarizable embedding with a multiconfiguration short-range density functional theory linear response method[END_REF][START_REF] Hedegård | Density matrix renormalization group with efficient dynamical electron correlation through range separation[END_REF][START_REF] Hedegård | Assessment of oscillator strengths with multiconfigurational short-range density functional theory for electronic excitations in organic molecules[END_REF][START_REF] Hedegård | Polarizable embedding density matrix renormalization group[END_REF][START_REF] Hubert | Investigation of Multiconfigurational Short-Range Density Functional Theory for Electronic Excitations in Organic Molecules[END_REF][START_REF] Hubert | Excitation Spectra of Nucleobases with Multiconfigurational Density Functional Theory[END_REF] , stringently avoids this double counting by means of a range separation of the two-electron repulsion operator [START_REF] Savin | Density functionals for the Yukawa electron-electron interaction[END_REF][START_REF] Savin | On degeneracy, near degeneracy and density functional theory[END_REF] . The long-range part is then described by WFT, whereas the short-range part is described by a tailored short-range DFT functional. [START_REF] Toulouse | Short-Range Exchange-Correlation Energy of a Uniform Electron Gas with Modified ElectronElectron Interaction[END_REF] When employing an MC(SCF) wave function as the long-range component, we have dubbed this method MC-srDFT. For the common choice of CAS for the MCSCF part, this becomes CAS-srDFT. However, it should be emphasized, that the method of range separation also allows for other wave-function ansätze: CI-srDFT [START_REF] Leininger | Combining long-range Configuration Interaction with short-range density functionals[END_REF][START_REF] Pedersen | Description of Correlation and Relativistic Effects in Calculations of Molecular Properties[END_REF] , MP2-srDFT [START_REF] Ángyán | van der Waals forces in Density Functional Theory: Perturbational long-range electron-interaction corrections[END_REF][START_REF] Fromager | Self-consistent many-body perturbation theory in range-separated Density-Functional Theory: A one-electron reduced-density-matrix-based formulation[END_REF] , RPA-srDFT [START_REF] Toulouse | Adiabatic-connection fluctuation-dissipation density-functional theory based on range separation[END_REF]43 , CC-srDFT [START_REF] Goll | A short-range gradient-corrected density functional in long-range coupled-cluster calculations for rare gas dimers[END_REF] , RAS-srDFT, NEVPT2-srDFT [START_REF] Fromager | Merging multireference perturbation and Density-Functional Theories by means of range separation: Potential curves for Be 2 , Mg 2 , and Ca 2[END_REF] , and DMRG-srDFT [START_REF] Hedegård | Density matrix renormalization group with efficient dynamical electron correlation through range separation[END_REF] methods have all been implemented.

Our MC-srDFT implementation in a development version of DALTON [START_REF]DALTON, a molecular electronic-structure program[END_REF][START_REF] Aidas | The Dalton quantum chemistry program system[END_REF] employs a general implementation for the long-range MCSCF wave function that allows any spin multiplicity. However, the implementation has until now been restricted to short-range functionals depending only on the total density and not on the spin density. Calculations have accordingly been done on closed-shell systems only. It should in this regard be noted that introducing a dependence on the spin density in the DFT part is not necessary from a theoretical point of view where the energy is solely defined from the total density. Yet, approximate semi-local functionals must depend on the spin density to be accurate for openshell systems, and a semi-local srDFT implementation based solely on the total density in practice excludes useful calculations on many open-shell target systems where we expect MC-srDFT to be able to perform well. For example, many transition metals form openshell complexes, and these are notorious for displaying multireference character.

Since the original formulation of MC-srDFT [START_REF] Pedersen | Description of Correlation and Relativistic Effects in Calculations of Molecular Properties[END_REF] , several srDFT functionals have been generalized to include the spin density. [START_REF] Paziani | Local-spin-density functional for multideterminant density functional theory[END_REF][START_REF] Goll | A short-range gradient-corrected spin density functional in combination with long-range coupled-cluster methods: Application to alkali-metal rare-gas dimers[END_REF] In this paper, we extend the MC-srDFT method with these spin-dependent functionals, both for a local-density approximation (LDA) and a generalized-gradient approximation (GGA). This enables meaningful calculations for openshell systems, which in particular will be beneficial for our endeavors to describe transitionmetal complexes with CAS-srDFT or GAS-srDFT as a computationally cheaper alternative to CASPT2, RASPT2, and NEVPT2.

As a first test of our novel implementation we investigate two prototypical open-shell molecules: the O 2 molecule and the transition metal complex [Fe(H 2 O) 6 ] 3+ . We focus here on rather simple systems which allow for a thorough investigation, and where literature data are available. For both the 3 Σ - g ground-state of O 2 and the spin-state splitting in [Fe(H 2 O) 6 ] 3+ we also carry out a small study on the dependence of CAS-srDFT on the range-separation parameter µ, similarly to what has previously been done for a number of small, closed-shell p-block molecules and atoms. [START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF] It should be emphasized that although O 2 is a prototypical open-shell system for which a wealth of literature data is available [START_REF] Schweitzer | Physical Mechanisms of Generation and Deactivation of Singlet Oxygen[END_REF] , the spin-state splitting between the 3 Σ - g ground-state and the first singlet state, 1 ∆ g , is surprisingly demanding in terms of the theoretical treatment. [START_REF] Klotzchristel | Calculation of Spin-Forbidden Radiative Transitions Using Correlated Wave Funnctions: Lifetimes of b 1 Σ + , a 1 ∆ States in O 2 , S 2 and SO[END_REF][START_REF] Gadzhiev | Structure, Energy, and Vibrational Frequencies of Oxygen Allotropes O n (n 6) in the Covalently Bound and van der Waals Forms: Ab Initio Study at the CCSD (T) Level[END_REF][START_REF] Zen | Properties of reactive oxygen species by quantum Monte Carlo Properties of reactive oxygen species by quantum Monte Carlo[END_REF] Likewise, the spin-state splitting between the lowest electronic states, [START_REF] Chan | The Density Matrix Renormalization Group in Quantum Chemistry[END_REF] A g and 4 T 1g (assuming idealized, octahedral symmetry) in [Fe(H 2 O) 6 ] 3+ has turned out to be challenging; there is currently no consensus on which method is the most accurate, as both CASPT2 and CCSD(T) have been shown to overestimate the spin-state splitting somewhat. [START_REF] Ghosh | High-level ab initio calculations on the energetics of low-lying spin states of biologically relevant transition metal complexes: A first progress report[END_REF] Accordingly, both O 2 and [Fe(H 2 O) 6 ] 3+ have low-lying excited states with different spin multiplicity than the ground-state, and for both a high-level method is evidently needed to obtain correct results.

The paper is organized as follows. In Section II we describe the necessary theory for the extension of MC-srDFT to srDFT functionals that include the spin density as well as the gradient of the spin density, and in Section III we describe implementation details. The theory and implementation details for the closed-shell singlet case is obviously a special case, which previously only has been published in the PhD dissertation by J. K. Pedersen. [START_REF] Pedersen | Description of Correlation and Relativistic Effects in Calculations of Molecular Properties[END_REF] Next, we provide the computational details that were employed for the test calculations (Section IV), and in Section V we discuss the results. Finally, we give a brief conclusion and outlook in Section VI.

II. RANGE-SEPARATED MULTICONFIGURATIONAL DENSITY FUNCTIONAL THEORY

The MC-srDFT method employs a range separation of the two-electron operator [START_REF] Savin | Density functionals for the Yukawa electron-electron interaction[END_REF][START_REF] Savin | On degeneracy, near degeneracy and density functional theory[END_REF] ĝ(1, 2) = ĝlr (1, 2) + ĝsr (1, 2), (1) in which the exact definition of ĝlr (1, 2) and ĝsr (1, 2) can differ [START_REF] Toulouse | Long-range-short-range separation of the electron-electron interaction in Density-Functional Theory[END_REF] . In this work we exclusively use the error function for the range separation:

ĝlr (1, 2) = erf(µr 12 ) r 12 and ĝsr (1, 2) = 1 -erf(µr 12 ) r 12 , (2)
where r 12 = |r 1 -r 2 | and µ is the range-separation parameter, given in reciprocal bohr.

Eq. ( 1) is given in atomic units which we will employ throughout this paper. The energy of a range-separated WFT-DFT hybrid is then given by

E(λ) = E lr (λ) + E sr-H [ρ C ] + E sr-xc [ρ C , ρ S ] = Ψ lr (λ)| Ĥlr |Ψ lr (λ) + E sr-H [ρ C ] + E sr-xc [ρ C , ρ S ], (3)
where λ are the wave-function parameters, to be defined below. We denote the electron charge density as ρ C and the spin density as ρ S , and in order not to overload the notation we let it be implicit that they depend on the wave-function parameters λ:

ρ C ≡ ρ C (r, λ)
and ρ S ≡ ρ S (r, λ). We define the short-range Hartree and exchange-correlation functionals,

E sr-H [ρ C ] and E sr-xc [ρ C , ρ S ],
with the standard one-electron integrals h pq and nuclei-nuclei interaction energy V nn , but with the two-electron Coulomb integrals replaced by integrals with the modified long-range interaction in Eq. ( 2),

g lr pq,rs = φ p φ r |ĝ lr (1, 2)|φ q φ s , (5)
in an orthonormal basis of spatial orbitals {φ p }. The second-quantized operators in Eq. ( 4)

have their usual meaning: Êpq = â † pα âqα + â † pβ âqβ is the singlet excitation operator and êpq,rs = Êpq Êrs -δ qr Êps . The indices p, q, r, s denote general orbitals. For the MC(SCF) wave function parameters we use the notation λ T = (d T , κ T ) where the row vectors d T = {d j } and κ T = {κ pq } designate the configuration coefficients and the orbital rotation parameters, respectively. The real-valued MC wave function is parameterized as

|Ψ lr (λ) = e -κ   |0 + P| d 1 + d | P| d   , (6)
where |0 denotes a normalized reference state

|0 = j c j |j , (7)
while

|d = j d j |j , (8)
is a configuration correction and P = 1 -|0 0| is the projection operator onto the complement to the reference state |0 . The κ operator in Eq. ( 6) is the usual antisymmetric real singlet orbital-rotation operator

κ = pq κ pq Êpq = p>q κ pq Êpq -Êqp ≡ p>q κ pq Ê- pq . (9)
The charge-density-and spin-density-dependent terms in Eq. ( 3) can now be expressed in a second-quantization formulation [START_REF] Saue | Four-component relativistic Kohn-Sham theory[END_REF][START_REF] Salek | Density-Functional Theory of linear and nonlinear time-dependent molecular properties[END_REF] in terms of their associated density operators

ρX (r) = pq Ω pq (r) ÔX pq , (10)
with Ω pq (r) = φ * p (r)φ q (r). We have in Eq. ( 10) introduced a nomenclature that will turn convenient in the next sections for equations which are otherwise identical for charge-density and spin-density operators. In this nomenclature, X = C is used for the regular chargedensity operator ( ÔC pq ≡ Êpq ), while X = S denotes the spin-density operator

ÔS pq ≡ Tpq = â † pα âqα -â † pβ âqβ , (11)
which is required to describe spin-density effects (in a non-relativistic framework). We note that for spin-restricted models as considered here, the charge-density operator is a singlet operator and the spin-density operator is a triplet operator of M S = 0 type. The wave function can be of any spin symmetry S, however, as in Kohn-Sham DFT, in practice the wave function should correspond to a spin component M S = S or M S = -S to give an appropriate spin density for the approximate spin-polarized exchange-correlation functional.

The electron charge density ρ C (r, λ) and electron spin density ρ S (r, λ) are obtained as expectation values according to

ρ X (r, λ) = Ψ lr (λ)|ρ X (r)|Ψ lr (λ) = pq Ω pq (r)D X pq (λ), (12)
where D X pq (λ) is the (p, q) element of the one-electron reduced (spin-)density matrix

D X pq (λ) = Ψ lr (λ)| ÔX pq |Ψ lr (λ) . (13)
The E sr-H [ρ C ] and E sr-xc [ρ C , ρ S ] terms in Eq. ( 3) can now be defined. The former is the short-range Hartree energy which depends only on the total density matrix

E sr-H [ρ C ] = 1 2 pq,rs D C pq (λ) g sr pqrs D C rs (λ) ≡ 1 2 pq D C pq (λ) j sr pq (λ), (14)
where g sr pqrs are the short-range two-electron integrals (defined with the short-range interaction in Eq. ( 2)). The final term, E sr-xc [ρ C , ρ S ], in Eq. ( 3) is the short-range exchangecorrelation functional. This term has an explicit dependence on both the charge and spin densities

E sr-xc [ρ C , ρ S ] = e sr-xc (ρ C (r, λ), ρ S (r, λ))dr, (15)
in accordance with Ref. 56, where e sr-xc (ρ C (r, λ), ρ S (r, λ)) is the short-range exchangecorrelation energy density. It should be noted that the expression in Eq. ( 15) assumes the local-density approximation of the short-range exchange-correlation functional (srLDA) [START_REF] Paziani | Local-spin-density functional for multideterminant density functional theory[END_REF] .

For brevity we only consider srLDA here, but we have also implemented the code for functionals based on a generalized-gradient approximation (GGA). The spin-dependent GGAs additionally depend on the electron charge density gradient and the electron spin density gradient. The additional terms are structurally similar to the srLDA terms and will be given in Appendix A.

The MC-srDFT wave function is optimized with the restricted-step second-order MCSCF optimization algorithm [START_REF] Aa | A direct approach to second order MCSCF calculations using a Norm-Extended Optimization scheme[END_REF][START_REF] Aa | MCSCF optimization using the direct, restricted step second-order norm-extended optimization method[END_REF][START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF][START_REF] Aa | Efficient optimization of large scale MCSCF wave functions with a restricted-step algorithm[END_REF][START_REF] Aa | Electron correlation in molecules using direct second order MCSCF[END_REF] as implemented in DALTON. [START_REF]DALTON, a molecular electronic-structure program[END_REF][START_REF] Aidas | The Dalton quantum chemistry program system[END_REF] This algorithm is based on a second-order Taylor expansion of the electronic energy in the wave-function parameters, λ,

around λ = 0 E(λ) = E 0 + g T λ + 1 2 λ T Hλ + • • • . (16)
The electronic gradient, g, and electronic Hessian, H, are blocked according to the configurational and orbital parameters of the wave function. Thus the gradient reads

g =   g c g o   , (17)
and each gradient element (configurational or orbital) has both WFT (lr) and DFT (sr) contributions

g i = g lr,i + g sr-H,i + g sr-xc,i = ∂E lr ∂λ i + ∂E sr-H [ρ C ] ∂λ i + ∂E sr-xc [ρ C , ρ S ] ∂λ i . ( 18
)
The electronic Hessian in Eq. ( 16) is evaluated as 62 H = P KP , where P denotes the matrix representation of the P operator in Eq. ( 6) and K has WFT and DFT contributions

K ij = K lr,ij + K sr-H,ij + K sr-xc,ij = ∂ 2 E lr ∂λ i ∂λ j + ∂ 2 E sr-H [ρ C ] ∂λ i ∂λ j + ∂ 2 E sr-xc [ρ C , ρ S ] ∂λ i ∂λ j . (19)
In the actual implementation the Hessian matrix is never constructed explicitly. Instead, Hessian contributions to the wave-function optimization process are obtained in a direct fashion based on trial vectors

σ n = P KP b n = P K b n = P   K cc K co K oc K oo     b c n b o n   , (20)
where the individual contributions to σ n can be written in terms of modified Fock matrices (see Section III) and the second equal sign stems from that we require each trial vector to A. The sr-DFT contributions to the electronic gradient

The short-range exchange-correlation contributions to the gradient in Eq. ( 18) can be obtained as

g sr-xc,i = ∂E sr-xc [ρ C , ρ S ] ∂λ i = ∂e sr-xc (ρ C (r), ρ S (r)) ∂λ i dr. ( 21
)
The derivative within the kernel in Eq. ( 21) is obtained through the chain rule

∂e sr-xc (ρ C (r), ρ S (r)) ∂λ i = ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ C ∂ρ C (r) ∂λ i + ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ S ∂ρ S (r) ∂λ i = pq ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ C Ω pq (r) ∂D C pq (λ) ∂λ i + ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ S Ω pq (r) ∂D S pq (λ) ∂λ i , (22)
where ρ C (r) and ρ S (r) from Eq. ( 12) have been inserted. As usual in MCSCF schemes, the orbital and configurational parts of the gradient leads to different computational expressions.

The orbital part of the short-range exchange-correlation gradient becomes

∂E sr-xc [ρ C , ρ S ] ∂κ rs = g C sr-xc,rs + g S sr-xc,rs , (23)
where the two gradient terms are

g C sr-xc,rs = 0| [ Ê- rs , V C,g sr-xc ] |0 = 2 0| [ Êrs , V C,g sr-xc ] |0 , (24a)
g S sr-xc,rs = 0| [ Ê- rs , V S,g sr-xc ] |0 = 2 0| [ Êrs , V S,g sr-xc ] |0 , (24b)
both defined in terms of an effective operator (X = C and X = S, respectively)

V X,g sr-xc = pq ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ X Ω pq (r)dr ÔX pq ≡ pq V X,g sr-xc,pq ÔX pq . ( 25
)
We note that the g C sr-xc,rs term in Eq. ( 24a) is the term already present in the closed-shell formalism, although in that case the functional kernel in the effective operator does not depend on ρ S . The g S sr-xc,rs term in Eq. (24b) exclusively occurs in the open-shell formalism. The two gradient terms are commonly denoted g X sr-xc . The configurational part of the gradient can be defined similarly

∂E sr-xc [ρ C , ρ S ] ∂d j = g C sr-xc,j + g S sr-xc,j , (26)
where

g X sr-xc,j = 2 j| V X,g sr-xc |0 -c j 0| V X,g sr-xc |0 , (27)
and we again have employed the effective operators defined in Eq. ( 25).

Finally, we also briefly describe the gradient contribution from the short-range Hartree term, E sr-H [ρ C ], in Eq. ( 14)

g sr-H,i = ∂E sr-H [ρ C ] ∂λ i = pq j sr pq (λ) ∂D C pq (λ) ∂λ i . ( 28
)
The configurational and orbital gradient contributions are identical to Eqs. (24a) and ( 27), except that the effective operator

V g sr-H = pq j sr pq Êpq , (29)
replaces V C,g sr-xc . Combined, we obtain the effective operators for the total short-range Hartreeexchange-correlation functional as

V X,g sr-Hxc =    V g sr-H + V C,g sr-xc , X = C V S,g sr-xc , X = S (30)
The total gradient vector (from the short-range DFT part) can thus be obtained from

g sr-Hxc = g C sr-Hxc + g S sr-xc . (31)

B. The sr-DFT contributions to the electronic Hessian sigma vectors

As for the gradient in the previous subsection, we focus on the exchange-correlation term since the second derivative of E lr in Eq. ( 3) is equivalent to the terms from regular MCSCF [START_REF] Aa | A direct approach to second order MCSCF calculations using a Norm-Extended Optimization scheme[END_REF][START_REF] Aa | MCSCF optimization using the direct, restricted step second-order norm-extended optimization method[END_REF][START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF][START_REF] Aa | Efficient optimization of large scale MCSCF wave functions with a restricted-step algorithm[END_REF] , and the Hartree term can be obtained from the total density matrix and has no explicit spin-dependent terms. The Hessian contributions from the exchange-correlation functional are given by

K sr-xc,ij = ∂ 2 E sr-xc [ρ C , ρ S ] ∂λ i ∂λ j = ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂λ i ∂λ j dr, (32)
and from the chain rule

∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂λ i ∂λ j = ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ 2 C ∂ρ C (r) ∂λ i ∂ρ C (r) ∂λ j + ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ C ∂ρ S ∂ρ S (r) ∂λ i ∂ρ C (r) ∂λ j + ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ S ∂ρ C ∂ρ C (r) ∂λ i ∂ρ S (r) ∂λ j + ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ 2 S ∂ρ S (r) ∂λ i ∂ρ S (r) ∂λ j + ∂e xc ρ C (r), ρ S (r)) ∂ρ C ∂ 2 ρ C (r) ∂λ i ∂λ j + ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ S ∂ 2 ρ S (r) ∂λ i ∂λ j . (33)
In practice we utilize a direct Hessian technique [cf. Eq. ( 20)] where the quantity in Eq. ( 32)

is contracted with configurational (b c j ) or orbital (b o pq ) trial vectors. The direct Hessian contributions from the E sr-xc [ρ C , ρ S ] term thus consist of four different types j ∂ 2 E sr-xc [ρ C , ρ S ] ∂d i ∂d j b c j = σ C;c sr-xc,i + σ S;c sr-xc,i (34a) j ∂ 2 E sr-xc [ρ C , ρ S ] ∂κ pq ∂d j b c j = σ C;c sr-xc,pq + σ S;c sr-xc,pq (34b)
r>s

∂ 2 E sr-xc [ρ C , ρ S ] ∂d i ∂κ rs b o rs = σ C;o sr-xc,i + σ S;o sr-xc,i (34c)
r>s

∂ 2 E sr-xc [ρ C , ρ S ] ∂κ pq ∂κ rs b o rs = σ C;o sr-xc,pq + σ S;o sr-xc,pq . (34d)
The expressions in Eqs. (34a)-(34d) are obtained by inserting ρ C (r) and ρ S (r) from Eq. ( 12)

into Eq. ( 33), and inserting the resulting expression into Eq. ( 32) and contracting with configurational or orbital trial vectors, respectively. We sketch the two steps below; the insertion of the densities gives

∂ 2 E sr-xc [ρ C , ρ S ] ∂λ i ∂λ j = X,Y =C,S pq,rs ∂D Y rs (λ) ∂λ j ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ Y ∂ρ X Ω pq (r)Ω rs (r)dr ∂D X pq (λ) ∂λ i + X=C,S pq ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ X Ω pq (r)dr ∂ 2 D X pq (λ) ∂λ i ∂λ j . ( 35
)
By contraction of the above equation with the trial vectors, the configurational part becomes

σ X;c sr-xc,i = Y =C,S pq,rs ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ Y ∂ρ X j ∂D Y rs (λ) ∂d j b c j Ω rs (r) Ω pq (r)dr ∂D X pq (λ) ∂λ i + pq ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ X Ω pq (r)dr j ∂ 2 D X pq (λ) ∂λ i ∂d j b c j , (36)
while the orbital part becomes

σ X;o sr-xc,i = Y =C,S pq,rs ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ Y ∂ρ X t>u ∂D Y rs (λ) ∂κ tu b o tu Ω rs (r) Ω pq (r)dr ∂D X pq (λ) ∂λ i + pq ∂e sr-xc (ρ C (r), ρ S (r)) ∂ρ X Ω pq (r)dr t>u ∂ 2 D X pq (λ) ∂λ i ∂κ tu b o tu . (37)
The first term in both Eqs. ( 36) and ( 37) occur due to the non-linearity of the srDFT energy functional. These terms can be formulated as transformed gradient terms with effective operators as summarized below. The last term in both Eqs. (36) and (37) are similar to the one-electron Hessian part for a regular MCSCF calculation, the only difference being that they employ different integrals (over the first-order derivative of the srDFT functional). To define the effective operators employed in the first term of the two equations, we define the transition density matrices for configurational trial vectors

D Y (1c) rs ≡ j ∂D Y rs (λ) ∂d j b c j = j b c j j| ÔY rs |0 + 0| ÔY rs |j b c j = B| ÔY rs |0 + 0| ÔY rs |B , ( 38
)
where |B is the state vector generated by the current trial vector: |B = j b c j |j . For the orbital trial vectors, we define the one-index transformed density matrix

D Y (1o) rs ≡ t>u ∂D Y rs (λ) ∂κ tu b o tu = t D Y ts b o rt + D Y rt b o st . (39)
The transformed configurational and orbital (spin-)density matrices, D Y (1c) and D Y (1o) , are collectively denoted by D Y (1λ) . We can now define the integrals that comprise the effective operators in terms of the linearly transformed densities

V X[1λ] sr-xc,pq [D Y (1λ) ] = ∂ 2 e sr-xc (ρ C (r), ρ S (r)) ∂ρ Y ∂ρ X rs D Y (1λ) rs Ω rs (r) Ω pq (r)dr, (40)
which is a generic term for the integrals in the curvy brackets of the first term in Eqs. ( 36) and (37). With this generic definition of the integrals, we can define the two effective operators of either singlet or triplet types

V C[1λ] sr-xc = pq V C[1λ] sr-xc,pq [D C(1λ) ] + V C[1λ] sr-xc,pq [D S(1λ) ] Êpq (41a) V S[1λ] sr-xc = pq V S[1λ] sr-xc,pq [D C(1λ) ] + V S[1λ] sr-xc,pq [D S(1λ) ] Tpq . ( 41b
)
The direct Hessian terms in Eqs. (34a)-(34d) can now be written in the generic form

σ X;c sr-xc,i = 2 i| V X[1c] sr-xc |0 -c i 0| V X[1c] sr-xc |0 + 2 i| V X,g sr-xc |B -0| V X,g sr-xc |0 b c j (42a) σ X;c sr-xc,pq = 2 0| [ Êpq , V X[1c] sr-xc ] |0 + 2 0|[ Ê- pq , V X,g sr-xc ]|B (42b) σ X;o sr-xc,i = 2 i| V X[1o] sr-xc |0 -c i 0| V X[1o] sr-xc |0 + 2 i| V X,g sr-xc |0 (42c) σ X;o sr-xc,pq = gX sr-xc,pq + 0|[ Ê- pq , V X[1o] sr-xc + V X,g sr-xc ]|0 + 1 2 t g X sr-xc,tp b o qt -g X sr-xc,tq b o pt . (42d)
The orbital gradient elements, g X sr-xc,pq , were defined in Eqs. (24a)-(24b), the V X,g sr-xc operators were given in Eq. ( 25), and V X,g sr-xc is the one-index transformed form of the V X,g sr-xc operator in Eq. ( 25). The one-index transformed operators employ transformed integrals [START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF][START_REF] Aa | Electron correlation in molecules using direct second order MCSCF[END_REF] , e.g., V X,g sr-xc is equivalent to Eq. ( 25), but with the V X,g sr-xc,pq integrals replaced by

Ṽ X,g sr-xc,pq = t V X,g sr-xc,tq b o pt + V X,g sr-xc,pt b o qt . (43)
Similarly, the expression for the transformed gradient elements, gX sr-xc,pq , is equivalent to Eqs. (24a)-(24b), but with V X,g sr-xc replaced with V X,g sr-xc . The terms with V X,g sr-xc and V X,g sr-xc operators arise from the last terms of Eqs. (36) and (37), cf. the part within curly brackets. For these terms, the regular Hessian structure can be discerned.

At this point, it is noted that the σ S;c and σ S;o terms in Eqs. (42a)-(42d) enter exclusively for open-shell systems, and are thus part of the extension for this work. Accordingly, with respect to the Hessian this work concerns the implementation of the effective operator of triplet type, i.e. Eq. (41b), employed to define the σ S;c and σ S;o terms. In addition, the term that depends on D S(1λ) in Eq. (41a) also enters exclusively for open-shell systems and has been implemented here.

In addition to the contributions from the exchange-correlation functional, the direct Hessian terms also contain contributions from the short-range Hartree term, E sr-H [ρ C ], in Eq. ( 14)

∂ 2 E sr-H [ρ C ] ∂λ j ∂λ i = pq,rs ∂D C rs (λ) ∂λ j g sr pqrs ∂D C pq (λ) ∂λ i + pq,rs D C rs (λ)g sr pqrs ∂ 2 D C pq (λ) ∂λ j ∂λ i . ( 44
)
As for the short-range exchange-correlation contributions, the direct Hessian short-range Hartree contributions can be obtained by contraction of configurational and orbital trial vectors, respectively. The terms have the same structure as Eqs. (42a)-(42d), but without any spin-dependence, and can be defined employing the effective operator

V C[1λ],c sr-H = pq j C(1λ) sr,pq [D C(1λ) ] Êpq with j C(1λ) sr,pq [D C(1λ) ] = rs D C(1λ) rs g sr rspq . (45)
We can combine the operators utilized for the short-range Hartree and short-range exchangecorrelation terms in a common, effective operator

V X[1λ] sr-Hxc =    V C[1λ] sr-H + V C[1λ] sr-xc , X = C V S[1λ] sr-xc , X = S , (46)
which replaces V X[1λ] sr-xc in Eqs. (42a)-(42c).

III. IMPLEMENTATION

The configurational and orbital gradients for MC-srDFT were given in Eqs. ( 23)-( 31) and depend on V C,g sr-xc and V S,g sr-xc operators [Eq. (25)]. The implementations of the equations for the gradient and (direct) Hessian in regular MCSCF within DALTON has been described previously. [START_REF] Aa | A direct approach to second order MCSCF calculations using a Norm-Extended Optimization scheme[END_REF][START_REF] Aa | MCSCF optimization using the direct, restricted step second-order norm-extended optimization method[END_REF][START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF] A central part of this implementation builds on the generation of generalized Fock-type matrices

f pq = r D C pr h qr + rst P pr,st g qr,st , (47)
where P is the reduced two-electron density matrix with elements P pr,st = Ψ|ê pr,st |Ψ . In the following, we will in addition to the general orbitals (with indices p, q, r, s) need to denote both inactive and active orbitals, for which we use the indices i, j and u, v, x, y, respectively.

In terms of the generalized Fock matrix, the orbital gradient becomes

g o rs = 2(f rs -f sr ). ( 48
)
The direct Hessian is constructed in terms of modified gradient expressions, which in Ref. 60 are denoted transition and one-index transformed gradient expressions. These can again be constructed from Eq. ( 47) by employing a one-index transformed density matrix or a transition density matrix instead of the reference density matrix. [START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF][START_REF] Siegbahn | A new direct CI method for large CI expansions in a small orbital space[END_REF] The only non-zero long-range contributions to the Fock matrix elements are 60

f lr iq = 2(f I,lr iq + f A,lr iq ) ( 49
)
f lr vq = u D C uv f I,lr uq + Q lr vq , (50)
where Q lr vq = uxy P vu,xy g lr qu,xy . The Fock matrices f I,lr and f A,lr are called the inactive and active Fock matrix, respectively. They are defined as in Eqs. (3.11) and (3.12) of Ref. 60, but with long-range two-electron integrals. For an MC-srDFT description, the generalized Fock matrix elements must be altered to accommodate the potential contributions from the short-range DFT functional. Eqs. ( 49) and ( 50) become

f iq = f lr iq + V C,g sr-Hxc,iq + V S,g sr-xc,iq (51)
f vq = u D c uv f I,lr uq + V C,g sr-Hxc,uq + V S,g sr-xc,uq + Q lr vq , (52)
where V X,g sr-Hxc,pq are the matrix elements of the operator in Eq. (30). Depending on the employed (spin-)density matrix, Eqs. ( 51) and ( 52) can be either a regular (spin-)Fockmatrix, a transition-(spin-)Fock matrix or an one-index transformed (spin-)Fock matrix.

The auxiliary routines that construct Fock matrices were modified to accommodate shortrange functionals in connection with previous work [START_REF] Pedersen | Description of Correlation and Relativistic Effects in Calculations of Molecular Properties[END_REF] . These auxiliary routines have been extended with the construction of the spin-density matrix required for the spin-dependent terms of V X,g sr-Hxc,pq . The direct Hessian [Eqs. (42a)-(42d)] is comprised of terms that include D X(1λ) (linear transformed terms) and "regular" Hessian terms. The former terms required slight modification for this work, whereas the latter terms are automatically included after the above-mentioned modifications of the gradient routines.

Finally, we have added code for the spin-dependent short-range exchange-correlation functionals from Ref. 48 and the spin-dependent gradient-corrected functional of Ref. 49 (see next section). This includes code for the numerical evaluation of the gradient and Hessian of the short-range functional kernel with respect to the spin density.

IV. COMPUTATIONAL DETAILS

The theory presented in the previous sections was implemented in a development version of the DALTON 2016 program. [START_REF]DALTON, a molecular electronic-structure program[END_REF][START_REF] Aidas | The Dalton quantum chemistry program system[END_REF] All CAS-srDFT calculations were carried out with this new implementation. We have employed CAS-srDFT for two molecules, namely O For O 2 we focus on the 3 Σ - g and 1 ∆ g lowest-energy states. The calculations on O 2 were carried out at the experimental ground-state equilibrium bond distance (R eq = 1.207 Å) and at a stretched bond distance (R stretch = 2.0 Å) with the cc-pVTZ basis set. [START_REF] Dunning | Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen[END_REF][START_REF] Dunning | Gaussian Basis Functions for Use in Molecular Calculations. I. Contraction of (9s5p) Atomic Basis Sets for the First-Row Atoms[END_REF] For O 2 , we have employed three CAS(n, m) spaces with n electrons in m orbitals; CAS (8,6), CAS (12,8) and CAS (12,16). These spaces correspond to including the orbitals 3σ

2 g 1π 4 u 1π * 2 g 3σ * 0 u , 2σ 2 g 2σ * 2 u 3σ 2 g 1π 4 u 1π * 2 g 3σ * 0 , and 2σ 2 g 2σ * 2 u 3σ 2 g 1π 4 u 1π * 2 g 3σ * 0 u 4σ 0 g 4σ * 0 u 5σ 0 g 2π 0 u 2π * 0 g 5σ * 0
u , respectively, where the occupations are for the reference configuration.

For [Fe(H 2 O) 6 ] 3+ we focus on the 6 A g and 4 T 1g lowest-energy states (using a nomenclature that assumes octahedral symmetry), and we used the structure from Ref. 66, which was optimized for the 6 A g ground state within a cluster including the second solvation shell. We have removed this second solvation shell, and our calculations only concern the [Fe(H 2 O) 6 ] 3+ core. Note that the CASPT2 value from Ref. 66 to which we compare was obtained exactly this way (cf. Table 7 of Ref. 66). It was shown that such a calculation includes a large part of the solvation effect. For this complex, we employed the def2-TZVPP basis set [START_REF] Weigend | Balanced basis sets of split valence, triple zeta valence and quadruple zeta valence quality for H to Rn: Design and assessment of accuracy[END_REF] for iron and oxygen, and the def2-SV(P) basis set for hydrogen. Following Ref. 66 we employed an active space that includes the five d orbitals, two additional σ Fe-O bonding orbitals, and a second shell of d orbitals yielding a CAS(9,12) active space. The orbitals comprising the active space are shown in Figure 1 As short-range exchange-correlation functionals we employed the spin-dependent srLDA functional from Ref. 48 and the spin-dependent gradient-corrected functional of Ref. 49, denoted by srPBE. The name of the latter reflects that it originates from the Perdew-Burke-Ernzerhof (PBE) functional [START_REF] Perdew | Generalized Gradient Approximation made simple[END_REF] , and have been extended to the short-range interaction in Ref. 69, and further modified by Goll et al. [START_REF] Goll | A short-range gradient-corrected density functional in long-range coupled-cluster calculations for rare gas dimers[END_REF][START_REF] Goll | A short-range gradient-corrected spin density functional in combination with long-range coupled-cluster methods: Application to alkali-metal rare-gas dimers[END_REF] .

V. DEMONSTRATION OF THE OPEN-SHELL CAS-SRDFT METHOD

We here discuss the results from CAS-srLDA and CAS-srPBE calculations on two prototype open-shell systems.

A. The dioxygen molecule

For dioxygen we start with a discussion of the optimal value of the range-separation parameter µ. This matter has been discussed for CAS-srDFT [START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF][START_REF] Hubert | Investigation of Multiconfigurational Short-Range Density Functional Theory for Electronic Excitations in Organic Molecules[END_REF] and also for other rangeseparated methods. Note that with an exact srDFT functional and a full CI long-range wave function (in a complete basis) the exact result would be obtained independent of the µ value. The µ value determines the mixture between the wave function and srDFT parts, and is thus a parameter that shifts correlation effects between the two parts. When we refer to an "optimal µ value" in the following, we mean the value giving accurate results with the least possible computational effort. As the computational time for the srDFT part is basically independent of µ while the time for accurate long-range wave function contributions grows rapidly with µ, the optimal µ value is the smallest µ value for which the non-local correlations not treated properly by semi-local srDFT functionals are not in the realm of the srDFT part. For multireference models, one recipe has been to compare HF-srDFT and CAS-srDFT. From this, it was found that a value of µ of 0.4 bohr -1 allocates the main part of the dynamical correlation to the srDFT functional, while the long-range wave function still incorporates a substantial part of static correlation in multiconfigurational systems. [START_REF] Fromager | On the universality of the long-/short range separation in multiconfigurational Density-Functional Theory[END_REF] The µ = 0.4 value was recently confirmed for excitation energies obtained from linear-response theory. [START_REF] Hubert | Investigation of Multiconfigurational Short-Range Density Functional Theory for Electronic Excitations in Organic Molecules[END_REF] Here, we investigate the optimal value of µ for the O 2 molecule in its triplet ground state 3 Σ - g . It is expected that this state is dominated by the electronic configuration

|Ψ O 2 Req ( 3 Σ - g ) = C 0 3 |1σ 2 g 1σ * 2 u 2σ 2 g 2σ * 2 u 3σ 2 g 1π 2 u,x 1π 2 u,y 1π * 1 g,x 1π * 1 g,y 3σ * 0 u + • • • ≡ C 0 3 |π 2 u,x π 2 u,y π * 1 g,x π * 1 g,y + • • • (53)
at the equilibrium bond distance. The 3 | . . . notation includes the spin multiplicity of the configuration state function, and we have in the last equality left out the σ electrons for brevity. This dominant configuration can be described by a single-determinant high-spin wave function (M S = 1), which means that HF-srDFT and CAS-srDFT should provide similar energies. In order to have a system that displays multiconfigurational character, and is directly comparable to the above system, we also consider a bond distance corresponding to a stretched bond (R stretch = 2.0 Å). In this case, we expect a splitting of HF-srDFT and CAS-srDFT energies, and the parameter µ should be chosen to reflect this.

The HF-srDFT and CAS-srDFT energies as functions of the value of µ are shown in Figure 2, both for R eq and R stretch . For the equilibrium bond distance, the HF-srDFT and MC-srDFT energies are identical (within 10 -3 hartree) until we reach a value of µ ≈ 0.4.

This value is obtained for both of the employed srDFT functionals, and the result is thus very similar to what was found in Ref. 27. We have tested three different active spaces of increasing size. The smallest CAS (8,6) contains only σ and π orbitals from the oxygen 2p orbitals, whereas the CAS (12,8) active space additionally includes the σ orbitals with origin in the 2s orbitals and thus constitutes the full valence space. Finally, the CAS (12,16) is an extended active space also including the oxygen 3s and 3p orbitals. For R eq the splitting occurs at around the same value of µ (cf. Figure 2, left), independent of the employed active space.

For an elongated bond, the system becomes multiconfigurational as can be seen from the fact that the HF-srDFT and MC-srDFT energies are clearly different, already at small values of µ (cf. We further comment on a few more observations from Figure 2. First, we note that for 0.5 the srPBE functional, the µ dependence is benign for small µ values in case of O 2 . For this small molecule the long-range CASSCF wave function can cover the correlation up to about µ = 0.6 -0.7, as seen from the almost horizontal energy curves in this range of values of µ.

As expected, this is in contrast with the CAS-srLDA energy, which is too high for µ = 0 and therefore decreases with increasing µ, since the long-range wave function with increasing µ can describe an increasing part of the correlation. That this happens before the HF-srLDA and MC-srLDA curves separate must mean that this effect is mostly in the exchange part.

Further, that the CAS-srLDA energies decrease with increasing µ while CAS-srPBE energies increase indicates that pure CASSCF is better than any of the CAS-srLDA models, while CAS-srPBE contains some correlation effects not described by the pure CASSCF.

It is also interesting to study the first singlet excited state 1 ∆ g . In D 2h symmetry this electronic state can either be characterized by a dominant configuration given by a linear combination of the Slater determinants in which the π * g orbitals are doubly occupied, i.e.

|Ψ O 2

Req

( 1 ∆ x 2 -y 2 g ) = C 0 1 √ 2 ( 1 |π 2 u,x π 2 u,y π * 2 g,x π * 0 g,y -1 |π 2 u,x π 2 u,y π * 0 g,x π * 2 g,y ) + • • • (54)
or by a dominant configuration in which the π * g,x and π * g,y orbitals are spin-singlet coupled

|Ψ O 2 Req ( 1 ∆ xy g ) = C 0 1 |π 2 u,x π 2 u,y π * 1 g,x π * 1 g,y + • • • . (55)
The states in Eqs. ( 54) and ( 55) are exactly degenerate in D ∞h but belong to different irreducible representations in the D 2h subgroup: A 1g and B 1g . None of the dominant configurations in these two wave functions can be represented by a spin-restricted single determinant, and thus cannot be obtained in standard spin-restricted Kohn-Sham DFT. By contrast, the long-range wave function of MC-srDFT can describe the 1 ∆ g state in both the D 2h irreducible representations of Eqs. ( 54) and (55). The method correctly predicts that the 1 ∆ x 2 -y 2 g and 1 ∆ xy g states are degenerate. As expected, the dominant configurations in the long-range CAS wave functions are indeed given as the ones in Eqs. ( 54) and (55). In Figure 3 (left), we have reported the value of |C 0 | 2 for the dominant configuration in the long-range wave function as a function of µ.

The 3 Σ - g , 1 ∆ x 2 -y 2 g
, and 1 ∆ xy g states behave remarkably similarly, in agreement with textbook qualitative MO theory, and the same trend is observed: A steadily decreasing value of |C 0 | 2 as we increase µ, reflecting that the longrange wave function must become more multiconfigurational to treat the electron correlation effects moved from the short-range to the long-range parts with increasing µ. The value of |C 0 | 2 for the 1 ∆ states decreases fastest with increasing mixture of long-range wave function.

Yet, at small values of µ (including µ = 0.4) we can expect that the singlet-triplet splitting will be determined to a large degree by the leading configuration for each state, and by how the short-range functionals translate the density obtained from the leading configuration into an energy. Any possible error is therefore likely due to the srDFT functional (as will be further discussed below). As a remark, we note that the 1 ∆ x 2 -y 2 g and 1 ∆ xy g energies as well as their respective |C 0 | 2 values are indistinguishable at all values of µ. This shows that the method treats singlet-coupled open shells qualitatively correct.

With these results in mind, we can also investigate the singlet-triplet splitting ∆E( 3 Σ - g → 1 ∆ g ). We note however, that we cannot expect accurate results due to deficiencies in the employed short-range functionals. In particular, the short-range exchange does not correctly cancel the short-range Hartree self-repulsion for the singly occupied π orbitals in singlet states. In case of O 2 , this can be seen from the above analysis of the dominant configurations: the long-range wave function correctly describes a spin-singlet coupled state in which the total spin density is zero but which locally may become both positive and negative (perhaps as shown most obviously in the 1 ∆ xy g state, Eq. ( 55)). Yet, our present srDFT functionals are not equipped with the ability to employ local spin densities, as required for a spin-singlet coupled state. We therefore expect too high energy of the 1 ∆ g

( 3 Σ - g → 1 ∆ g ).
state and thus of the splitting (as indeed observed, see below). Our investigation of the singlet-triplet splitting therefore focus mainly on the µ-dependence whereas the values itself is only briefly discussed. The dependence on µ of the singlet-triplet splitting is shown in Figure 3 (right) for different active spaces and short-range functionals and given in Table I for for selected values of µ with CAS(12,8)-srDFT and CAS(12,16)-srDFT. In general, the dependence on µ is moderate (at most a few kJ/mol) in the interval µ = 0.1-0.5. The difference is larger (≈ 25-30 kJ/mol) between the results obtained with the srLDA and srPBE functionals; thus there is a significant effect of employing a gradient-corrected functional. Although the srLDA results are in fact closer to the experimental value than the srPBE results, both functionals have the above-mentioned flaws leading to destabilization of the singlet energy. Therefore, the good performance of srLDA must be considered fortuitous. We compare the splitting obtained with srLDA and srPBE to other results from the literature in Table II We now turn to a system where methods such as CCSDTQ are probably beyond reach.

The prototypical [Fe(H 2 O) 6 ] 3+ complex has been described in numerous theoretical works with a variety of methods. We will here focus on the 6 A g → 4 T 2g transition, where in our labeling we assume octahedral symmetry. The higher spin state ( 6 A g ) can qualitatively be described as an open-shell d-complex with five singly occupied d-orbitals, whereas the lower spin state ( 4 T 2g ) can be described as a doubly occupied orbital in the lowest-lying t 2g orbital level and a hole in the e g level. The lower spin-state thus does not involve a coupling between spatially separated orbitals, suggesting that we will not observe errors due to the lack of local spin dependence in the employed srDFT functionals as described in the previous section.

The experimental band for the 6 A g → 4 T 2g transition has a maximum at 12600 cm -1

(1.56 eV) in acidic, aqueous solution [START_REF] Jørgensen | Studies of Absorption Spectra IV. Some New Transition Group Bands of Low Intensity[END_REF] . However, previous results obtained with CASPT2 and CCSD(T) have showed remarkably large deviations (up to 1.5 eV) from this value [START_REF] Ghosh | High-level ab initio calculations on the energetics of low-lying spin states of biologically relevant transition metal complexes: A first progress report[END_REF][START_REF] Yang | Multireference Ab Initio Study of Ligand Field d-d Transitions in Octahedral Transition-Metal Oxide Clusters[END_REF] .

Large deviations were also found with the spectroscopically oriented configuration interaction (SORCI) method [START_REF] Neese | Advanced aspects of ab initio theoretical optical spectroscopy of transition metal complexes: Multiplets, spin-orbit coupling and resonance Raman intensities[END_REF] . On the other hand, both DFT [START_REF] Harris | Structure and Relative Spin-State Energetics of [Fe(H2O)6]3+: A Comparison of UHF, Møller-Plesset, Nonlocal DFT, and Semiempircal INDO/S Calculations[END_REF][START_REF] Hughes | Correcting Systematic Errors in DFT Spin-Splitting Energetics for Transition Metal Complexes[END_REF][77] and the semi-empirical INDO/S This leads to a CAS space of 9 electrons in 12 orbitals, whereas some earlier studies employed the smaller CAS(5,5) and CAS(5,10) active spaces. [START_REF] Yang | Multireference Ab Initio Study of Ligand Field d-d Transitions in Octahedral Transition-Metal Oxide Clusters[END_REF] . The effect was estimated to be around 0.24 eV on the transition energy. However, Radoń et al. [START_REF] Radoń | Spin-State Energetics of Fe(III) and Ru(III) Aqua Complexes: Accurate ab Initio Calculations and Evidence for Huge Solvation Effects[END_REF] also stressed that their study was not directly comparable to some of the previous studies as they (unlike Ref. 73) employed an IPEA-shifted zeroth-order Hamiltonian and also larger basis sets of which the combined effect can be estimated to be around 0.3 eV on the transition energy.

Here, we start by investigating the dependence on µ of the spin-state splitting ∆E(

VI. CONCLUSIONS

In this paper, we have presented the theory and implementation of a generalization of the MC-srDFT method to employ spin-dependent short-range density functionals. We have applied the method to two well-known cases, namely the O 2 molecule which has a triplet ground state, and the transition metal complex [Fe(H 2 O) 6 ] 3+ .

For the O 2 molecule, a wealth of literature data exist for the spin-state splitting between the 3 Σ - g ground state and the lowest-lying singlet state 1 ∆ g , and QMC or CCSDTQ methods provide a value in very good agreement with experiment. The spin-state splitting has here been investigated with a range of values for the range-separation parameter µ, but the splitting is not strongly dependent on µ. Compared to the experimental value, CAS-srPBE lies between CCSD(T) and CCSDT in accuracy. This is already impressive, but there is also room for improvement, as the 1 ∆ g state is estimated too high in energy by both CCSD(T), CCSDT and CAS-srDFT. The former two methods can be improved by higher-order cluster expansions, but this comes at a high computational cost. An improvement (that will be computationally cheap) for CAS-srDFT is to allow dependence of local spin densities in the employed srDFT functionals, which is expected to lower the energy of the spin-coupled 1 ∆ g state.

Next, we investigated the spin-state splitting ∆E( 6 A g → 4 T 1g ) in [Fe(H 2 O) 6 ] 3+ . In this case, CCSD(T) and CASPT2 are both significantly off. A method such as CCSDTQ is too computationally demanding, while CAS-srPBE shows very good agreement with experiment.

It should also be emphasized that most other high-level methods require large basis set expansions with high-order angular momenta to describe the electron-electron Coulomb cusp accurately. The MC-srDFT method avoids this by replacing an explicit description of the Coulomb cusp with an effective density functional. Hence, the basis set requirements can be expected to be similar to regular DFT, which is known for its fast convergence with respect to basis set expansion.

The first studies shown here are promising, but it should be emphasized that further developments are required. An obvious extension is to employ srDFT functionals that include the kinetic energy density, which has been developed for range-separated methods. 79 Also, a good measure of local spin density, for example in terms of the on-top pair density [START_REF] Perdew | Escaping the symmetry dilemma through a pair-density of spin-density functional theory[END_REF] , needs to be implemented to satisfactorily describe for example open-shell singlets which exhibit zero mean spin density. Extension of the open-shell algorithm to molecular properties via response theory will also be an important extension of the method. Developments in these directions are ongoing.

In the future it would also be interesting to extend our state-averaged version of MC-srDFT [START_REF] Hedegård | Density matrix renormalization group with efficient dynamical electron correlation through range separation[END_REF] , in particular for unbiased location of conical intersections.

  below after first defining the quantities involved in the long-range part of the energy in Eq. (3). The long-range Hamiltonian Ĥlr takes a form similar to the regular electronic (non-or scalar-relativistic) Hamiltonian in second quantization Ĥlr = pq h pq Êpq + 1 2 pqrs g lr pq,rs êpq,rs + V nn ,

  fulfill P b n = b n . The explicit expressions for the new short-range DFT contributions to the individual gradient and σ n element types are derived below. The MCSCF long-range contributions are equivalent to the terms in the regular MCSCF method which together with the second-order optimization algorithm are documented in the original works describing the Dalton implementation.[START_REF] Aa | A direct approach to second order MCSCF calculations using a Norm-Extended Optimization scheme[END_REF][START_REF] Aa | MCSCF optimization using the direct, restricted step second-order norm-extended optimization method[END_REF][START_REF] Aa | A direct, restricted-step, second order MCSCF program for large scale ab initio calculations[END_REF][START_REF] Aa | Efficient optimization of large scale MCSCF wave functions with a restricted-step algorithm[END_REF][START_REF] Aa | Electron correlation in molecules using direct second order MCSCF[END_REF]

2

 2 and [Fe(H 2 O) 6 ] 3+ , exploiting point group symmetry: the D 2h symmetry sub-group of D ∞h for O 2 , and the C i symmetry group for [Fe(H 2 O) 6 ] 3+ .

Figure 2 ,

 2 right). Thus, a value of µ = 0.4 also allows the long-range CAS wave function to include static correlation for the present multiconfigurational open-shell system.

2 O 2

 22 FIG.3. Left: Weight |C 0 | 2 of the dominant configuration in the long-range CAS(12,8) wave function for O 2 in the 3 Σ - g and 1 ∆ g states at the equilibrium bond distance. The srLDA and srPBE functionals give rise to practically identical long-range wave functions. Right: Singlettriplet splitting ∆E( 3 Σ - g → 1 ∆ g ).

  Appendix A: Gradient contributions from short-range GGA functionals For an srGGA functional, the short-range exchange-correlation energy density e sr-xc also depends on the three gradient variables ξ CC (r) = ∇ρ C (r)•∇ρ C (r), ξ SS (r) = ∇ρ S (r)•∇ρ S (r), and ξ CS (r) = ∇ρ C (r) • ∇ρ S (r) in addition to the electron charge density ρ C (r) and the electron spin density ρ S (r):E sr-xc [ρ C , ρ S , ξ CC , ξ SS , ξ CS ] = e sr-xc (ρ C (r), ρ S (r), ξ CC (r), ξ SS (r), ξ CS (r)) dr. (A1)where we for brevity define the integrand as e sr-xc . The derivatives required for the electronic gradient then become

  ( 3 Σ -g state) at 1.2075 Å (R eq. ) Dependence on the range-separation parameter µ of the HF-srDFT and CAS-srDFT energies for O 2 in its 3 Σ - g ground state for the equilibrium bond distance R eq (left) and a stretched bond distance R stretch (right). All calculations are with the cc-pVTZ basis set.

		148.5 148.0	O 2 HF-srLDA CAS(8,6)-sLDA CAS(12,8)-srLDA CAS(12,16)-srLDA		HF-srPBE CAS(8,6)-srPBE CAS(12,8)-srPBE CAS(12,16)-srPBE	148.5 148.0		O 2 ( 3 Σ -g state) at 2.0 Å HF-srLDA CAS(8,6)-sLDA CAS(12,8)-srLDA CAS(12,16)-srLDA HF-srPBE CAS(8,6)-srPBE CAS(12,8)-srPBE CAS(12,16)-srPBE
	Energy (E h )	149.5 149.0			Energy (E h )	149.5 149.0
		150.0				150.0
		150.5	1.0 µ-value	1.5	2.0	150.5	0.5	1.0 µ-value	1.5	2.0
	FIG. 2.

  . Our best theoretical CAS(12,16)-srPBE result yields a splitting of 123 kJ/mol, which is in between the CCSD(T) and CCSDT results, but (as expected) too large compared to the experimental value of 94.6 kJ/mol[START_REF] Huber | Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules[END_REF] . The fact that even CCSDT and some of the MRCI results[START_REF] Klotzchristel | Calculation of Spin-Forbidden Radiative Transitions Using Correlated Wave Funnctions: Lifetimes of b 1 Σ + , a 1 ∆ States in O 2 , S 2 and SO[END_REF] overestimate the experiment value by 20 and 16 kJ/mol, respectively, emphasizes that the calculation of this singlet-triplet splitting energy is indeed delicate. Only the highly accurate CCSDTQ and QMC methods are able to achieve a result that is within 5 kJ/mol of the experimental value. A peculiar observation is that the small

TABLE I .

 I Vertical singlet-triplet splitting, ∆E( 3 Σ - g → 1 ∆ g ) (in kJ/mol), for O 2 at the equilibrium bond distance, calculated with CAS-srDFT employing different srDFT functionals and values of µ (in bohr -1 ). To improve spin-state energetics, we are planning to make a number of improvements on the srDFT functionals, e.g., employing local spin densities. One option could be to employ the on-top pair density[START_REF] Li Manni | Multiconfiguration Pair-Density Functional Theory[END_REF][START_REF] Perdew | Escaping the symmetry dilemma through a pair-density of spin-density functional theory[END_REF] for this purpose.B. The6 A g → 4 T 2g spin-state splitting in [Fe(H 2 O) 6 ] 3+

	Method	µ = 0.0 µ = 0.1 µ = 0.3 µ = 0.4 µ = 0.5 µ = 1.0 µ → ∞
	CAS(12,8)-srLDA	96.96	95.11	95.18	95.82	96.86	102.24	92.35
	CAS(12,16)-srLDA 96.96	96.97	97.66	100.13	103.53	107.78 110.44
	CAS(12,8)-srPBE 127.62	123.79	121.24	120.36	119.88	118.92	92.75
	CAS(12,16)-srPBE 127.62	126.91	123.83	122.99	122.80	125.74 110.44
	Exp 70				94.6
	active space CAS(12,8) is very close to the experimental result, and CAS(12,8)-srLDA is
	almost identical to experiment. Both results are due to fortuitous error cancellation, as
	extending the active space leads to deterioration. Yet, CAS(12,16) and CAS(12,16)-srLDA
	are in fact still closer to experiment than both CCSDT and MRCI, showing that assessment
	of the functionals can require some consideration of the active space employed. Finally, we
	note that approaches based on perturbation theory (e.g. CASPT2 or NEVPT2) are likely
	also to provide a reasonable spin-state splitting. Indeed, we performed a NEVPT2(12,8)
	calculation which gives a singlet-triplet splitting of 92.35 kJ/mol. However, an in-depth
	discussion of CASPT2/NEVPT2 results is out of the scope for our current paper, where we
	focus on MC-srDFT.

TABLE II .

 II Singlet-triplet splitting, ∆E( 3 Σ - g → 1 ∆ g ) (in kJ/mol), for O 2 at the equilibrium bond distance. The CAS-srLDA and CAS-srPBE results are for µ = 0.4 bohr -1 . Note that Refs. 52,53 report the adiabatic energy difference, but with a difference of 0.008 Å between R e for the two states, this is expected to be of minor importance.

	Method	∆E( 3 Σ -g → 1 ∆ g )
	Our work (all with cc-pVTZ):
	CAS(12,8)	92.35
	CAS(12,8)-srLDA	95.82
	CAS(12,8)-srPBE	120.36
	CAS(12,16)	110.44
	CAS(12,16)-srLDA	100.13
	CAS(12,16)-srPBE	122.99
	Literature:
	CCSD(T)/cc-pVTZ 52	124.66
	CCSDT/cc-pVTZ 52	114.34
	CCSDTQ/cc-pVTZ 52	99.48
	MRCI 51	111.61
	MRCI/cc-pVTZ 52	98.40
	QMC 53	96.58
	Exp. 70	94.6

TABLE III .

 III Vertical quartet-sextet splitting (in eV) of [Fe(H 2 O) 6 ] 3+ , calculated with CAS-srPBE employing different values of µ (in bohr -1 ). method seem to be reasonably accurate[START_REF] Harris | Structure and Relative Spin-State Energetics of [Fe(H2O)6]3+: A Comparison of UHF, Møller-Plesset, Nonlocal DFT, and Semiempircal INDO/S Calculations[END_REF]78 . A recent study by Radoń et al.[START_REF] Radoń | Spin-State Energetics of Fe(III) and Ru(III) Aqua Complexes: Accurate ab Initio Calculations and Evidence for Huge Solvation Effects[END_REF] pointed out some of the possible failures of accurate ab initio methods. First, it was argued that the second solvation shell was important, mainly for the molecular structure of the [Fe(H 2 O) 6 ] 3+ core. If this effect is taken into account, even calculations on [Fe(H 2 O) 6 ] 3+ are in significant better agreement with experiment. Further, Radoń et al.[START_REF] Radoń | Spin-State Energetics of Fe(III) and Ru(III) Aqua Complexes: Accurate ab Initio Calculations and Evidence for Huge Solvation Effects[END_REF] argued that in addition to the d orbitals, the employed CAS space should at least contain two σ Fe-O bonding orbitals.[START_REF] Yang | Multireference Ab Initio Study of Ligand Field d-d Transitions in Octahedral Transition-Metal Oxide Clusters[END_REF]

	Method	µ = 0.01 µ = 0.1 µ = 0.33 µ = 0.4 µ = 0.5 µ = 0.75 µ = 1.0 µ = ∞
	CAS(9,12)-srPBE	1.15	1.15	1.35	1.44	1.61	1.74	2.10	2.24

  For values in the interval µ = 0.33 to µ = 0.50 bohr -1 expected a priori to be optimal, we do obtain good agreement with the experimental value of 1.56 eV. The "end-points" in Table

	III corresponds to a pure CAS(9,12) (µ = ∞) and pure PBE (µ = 0.01); these respectively
	overestimate and underestimate the experimental value significantly. It is gratifying the
	range of µ-values expected to provide the best results, indeed all are the best ones obtained
	(an in fact, all closer to experiment than CASPT2 with srPBE). Yet, before more general

[START_REF] Chan | The Density Matrix Renormalization Group in Quantum Chemistry[END_REF]

A g → 4 T 1g ). Results for selected values of µ are shown in Table

III

. Compared to the spin-state splitting in O 2 (see previous section), the ∆E( 6 A g → 4 T 1g ) spin-state splitting in the range µ = 0.4 -1.0 bohr -1 varies somewhat more with 1.15-2.10 eV for CAS(9,12)-srPBE. conclusions concerning the performance can be made, a larger span of systems must be explored. In Table

IV

, we compare the CAS

(9,12)

-srPBE result at µ = 0.4 bohr -1 with the CAS(9,12)-srLDA one, as well as results obtained with a number of other theoretical methods. As mentioned above, the literature results show a rather large scatter, and all

TABLE IV .

 IV Vertical quartet-sextet splitting (in eV) for [Fe(H 2 O) 6 ] 3+ . The CASPT2, CAS(9,12)-srLDA, and CAS(9,12)-srPBE calculations are for the [Fe(H 2 O) 6 ] 3+ core with a molecular structure that included the second solvation shell in the optimization, see Ref.66. The CAS(9,12)-srLDA and CAS(9,12)-srPBE calculations are for a range-separation parameter of µ = 0.4 bohr -1 .

	Method	∆E( 6 A g → 4 T 1g )
	[Fe(H 2 O) 6 ] 3+ optimized in vacuum
	Literature:
	CCSD(T) 54	2.47
	SORCI 74	2.07
	CASPT2 54	2.64
	[Fe(H 2 O) 6 ] 3+ optimized with a second solvation sphere
	Literature:
	CASPT2 66	1.96
	Our work:
	CAS(9,12)-srLDA	1.14
	CAS(9,12)-srPBE	1.44
	Exp. 72	1.56
	the values for the quartet-sextet splitting are too large. Our result with CAS(9,12)-srLDA
	shows the opposite pattern, it underestimates the experimental value significantly. The
	corresponding values for µ = 0.33 and µ = 0.5 are 1.04 eV and 1.29 eV and are thus still
	significantly underestimated. However, the gradient-corrected srDFT functional significantly
	improves the spin-state splitting energy for the value µ = 0.4. This shows that a GGA srDFT
	functional is important, as it also was for O 2 . We found in previous section that µ = 0.4
	was a good compromise between flexibility of the long-range wave function to include static
	correlation and ability of the srDFT functional to recover dynamical correlation. This value
	also provides a very accurate estimate for the quartet-sextet splitting ∆E( 6 A

g → 4 T 1g ) with CAS

(12,9)

-srPBE. The result is 1.44 eV, only 0.12 eV below the experimental value. In comparison, a calculation with CAS(9,12)PT2 on the same [Fe(H 2 O) 6 ] 3+ molecular structure gives a spin-state splitting of 1.96 eV, which is 0.44 eV too high. We note that µ = 0.33 and µ = 0.5 also give splittings better than the CAS(9,12)PT2 values.

  and the derivatives required for the electronic Hessian become∂ 2 e sr-xc ∂λ i ∂λ j = X=C,S ∂e sr-xc ∂ρ X ∂ 2 ρ X (r) ∂λ i ∂λ j + ∂e sr-xc ∂ξ XX ∂ 2 ξ XX (r) ∂λ i ∂λ j + ∂e sr-xc ∂ξ CS ∂ 2 ξ CS (r) ∂λ i ∂λ j

	+	X,Y =C,S			∂ 2 e sr-xc ∂ρ Y ∂ρ X	∂ρ X (r) ∂λ i	∂ρ Y (r) ∂λ j	+	∂ 2 e sr-xc ∂ρ Y ∂ξ XX	∂ξ XX (r) ∂λ i	∂ρ Y (r) ∂λ j
	+	∂ 2 e sr-xc ∂ξ Y Y ∂ρ X	∂ρ X (r) ∂λ i	∂ξ Y Y (r) ∂λ j	+	∂ 2 e sr-xc ∂ξ Y Y ∂ξ XX	∂ξ XX (r) ∂λ i	∂ξ Y Y (r) ∂λ j
	+	X=C,S	∂ 2 e sr-xc ∂ξ CS ∂ρ X	∂ρ X (r) ∂λ i	∂ξ CS (r) ∂λ j	+	∂ 2 e sr-xc ∂ξ CS ∂ξ XX	∂ξ XX (r) ∂λ i	∂ξ CS (r) ∂λ j
	+	∂ 2 e sr-xc ∂ξ 2 CS	∂ξ CS (r) ∂λ i	∂ξ CS (r) ∂λ j	.	(A3)
	The first-order derivatives of the three gradient terms are
	∂ξ CC (r) ∂λ i	= 2				∂D C pq ∂λ i	,	(A4)
	∂ξ SS (r) ∂λ i	= 2				∂D S pq ∂λ i	,	(A5)
	∂ξ CS (r) ∂λ i	=					∂D S pq ∂λ i	+ ∇ρ S (r)	∂D C pq ∂λ i	,	(A6)
	and the second-order derivatives are
	∂ 2 ξ CC (r) ∂λ i ∂λ j	= 2	pq,rs	(∇Ω pq (r) • ∇Ω rs (r))	∂D C rs ∂λ j	∂D C pq ∂λ i
			+ 2		pq	(∇Ω pq (r) • ∇ρ C (r))	∂ 2 D C pq ∂λ i ∂λ j	,	(A7)
	∂ 2 ξ SS (r) ∂λ i ∂λ j	= 2	pq,rs	(∇Ω pq (r) • ∇Ω rs (r))	∂D S rs ∂λ j	∂D S pq ∂λ i
			+ 2		pq	(∇Ω pq (r) • ∇ρ S (r))	∂ 2 D S pq ∂λ i ∂λ j	,	(A8)
	∂e sr-xc ∂λ i ∂ 2 ξ CS (r) ∂λ i ∂λ j	= + = +	∂e sr-xc ∂ρ C ∂e sr-xc ∂ξ CC pq,rs	∂ρ C (r) ∂λ i ∂ξ CC (r) + ∂λ i + ∂e sr-xc ∂ρ S ∂e sr-xc ∂ρ S (r) ∂λ i ∂ξ SS ∂ξ SS (r) ∂λ i ∂D C rs ∂λ j ∂D S ∂e sr-xc pq ∂λ i ∂D S rs + ∂λ j ∂ξ CS (r) ∂D C pq ∂λ i + ∂ξ CS ∂λ i
										(A2)

pq (∇Ω pq (r) • ∇ρ C (r)) pq (∇Ω pq (r) • ∇ρ S (r)) pq ∇Ω pq (r) • ∇ρ C (r) (∇Ω pq (r) • ∇Ω rs (r)) pq ∇Ω pq (r) • ∇ρ C (r) ∂ 2 D S pq ∂λ i ∂λ j + ∇ρ S (r) ∂ 2 D C pq ∂λ i ∂λ j . (

A9

)

ACKNOWLEDGMENTS

EDH thanks the Carlsberg foundation and the European Commission for post-doc stipends. JT and HJJ acknowledge support from the French government via the French-Danish collaboration project "WADEMECOM.dk".


Nadine Malhouroux-Gaffet

Patrice Pajusco
email: patrice.pajusco@telecom-bretagne.eu

Roxana Burghelea

Christian Leray
email: christian.leray@time-reversal-communications.fr

Nadine Malhouroux Gaffet

Capacity Gain of MIMO Systems with Micro-structured Antenna Arrays

Keywords: antenna, propagation, measurement, MIMO. I

published or not. The documents may come

INTRODUCTION

MIMO communication systems (multiple input multiple output) have created a technological breakthrough allowing a significant capacity increase of the transmission channels. In this context, the antennas that provide the interface between the circuit and the propagation medium are essential but complex due to physical constraints (size, electronic circuitry integration), technologic limitations (radiation, efficiency and bandwidth) or economic restraints (cost). New antenna design, using the hybridization band gap concept offers interesting prospects for miniaturization or large bandwidth sensors.

Such antennas are very compact and allow integrating a large number of highly uncorrelated sensors in a volume corresponding to WLAN equipment or to FEMTO equipment. The objective of this work is to evaluate the capacity gain of this type of antenna for MIMO systems.

II. DESCRIPTION OF THE ANTENNAS ARRAYS.

A. Hybridization Band Gap Based Antenna array

The antenna array [START_REF] Chu | Physical Limitations of Omnidirectional Antennas[END_REF] developed by the company Time Reversal Communications is shown in Fig. 1 and Fig. 2 (four port array). The antenna arrays consist of a set of close and small monopoles, which are embedded in a dense arrangement of parasitic radiators. These last have two functions: efficiently radiate the field produced by the small monopoles and ensure a good isolation between them. In such a case, the field radiated by each port of the antenna is de-correlated. In other words, the radiation patterns of each port are orthogonal. To obtain this result, this antenna design uses the property of the resonant behavior of the parasitic radiators that leads to efficiently convert evanescent waves (sub-wavelength field fluctuations) into propagative ones. This result is in accordance with the Chu's theorem which implies that an efficient small antenna should be resonant [START_REF] Lerosey | Hybridization band gap based smart antennas: Deep subwavelength yet directional and strongly decoupled MIMO antennas[END_REF].

Two micro-structured antenna arrays have been investigated in this study: At first time a four element array and in a second time an eight element array. The four port antenna structure are spaced of about λ/10 on a FR4 substrate [START_REF] Balanis | Antenna Theory: Analysis and Design[END_REF]. Fig. 3 shows the variation of the reflection coefficient S 11 as a function of frequency. Each port of the antenna has been measured and shows that it is very well adapted in the 2.4-2.5 GHz band.

To characterize and evaluate the 3D radiation pattern of these antennas, a measurement campaign has been carried out in near field condition which is presented on Fig. 4b. The antenna pattern is measured at 2.45 GHz for four ports.

B. Reference Antennas array

The reference antenna array (Fig. 6) consists in two dual polarization square patches printed on a FR4 substrate. Each patch is fed by two microstrip lines supplying the two orthogonal polarizations (V/H). A line of quarter-wave length provides impedance matching. This structure is duplicated on the other sides. This solution enables a quite compact eight element array.

The antennas are adapted on the 2.4-2.5 GHz frequency band as shown by the S 11 coefficient measured in Fig. 7. Radiation patterns of the reference antenna array have been measured in the xOz plane of the Satimo SG32 (Fig4-a). Taking into account the symmetries to the ground plane, Fig. 7 shows the 3D radiation diagrams obtained from the measurements at the central frequency of 2.45 GHz for the four ports on the front face of the antenna array. Pattern and gains are almost similar; the only difference comes from the direction of polarization.

C. Comparison of the two types of antenna array

Table I summarizes the performances of these two antenna prototypes measured at 2.45 GHz.The performances of the two antenna arrays are similar in terms of bandwidth, isolation and correlation of ports. Only the gain of the reference patch array is slightly smaller than the band gap antenna. This can be explained by the difference of loss in the substrate and by the lines supplying the ports which are printed for the reference array and composed of little cables for the micro-structured antennas. Figure 9 shows the equipment carried out for this study. A vector network analyzer (VNA -ZVB20) is used allowing the measurement of the propagation channel over a wideband. The ports 1 and 2 are connected to the transmitter and to the receiver respectively. The impulse responses are obtained measuring the S 21 parameter channel (including antenna arrays). For this study, the channel measurements are performed in an indoor environment. Two SP8T RF switches (Radio Frequency) are used both for the transmitter and the receiver allowing the selection of Tx and /or Rx antenna ports. Thus, all possible combinations between the transmitting and receiving antennas can be measured. All elements of the acquisition chain are controlled by a computer and the measurement process is set through a graphic user interface developed in Matlab ™ . This interface provides calibration process, the switching of the antennas, the channel measurement, and the post processing of data.

This experimental set-up is designed for the measurement of 8×8 MIMO channels. The channel can be measured with different antenna arrays and thus their influence on the channel capacity can be estimated. A validation of the equipment was performed in an anechoic chamber.

IV. DESCRIPTION OF THE INDOOR ENVIRONMENT

MIMO channel measurements have been achieved in a typical indoor environment at Orange Labs premise (Belfort) that is statistically representative of an apartment. The layout and the measurement locations are shown in Fig. 10. The transmitter was located in 5 different locations, while the receiver was placed at a single location (Table II). The reception and transmission antenna arrays have been positioned at the same height (1 meter). Four different orientations of the arrays were considered: 0 °, 45 °, 90 ° and 135 ° for the receiver and the transmitter. Finally 10,000 impulse responses in LOS (Line Of Sight) and NLOS (No Line Of Sight) configurations have been recorded. T 5 corresponds to a delay max of 400 ns and a resolution of 0.1 ns.

V. ANALYSIS OF THE RESULTS

The MIMO analysis has been performed on the 4×4 and 8×4 configurations. Sub-antenna arrays have been created in order to maximize the number of channels:

• Micro-structured (TRCMA ) antennas at the transmitter and at the receiver. • Reference patch array at the transmitter and the receiver both located on the front side (P1, P2, P3, P4 ports). • Reference patch array at the transmitter and the receiver. Two patches on the front side and two patches on the back side (P1, P2, P5, P6 ports). • Reference patch array at the transmitter (P1, P2, P3, P4 ports) and micro-structured array at the receiver (TRCMA) • Reference patch array at the transmitter (P1, P2, P5, P6 ports) and micro-structured array at the receiver (TRCMA) • Reference patch array at the transmitter (P1, P2, P3, P4, P5, P6, P7, P8 ports) and micro-structured array at the receiver (MA8) • Reference patch array at the transmitter and at the receiver (P1, P2, P3, P4, P5, P6, P7, P8 ports).

A. MIMO 4×4

Water Filling capacities were calculated for each measurement configuration to quantify the contribution of micro-structured antenna for MIMO communication systems. The capacity is plotted at the center frequency of 2.4 GHz in function of the SNR (signal to noise ratio).

The comparison of the average water-filling capacity of the different configurations 4×4 is presented in the Fig. 11. An average of the different orientations of the antenna array is realized for each configuration and the measurements were normalized by the values of path loss. To improve the comparison, a reference plot has been added: the capacity of a Rayleigh channel 4×4. The curves are similar, including slopes that are substantially identical. The two types of antennas therefore are able to measure four degrees of freedom of the channel. However, these results highlight the increased channel capacity with micro-structured antennas that can be estimated of 3bits/s/Hz compared to the capacity obtained with the reference antennas. This gain corresponds to a value of 3 dB for the link budget of such configuration. It substantially corresponds to the efficiency ratio between the microstructured antenna array and the antenna reference array.

The variation of the channel capacity of all the measurements (minimum, maximum and average) is given in Fig. 12 for each 4×4 configuration. The dispersion of the curves is higher for the reference antenna array than for the micro-structured antenna.

B. MIMO 8×8

The capacity of 8×8 MIMO matrix measurement are shown in Fig. 13. The dispersion of the curves is consistent with the 4×4 systems results. An analysis of the results has been realized for different locations of the transmitter. Two typical results are presented in Fig. 14 and Fig. 15. In the first case, the transmitter is located in a room with two walls and furniture through the radio link. In this case, the gain of the micro-structured antennas is obvious. In a different configuration, when the configuration is Obstructed Line of Sight (OLOS), the gain of the micro-structured antenna is not observable (Fig. 15). These results are consistent with the number of multipath of the two radio link configurations. Indeed, some previous works [START_REF] Malhouroux-Gaffet | Effect of propagation phenomena on MIMO capacity of wireless systems between 3 and 10 GHz[END_REF] have shown that the number of "pertinent" multipath is higher for the T2 transmitter than for the T3 transmitter. So the microstructured antennas are well adapted for multipath channel.

VI. CONCLUSION

The purpose of this study was to evaluate experimentally the micro-structured antennas developed under the project TRIMARAN by Time Reversal Communications Company. A MIMO channel sounder has been developed and used to record MIMO measurements in indoor environment. Initial results confirm that the micro-structured antennas provide a diversity similar to a reference antenna array but in a smaller volume. The MIMO capacity is less dispersive with micro-structured antenna array and better in configuration with many multipath.

Fig. 1 .

 1 Fig. 1. Microstructured Antenna array.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Four ports and eight ports micro-structured antenna .

Fig. 4 .

 4 Fig. 4. Antenna array measurement and associated coordinated system.

Fig. 5 .

 5 Fig. 5. Radiation pattern of the micro-structured antenna array at the central resonnance frequency f=2.45 GHz.

Fig. 6 .

 6 Fig. 6. Reference Antenna.

Fig. 7 .

 7 Fig. 7. Reflection Coefficient of the reference antenna array.

Fig. 8 .

 8 Fig. 8. Radiation pattern of the reference antenna array at the central resonnance frequency f=2.45 GHz.

Fig. 9 .

 9 Fig. 9. Set-up of the MIMO channel sounder.

Fig. 10 .

 10 Fig. 10. Indoor environment with transmitter and receiver locations

Fig. 11 .

 11 Fig. 11. MIMO Waterfilling Capacity for 4×4 Systems-average for all transmitter positions and orientation and all receiver orientations

  Fig. 13. MIMO Waterfilling Capacity for 8×8 Systems-average, max,min for all transmitter positions and orientation and all receiver orientations

Fig

  Fig. 14. MIMO Waterfilling Capacity for 8×8 Systems-average, max,min for T2 transmitter positions and orientation and all receiver orientations .

TABLE I

 I

			.	ANTENNA MEASURED PERFORMANCES
		Measurement types	Antenna characteristics Reference patch Hybridization band antenna gap antenna
		Gain max			1.2-2.2 dB	1.3-2.7 dB
		Efficiency			29-33%	45-60%
		Correlation			<0.15	<0.15
		S11 (100 MHz)		-6 to -20 dB	-8 to -25 dB
		Isolation min			-13 dB	-12 dB
						Analysis and storage
		USB cable				USB cable
					GPIB	Measurement data
	DeviceI/O		RF cable	VNA		DeviceI/O
				Port-1	Port-2
	USB						USB
	cable	Switch	Tx			Rx	Switch	cable
				Radio channel
		Transmitter			Receiver

III. MIMO CHANNEL SOUNDER

A MIMO channel sounder has been developed for this experimentation in order to measure the MIMO matrix of this 4×4 or 4×8 antennas system.

TABLE II .

 II TRANSMITTER AND RECEIVER LOCATIONS COORDINATES

			12,45m
					135°	90°
			R 0		45°
					0°
	7,27m	T 1	T 3
		T 2			T 4	Rx
					Tx
		Locations X(m) Y(m) Z(m)
		R0	6.46	5.88	1
		T1	1.40	3.86	1
		T2	5.28	1.08	1
		T3	5.35	2.74	1
		T4	10.53	1.48	1
		T5	10.70	3.98	1
	The transmission channel has been characterized on 401
	frequencies between 2.2 and 2.7 GHz. This configuration

ACKNOWLEDGMENT

The authors thank the National Research Agency and the Pole Images and Networks for funding of TRIMARAN [5] project in which this work was conducted.