Jacques K Desmarais

Alessandro Erba

Jean-Pierre Flament

Bernard Kirtman

Perturbation Theory Treatment of Spin-Orbit Coupling. II: Coupled Perturbed Kohn Sham Method

A non-canonical coupled-perturbed Kohn-Sham density functional theory (KS-DFT)/Hartree-Fock (HF) treatment of spin-orbit coupling (SOC) is provided. We take the scalar-relativistic KS-DFT/HF solution, obtained with a relativistic effective core potential, as the zeroth-order approximation. Explicit expressions are given for the total energy through 4th-order, which satisfy the 2n + 1 rule. Second order expressions are provided for orbital energies and density variables of spin-current DFT. Test calculations are carried out on the halogen homonuclear diatomic and hydride molecules, including 6p and 7p elements, as well as open-shell negative ions. The computed properties through second-or third-order match well with those from reference two-component selfconsistent field calculations for total and orbital energies as well as spin-current densities. In only one case (At - 2 ) did a significant deviation occur for the remaining density variables. Our coupled perturbation theory approach provides an efficient way of adding the effect of SOC to a scalar relativistic single-reference KS-DFT/HF treatment, in particular because it does not require diagonalization in the two-component spinor basis, leading to saving factors on the number of required floating point operations that may exceed one order of magnitude.

I. INTRODUCTION

Relativistic effects are typically accounted for in modern electronic structure programs through self-consistent field (SCF) treatments, either in a two-or fourcomponent spinor basis (2c-SCF or 4c-SCF). 1-18 Such variational treatments have been made efficient, in particular for closed-shell systems, within a Kramersrestricted framework. In that context approaches have been developed to exploit time-reversal and double-group symmetries, which greatly reduce the computational task of diagonalizing the associated secular equations. [START_REF] Saue | Quaternion symmetry in relativistic molecular calculations: The Dirac-Hartree-Fock method[END_REF][START_REF] Repisky | Relativistic spec-troscopy DFT program package[END_REF][START_REF] Peng | On the construction of Kramers paired double group symmetry functions[END_REF][START_REF] Yanai | New implementation of molecular double pointgroup symmetry in four-component relativistic Gaussiantype spinors[END_REF][START_REF] Armbruster | Quaternionic formulation of the twocomponent Kohn-Sham equations and efficient exploitation of point group symmetry[END_REF] In practice, however, the symmetries are usually not exploited for diagonalization, due in part to a lack of efficient algebraic routines.

An alternative to the 2c-SCF or 4c-SCF treatment of relativistic effects is to include them from a perturbation theory approach. The all-electron four-component Direct perturbation theory (DPT) of Rutkowski and Kutzelnigg follows such a prescription. [START_REF] Rutkowski | Relativistic perturbation theory. I. A new perturbation approach to the Dirac equation[END_REF][START_REF] Rutkowski | Relativistic perturbation theory: II. Oneelectron variational perturbation calculations[END_REF][START_REF] Rutkowski | Relativistic perturbation theory. III. A new perturbation approach to the two-electron Dirac-Coulomb equation[END_REF][START_REF] Kutzelnigg | Stationary direct perturbation theory of relativistic corrections[END_REF][START_REF] Van Wüllen | A relativistic Kohn-Sham density functional procedure by means of direct perturbation theory[END_REF][START_REF] Stopkowicz | Direct perturbation theory in terms of energy derivatives: Fourth-order relativistic corrections at the Hartree-Fock level[END_REF][START_REF] Stopkowicz | A one-electron variant of direct perturbation theory for the treatment of scalar-relativistic effects[END_REF][START_REF] Cheng | Spin-free Dirac-Coulomb calculations augmented with a perturbative treatment of spin-orbit effects at the Hartree-Fock level[END_REF][START_REF] Cheng | Perturbative treatment of spin-orbit-coupling within spin-free exact two-component theory using equation-of-motion coupled-cluster methods[END_REF] One advantage of such a perturbation treatment is the possibility of accounting for multi-reference character based on quasidegenerate perturbation theory (QDPT), which could provide a framework for developing a cost effective computational procedure. [START_REF] Kutzelnigg | Relativistic MCSCF by means of quasidegenerate direct perturbation theory[END_REF][START_REF] Liu | Relativistic MC-SCF by means of quasidegenerate direct perturbation theory. II. Preliminary applications[END_REF][START_REF] Rutkowski | Effective Hamiltonian for near-degenerate states in direct relativistic perturbation theory. I. Formalism[END_REF][START_REF] Rutkowski | Effective Hamiltonian for near-degenerate states in relativistic direct perturbation theory. II. H 2+like systems[END_REF][START_REF] Kutzelnigg | Matrix formulation of direct perturbation theory of relativistic effects in a kinetically balanced basis[END_REF][START_REF] Kirtman | Simultaneous calculation of several interacting electronic states by generalized Van Vleck perturbation theory[END_REF][START_REF] Kirtman | Variational form of Van Vleck degenerate perturbation theory with particular application to electronic structure problems[END_REF][START_REF] Kirtman | Interchange Theorems in Degenerate Perturbation Theory[END_REF][START_REF] De Heer | Interchange Theorems, Accessible Correlation Expressions, and Physical Properties[END_REF] A particularly convenient representation of the Dirac equation from a computational perspective is provided through the relativistic effective core potential (RECP) approach, which allows for treating scalar-relativistic (SR) effects in a manner that is no more expensive than a one-component non-relativistic calculation (1c-SCF). [START_REF] Dolg | Relativistic pseudopotentials: their development and scope of applications[END_REF][START_REF] Ermler | AB initio effective core potentials including relativistic effects. A procedure for the inclusion of spinorbit coupling in molecular wavefunctions[END_REF] Such one-component SR RECP treatments, however, do not include spin-orbit coupling (SOC) effects. In Part I of this series we presented a double perturbation theory for adding SOC and dynamical electron correlation effects to the SR RECP solution. [START_REF] Desmarais | Perturbation Theory Treatment of Spin-Orbit Coupling. Part I: Double Perturbation Theory from a Single-Reference Initial Approximation[END_REF] This wavefunction treatment also provides a convenient starting point for the generalization to QDPT.

In this paper (Part II of the series) we develop a coupled perturbed Kohn-Sham (CPKS) treatment for including SOC effects within density functional theory (DFT). Our treatment follows along lines related to those developed for the response of molecular and periodic systems to external fields. [START_REF] Ferrero | The calculation of static polarizabilities of 1-3D periodic compounds. the implementation in the crystal code[END_REF][START_REF] Ferrero | Calculation of first and second static hyperpolarizabilities of one-to three-dimensional periodic compounds. Implementation in the CRYSTAL code[END_REF][START_REF] Ferrero | Coupled perturbed Hartree-Fock for periodic systems: the role of symmetry and related computational aspects[END_REF][START_REF] Ferrero | Coupled perturbed Kohn-Sham calculation of static polarizabilities of periodic compounds[END_REF][START_REF] Kirtman | Electric field polarization in conventional density functional theory: From quasilinear to two-dimensional and threedimensional extended systems[END_REF][START_REF] Maschio | Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis[END_REF][START_REF] Maschio | Ab initio analytical Raman intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method in an atomic orbital basis. II. Validation and comparison with experiments[END_REF][START_REF] Maschio | Ab initio analytical infrared intensities for periodic systems through a coupled perturbed Hartree-Fock/Kohn-Sham method[END_REF][START_REF] Maschio | Coupled perturbation theory approach to dual basis sets for molecules and solids. 1. General theory and application to molecules[END_REF][START_REF] Johnson | An implementation of ana-lytic second derivatives of the gradient-corrected density functional energy[END_REF][53] The perturbation theory treatment represents an efficient means of adding the effect of SOC to an SR single-reference KS-DFT/HF solution, in particular because it does not require diagonalization in the enlarged 2c-SCF spinor basis.

For SOC in a two-component spinor basis, the appropriate generalization is based on so-called spin-current DFT (SCDFT), which was developed originally for treating molecules in an external magnetic field by Vignale and Rasolt [START_REF] Vignale | Current-and spin-densityfunctional theory for inhomogeneous electronic systems in strong magnetic fields[END_REF][START_REF] Vignale | Density-functional theory in strong magnetic fields[END_REF] and subsequently generalized for SOC by Bencheikh. [START_REF] Bencheikh | Spin-orbit coupling in the spin-currentdensity-functional theory[END_REF] In SCDFT, the exchange-correlation (xc) functional depends, not only on the particle-number density, n, and z-component of the magnetization, m z , (as in the usual spin-DFT of von Barth and Hedin), [START_REF] Von Barth | A local exchange-correlation potential for the spin polarized case[END_REF] but also on the other Cartesian components of the magnetization, m x and m y , as well as the orbital-j and spincurrent J x , J y and J z densities. The dependence of the SCDFT functionals on such a large set of auxiliary density variables renders their parametrization a formidable task. Thus, SCDFT has seldom been used in actual applications to material systems. [START_REF] Pittalis | U (1)× SU (2) gauge invariance made simple for density functional approximations[END_REF][START_REF] Trushin | Spin-current density-functional theory for a correct treatment of spin-orbit interactions and its application to topological phase transitions[END_REF] However, some of the present authors have recently shown that j, J x , J y and J z , can be naturally included in xc functionals formulated within local-density and generalized-gradient approximations (LDA and GGA) of SCDFT by including a fraction of non-local Fock exchange in the single-particle Hamiltonian. [START_REF] Desmarais | Adiabatic connection in spin-current density functional theory[END_REF] This paper provides the corresponding perturbation theory formulation.

Our paper is organized as follows. In section II we outline the perturbation theory formulation including noncanonical solution of the perturbation equations, 2n + 1 rule energy expressions, and determination of orbitals energies as well as the density variables. Section III provides a step-by-step outline of the computational procedure whereas the specific details of our CPKS calculations are given in Section IV along with those for the preliminary 1c-SCF and reference 2c-SCF treatments. In section V we present results for the halogen homonuclear diatomic and hydride molecules. Section VI provides a comparison of computational requirements for the CPKS vs. 2c-SCF approaches. Conclusions and suggestions for future work are drawn in section VII. Additional mathematical details on the derivation of the order-by-order energy expressions are provided in the Appendix.

II. FORMALISM A. Statement of the Problem

Following Ref. 60, the adiabatic connection of the spincurrent density functional theory (SCDFT) allows us to write the exchange-correlation (xc) energy E xc in terms of fully-interacting E F I and non-interacting E N I contributions. In the local-density and generalized-gradient approximations (LDA and GGA) of the SCDFT, we have: [START_REF] Desmarais | Adiabatic connection in spin-current density functional theory[END_REF] E xc n, m, j, J x , J y , J z = E F I [n, m] +aE N I n, m, j, J x , J y , J z ,

where, n is the particle-number density, m is the Cartesian magnetization vector, j is the orbital-current density, J x , J y , J z are the three Cartesian spin-current densities and a is the dimensionless fraction of exact-exchange. The xc energy of the fully interacting system is, then, written as the sum of pure exchange, E x , and correlation, E c , contributions:

E F I [n, m] = (1 -a)E x [n, m] + E c [n, m] . (2)
As was the case in Part I, [START_REF] Desmarais | Perturbation Theory Treatment of Spin-Orbit Coupling. Part I: Double Perturbation Theory from a Single-Reference Initial Approximation[END_REF] the formalism presented here is based on the relativistic effective-core potential approximation (RECP). Both scalar-relativistic (SR) and spin-orbit coupling (SOC) effects are treated in a computationally convenient way by writing the Hamiltonian in a basis of two-component spinors. The two-component spinors |i , in turn, are expanded in the chosen basis as:

r|i = µ c α µ,i r|µ ⊗ |α + c β µ,i r|µ ⊗ |β , (3)
where the c α µ,i are the KS orbital coefficients and ⊗ represents the direct product with the spin functions |α and |β . Here, the r|µ = χ µ (r) are contracted atomcentered Gaussian functions = atomic orbitals (AOs).

In using the expansion of Eq. ( 3) in an AO basis with a size N B , we define bold matrices with two spin indices as being of dimension N B × N B , while bold matrices with only one spin index will have the dimension 2N B × N B . Then the orthonormality condition involving the spinblocks of the overlap matrix S σσ , for σ = α or β and those of the KS orbital coefficients c σ may be written:

c α c β † S αα 0 0 S ββ c α c β = I , ( 4
)
where I is the identity matrix. We will use an underlined notation to denote super-matrices that contain all of the spin-blocks, so that Eq. ( 4) can be written in a more compact way as:

c † S c = I . (5)
The corresponding SCDFT single-particle equation in matrix form reads:

H c = S c . (6)
or, equivalently:

H αα H αβ H βα H ββ c α c β = S αα 0 0 S ββ c α c β , (7)
where

H σσ = h σσ + C σσ -aK σσ + V σσ , (8)
in which h σσ contains the matrix elements that can be built from mono-electronic integrals:

h σσ = δ σ,σ [v + u AR ] + u σσ SO . (9)
Here, v consists of the electronic kinetic energy and electron-nuclear interaction terms, u AR and u σσ SO are, respectively, the averaged and spin-orbit relativistic effective potential (AREP and SOREP) matrices; and C σσ and K σσ are the usual Coulomb and exact-exchange terms. We also define V σσ as

V σσ = V σσ c + (1 -a)V σσ x , (10)
where V σσ c and V σσ x are the matrices that represent the DFT correlation and exchange potentials. The connection between the xc energy and potentials is discussed in detail in sections II H and the Appendix.

u σσ SO is the matrix (expressed in the basis of AOs) of the SOREP operator ûSO , which, in turn, is a sum of oneelectron operators ûSO,i that have the general form:

42 ûSO,i = ζi Lz,i Ŝz,i + 1 2 L+,i Ŝ-,i + 1 2 L-,i Ŝ+,i (11)
where ζi is a linear combination of one-electron operators with a radial part given by solid Gaussian functions multiplied by powers of the electron-core distance and an angular part given by projectors for each orbital angular-momentum component. Lz,i ( Ŝz,i ) is the z-component of the one-electron orbital (spin) angularmomentum operator, and L±,i ( Ŝ±,i ) are the rbital (spin) angular-momentum ladder operators.

From Part I we recall that (assuming real AOs) the diagonal spin-blocks of u SO in Eq. ( 9) are pure imaginary: 42

R [u σσ SO ] = 0 , (12)
whereas the off-diagonal spin-blocks are complex.

B. Expansion of the Hamiltonian Matrix in Orders of Perturbation Theory

In general Eq. ( 6) can be expanded as a power series in orders of the dimensionless perturbation-strength parameter λ:

H (0) + λH (1) + λ 2 H (2) + . . . c (0) + λc (1) + λ 2 c (2) + . . . = S c (0) + λc (1) + λ 2 c (2) + . . . (0) + λ (1) + λ 2 (2) + . . . . (13)
We take the SR Hamiltonian as the zeroth-order approximation:

H σσ (0) = δ σ,σ v + u AR + C σσ(0) -aK σσ(0) + V σσ(0) . (14a) In first-order H σσ (1) = C σσ (1) -aK σσ (1) + V σσ (1) + u σσ SO , (14b)
and for all orders greater than one (N > 1): N ) . N > 1 (14c) The Coulomb and exchange matrices in Eqs. (14a)-(14c) are defined as usual through the elements of the density matrix P σσ (N ) (for the definition of the latter see Eq. (34) below). Thus, for the Coulomb matrix:

H σσ (N ) = C σσ (N ) -aK σσ (N ) + V σσ (
C σσ (N ) µν = δ σ,σ τ ω R P αα(N ) ωτ + P ββ(N ) ωτ (µν|τ ω) , (15a
) which includes only the real part of the density matrix since the contribution of the imaginary part cancels due to the Hermiticity of P σσ(N ) and the fact that the AOs are real; while for the exchange matrix:

K σσ (N ) µν = τ ω P σσ (N ) τ ω (µτ |ων) , (15b)
where the bi-electronic integrals have been written in Mulliken notation.

In considering the contribution to Eqs. (14a)-(14c) from the xc potential, it is convenient to introduce the spin-density n σ(N ) :

n σ(N ) (r) = µν R P σσ(N ) µν χ µ (r) χ ν (r) , (16)
which may be expanded in powers of λ:

n σ (r) = n σ(0) (r) + λn σ(1) (r) + λ 2 n σ(2) (r) + . . . ( 17
)
For the specific case of a collinear formulation, Eqs. ( 16) and ( 17) allow us to define orders of perturbation theory for the matrix elements of the xc potential through the expansion:

V σσ µν = V σσ (0) µν + λV σσ (1) µν + λ 2 V σσ (2) µν + . . . = δ σσ µ| δE F I n α , n β δn σ |ν . (18)
The V σσ (N ) µν are, then, obtained by taking the derivative through order N of Eq. ( 18) and evaluating the result at λ = 0. This leads to expressions that have been developed elsewhere in the context of one-component SR theory, for the case in which the AOs are independent of the perturbation. [START_REF] Johnson | An implementation of ana-lytic second derivatives of the gradient-corrected density functional energy[END_REF][START_REF] Orlando | Calculation of the first static hyperpolarizability tensor of threedimensional periodic compounds with a local basis set: A comparison of LDA, PBE, PBE0, B3LYP, and HF results[END_REF] Here, we recall the results for the simplest case represented by the local spin-density approximation (LSDA). Given that

E F I n α , n β = F xc n α (r), n β (r) dr (19)
then

V σσ (0) µν = δ σ,σ µ| ∂F xc n α(0) , n β(0) ∂n σ(0) |ν , (20)
V σσ (1)

µν = δ σ,σ σ µ| ∂ 2 F xc n α(0) , n β(0) ∂n σ(0) ∂n σ (0) n σ (1) |ν , (21)
and:

V σσ (2) µν = δ σ,σ σ ,σ × µ| ∂ 3 F xc n α(0) , n β(0)
∂n σ(0) ∂n σ (0) ∂n σ (0) n σ (1) n σ (1) |ν

+ σ µ| ∂ 2 F xc n α(0) , n β(0) ∂n σ(0) ∂n σ (0) n σ (2) |ν . ( 22
)
Corresponding formulas can be obtained for noncollinear generalizations of E F I . [START_REF] Hobbs | Fully unconstrained noncollinear magnetism within the projector augmentedwave method[END_REF][START_REF] Kurz | Ab initio treatment of noncollinear magnets with the full-potential linearized augmented plane wave method[END_REF][START_REF] Kubler | Density functional theory of non-collinear magnetism[END_REF] This is achieved by, first, recasting E F I in terms of both the particle-number density n and the z-component of the magnetization m z :

n (N ) (r) = n α(N ) (r) + n β(N ) (r) (23a) m (N ) z (r) = n α(N ) (r) -n β(N ) (r) (23b)
Substitution of Eqs. (23a) and (23b) in Eq. ( 19) yields the energy of the fully interacting system E F I in terms of n and m z . Subsequently, a non-collinear generalization can be obtained by replacing m z with m = m 2 x + m 2 y + m 2 z = |m| (see Eqs. (68a) and (68b) below for definitions of m x and m y ). [START_REF] Hobbs | Fully unconstrained noncollinear magnetism within the projector augmentedwave method[END_REF][START_REF] Kurz | Ab initio treatment of noncollinear magnets with the full-potential linearized augmented plane wave method[END_REF][START_REF] Kubler | Density functional theory of non-collinear magnetism[END_REF] For functionals beyond the LSDA, we note that the development of numerically stable implementations that include noncollinear magnetizations is significantly more challenging. In particular, for GGA functionals, a careful treatment of otherwise unstable terms appearing in the xc energy and potential expressions is required, taking into account the need for proper reduction to the closed-shell and collinear limits to ensure rotational invariance in practical calculations. [START_REF] Desmarais | Spin-orbit coupling from a two-component selfconsistent approach. II. Non-collinear density functional theories[END_REF][START_REF] Komorovsky | Fourcomponent relativistic time-dependent density-functional theory using a stable noncollinear DFT ansatz applicable to both closed-and open-shell systems[END_REF][START_REF] Goings | Current development of noncollinear electronic structure theory[END_REF][START_REF] Petrone | An efficient implementation of two-component relativistic density functional theory with torque-free auxiliary variables[END_REF][START_REF] Scalmani | A new approach to noncollinear spin density functional theory beyond the local density approximation[END_REF][START_REF] Peralta | Noncollinear magnetism in density functional calculations[END_REF][START_REF] Bulik | Noncollinear density functional theory having proper invariance and local torque properties[END_REF] C. Order by Order Expressions for the Orthonormality Conditions

For further development we turn now to the perturbation expansion of the orthonormality condition given by Eq. (5). In this context it is useful to write the KS orbital coefficients c in terms of N B × N B spin-blocks with double spin-indices, so that Eq. ( 4) can be recast as:

c αα c αβ c βα c ββ † S αα 0 0 S ββ c αα c αβ c βα c ββ = I . (24)
In Eq. ( 24) and throughout this paper, the dagger † operates simultaneously on the spatial and spin parts of the orbital coefficients (i.e. the dagger swaps both the spatial and spin indices, so that

c βα i,ν † = c αβ ν,i *
). c βα , for example, represents the α 2c spinor coefficients that are generated by perturbation of the β MOs of the reference SR solution -see Eq. ( 28) below. Carrying out the multiplication of the spin-blocks in Eq. ( 24), we obtain

[c αα ] † S αα c αβ † S ββ c βα † S αα c ββ † S ββ c αα c αβ c βα c ββ = I . (25a)
whence

σ c σ σ † S σσ c σσ = I σ σ ⊗ δ σ ,σ , (25b)
Then, expanding Eq. (25b) as a power series in orders of λ gives

σ c σ σ † S σσ c σσ = σ c σ σ (0) † δ σ ,σ + λ c σ σ(1) † + λ 2 c σ σ(2) † + . . . S σσ × c σ σ (0) δ σ ,σ + λc σσ (1) + λ 2 c σσ (2) + . . . . (26)
which leads to

c σσ(0) † S σσ c σσ(0) = I σσ (27a) σ c σ σ(1) † S σσ c σ σ (0) δ σ ,σ + c σ σ (0) † δ σ ,σ S σσ c σσ (1) = 0 (27b) σ c σ σ(2) † S σσ c σ σ (0) δ σ ,σ + c σ σ (0) † δ σ ,σ S σσ c σσ (2) + c σ σ(1) † S σσ c σσ (1) = 0 (27c) . . . .
At this point it is convenient to introduce the N th order matrix U σ σ(N ) , defined through the relation:

c σ σ(N ) = c σ σ (0) U σ σ(N ) . ( 28
)
Substituting Eq. ( 28) in Eq. ( 27), we obtain the following constraints through fourth-order:

U σ σ(1) = -U σ σ(1) † , (29)
U σ σ (2) † + U σ σ (2) = - σ U σ σ(1) † U σσ (1) , (30)
U σ σ (3) † + U σ σ (3) = - σ U σ σ(2) † U σσ (1) + U σ σ(1) † U σσ (2) , (31)
and, finally

U σ σ (4) † + U σ σ (4) = - σ U σ σ(3) † U σσ (1) + U σ σ(1) † U σσ (3) + U σ σ(2) † U σσ (2) . ( 32
)
D. Order by Order Expressions for the Density Matrix

The single-particle density matrix may be written in terms of the KS orbital coefficients of Eq. ( 24), as - (33) where f σ is the diagonal matrix of (in the general case fractional) occupancies of the KS orbitals. For the ordinary gapped systems considered in the examples provided in section V, the elements of f σ are simply 1 for occupied orbitals and 0 for virtual orbitals. Expanding the orbital coefficients in Eq. (33) following Eq. ( 26), and using also Eq. ( 28), we obtain explicit expressions for the P σσ (N ) . In zeroth order:

P σσ = σ c σσ f σ c σ σ † = P σσ (0) + λP σσ (1) + λ 2 P σσ (2) + . . . ,
P σσ (0) = δ σ,σ c σσ(0) f σ c σσ(0) † , (34)
in first order:

P σσ (1) = c σσ(0) f σ U σσ (1) † c σ σ (0) † + c σσ(0) U σσ (1) f σ c σ σ (0) † , (35)
and in second order:

P σσ (2) = c σσ(0) f σ U σσ (2) † c σ σ (0) † + c σσ(0) U σσ (2) f σ c σ σ (0) † + σ c σσ(0) U σσ (1) f σ U σ σ (1) † c σ σ (0) † , (36)
The corresponding expressions in third and fourth order are provided in the Appendix.

E. The SCDFT Perturbation Equations

For convenience we introduce the Hermitian matrices G (N ) , corresponding to the perturbed Hamiltonian matrices in the KS orbital basis:

G (N ) = c (0) † H (N ) c (0) = G (N ) † . (37)
Returning to Eqs. ( 12) and (14b), we see that the contribution from u SO to the diagonal spin-blocks of G (1) are pure imaginary. As we will see later on, it turns out that the diagonal spin blocks of P (1) (and, hence H (1) ) are also pure imaginary, so that:

R G σσ(1) = 0 . ( 38
)
After multiplication of the left by c (0) † Eq. ( 13) may be recast in terms of the matrices G (N ) and U (N ) , which gives at order zero: 53

G (0) = (0) , (39)
in first order:

G (1) + (0) U (1) = U (1) (0) + (1) , (40)
in second order:

G (2) + (0) U (2) +G (1) U (1) = U (2) (0) +U (1) (1) + (2) , (41)
in third order: [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] , (42) and in fourth order:

G (3) + G (2) U (1) + G (1) U (2) + (0) U (3) = U (3) (0) + U (2) (1) + U (1) (2) + (
G (4) + G (3) U (1) + G (1) U (3) + G (2) U (2) + (0) U (4) = U (4) (0) + U (3) (1) + U (1) (3) + U (2) (2) + (4) . (43)
We note that (N ) and U (N ) , for N = 1 -4 are, to an extent, arbitrary since any definition is acceptable, provided that it is consistent with Eqs. ( 29)-( 32) and ( 40)- (43). As discussed below in section II F, here, we solve Eqs. ( 40)-( 43) making the simplest choice possible in terms of computation for (N ) and U (N ) . Our treatment is "non-canonical", in the sense that we obtain expressions for (N ) that are non-diagonal.

F. Non-Canonical Solution of the Perturbation Equations

Eqs. ( 39)-( 43) can be solved following the noncanonical treatment of Karna and Dupuis 53 as adapted for our purposes. In that treatment one takes advantage of the fact that the occupied-virtual (occ-virt) blocks of the perturbed Lagrange multiplier matrices vanish. That is to say, for m ∈ occ and p ∈ virt (N = 1 -4):

σσ (N ) mp = σ σ(N ) pm = 0 , (44)
Inserting Eq. ( 44) into Eq. ( 40) with N = 1, we obtain an explicit expressions for the occ-virt blocks of U (1) :

U σσ (1) mp = G σσ (1) mp σ (0) p - σ(0) m = -U σ σ(1) pm * . ( 45a
)
Returning to the first order orthogonality condition, we find that the anti-Hermiticity of the occ-virt blocks of U (1) in Eq. ( 45a) is consistent with Eq. ( 29). On the other hand, we make the (arbitrary) choice U

(1)

OO = U (1) OO † and U (1) V V = U (1) V V †
for the occ-occ and virtvirt blocks of U (1) , respectively, in Eq. ( 29), so that (for n ∈ occ and q ∈ virt):

U σσ (1) mn = 0 , (45b)
and:

U σσ (1) pq = 0 . ( 45c
)
From Eqs. (14b), (15a), (15b), ( 35), ( 37) and (45a), U (1) depends on G (1) , and G (1) depends on P (1) , which is itself built from U (1) . Therefore U (1) depends on itself and must be determined by an iterative solution of the first order perturbation equation.

If the self-consistent CPHF/KS procedure is initiated with a null perturbed density matrix as a starting guess P (1) guess = 0, then from Eqs. ( 12), (15a), (15b), ( 21) and (14b), R G σσ(1) guess = 0 at the first iteration. Hence, from Eqs. (45a), (45b) and (45c), R U σσ(1) guess = 0. Then, taking into account Eq. ( 35) as well, subsequent cycles of the CPHF/KS procedure will simply include nonzero elements of P (1) through the self-consistently calculated U (1) . The initial R U σσ(1) guess = 0 is carried on to subsequent cycles of the self-consistent procedure. We conclude that the diagonal spin-blocks of the first order perturbed-density matrix, as well as those of U (1) are pure imaginary:

R P σσ(1) = 0 (46)
and:

R U σσ(1) = 0 (47)
In the non-canonical formulation the occ-occ and virtvirt blocks of the Lagrange multiplier matrices are nondiagonal. The requisite expressions in first order are found by substituting Eqs. (45b) and (45c) in Eq. (40). For the occ-occ blocks

σσ (1) mn = G σσ (1) mn (48a)
and for the virt-virt blocks

σσ (1) pq = G σσ (1) pq (48b)
Substituting Eqs. ( 44)-(48b) in the second order perturbation equation, Eq. ( 41), we obtain explicit expressions for the occ-virt blocks of U (2) :

U σ σ(2) pm = 1 σ (0) p - σ(0) m σ occ n U σ σ (1) pn G σ σ(1) nm - σ virt q G σ σ (1) pq U σ σ(1) qm -G σ σ(2) pm . (49)
Using Eqs. (45b) and (45c) in the second-order orthogonality relation, Eq. ( 30), we find that:

U σ σ(2) pm = -U σσ (2) mp * . ( 50
)
Eq. ( 50) shows that the occ-virt blocks of U (2) are anti-Hermitian, which is, of course, consistent with Eq. ( 49).

To find an expression for the occ-occ and virt-virt blocks of U (2) , we return to Eq. ( 30) to obtain:

U σσ (2) mn = - 1 2 σ virt p U σ σ(1) pm * U σ σ (1) pn , (51)
for the occ-occ blocks and

U σσ (2) pq = - 1 2 σ occ m U σ σ(1) mp * U σ σ (1) mq . (52)
for the virt-virt blocks by making the (arbitrary) choice:

U σσ (2) mn = U σ σ(2) nm * , (53a)
and

U σσ (2) pq = U σ σ(2) qp * . (53b)
The above procedure can be extended to higher order to find expressions for U [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] and U [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] from the third and fourth order orthogonality conditions and perturbation equations. However, as we show below in section II H, only the occ-occ and virt-virt blocks of U [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] alone are necessary to obtain the energy through fourth-order. Using Eq. ( 31), together with Eqs. (45b) and (45c), we find:

U σσ (3) mn = U σ σ(3) nm * = - 1 2 σ virt p U σ σ(2) pm * U σ σ (1) pn + U σ σ(1) pm * U σ σ (2) pn , (54a)
for the occ-occ blocks and:

U σσ (3) pq = U σ σ(3) qp * = - 1 2 σ occ n U σ σ(2) np * U σ σ (1) nq + U σ σ(1) np * U σ σ (2) nq , (54b)
for the virt-virt blocks.

G. Order by Order Contributions to the Orbital Energies

Occupied and virtual orbital energies can be found, including first order corrections, by diagonalizing the occocc -

(0) O + (1)
OO -and virt-virt -

(0) V + (1)
V V -blocks of the matrix of perturbed Lagrange multipliers given by Eqs. (48a) and (48b). Second order contributions to the occ-occ and virt-virt blocks can be found by substituting Eqs. ( 44), (45a)-(45c) and (48a)-(48b) in the second order perturbation equation, Eq. ( 41), to obtain:

σ σ (2) nm = 1 2 σ (0) n -σ (0) m σ virt p U σ σ(1) np U σσ (1) pm + σ virt p σ(0) p -σ (0) n U σ σ(1) np U σσ (1) pm + G σ σ (2) nm , (55)
for the occ-occ blocks and

σ σ (2) qp = 1 2 σ (0) q -σ (0) p σ occ m U σ σ(1) qm U σσ (1) mp + σ occ m σ(0) m -σ (0) q U σ σ(1) qm U σσ (1) mp + G σ σ (2) qp . (56)
for the virt-virt blocks.

H. 2n+1 Rule Energy Expressions

The SCDFT total energy can be written using the density and Hamiltonian matrices of Eqs. ( 8) and ( 33):

E = 1 2 RTr (h + H) P = E (0) + λE (1) + λ 2 E (2) + . . . . (57)
By expanding these matrices, following Eqs. ( 9), (14a)-(14c) and ( 33), taking the derivative w.r.t. λ of Eq. ( 57) at order N , and evaluating the result at λ = 0, it is possible to obtain simplified expressions for the E (N ) that are consistent with Wigner's 2n+1 rule. More details on the derivation are provided in the Appendix. Here we simply quote the final results.

To write the energy contributions at low order, it proves useful to define the following matrix built from zeroth order quantities:

Θ σσ = δ σ,σ c σσ(0) † v + u AR + H σσ(0) c σσ(0) .
(58) At order zero the energy is obtained from a previous onecomponent SR calculation, which can be written in terms of Θ σσ OO as (see Appendix):

E (0) = 1 2 σ RTr [f σO Θ σσ OO ] . (59)
The first order energy contribution is (see Eqs. (A4) and (A5) of the Appendix):

E (1) = σ RTr f σO Θ σσ OV U σσ(1) V O = 0 . (60)
E (1) vanishes because Θ σσ OV is pure real and, from Eq. ( 47), U σσ(1) V O is pure imaginary. For the second-order contribution to the energy it is expedient to introduce the quantity:

Ξ σσ = c σσ(0) † u σσ SO c σ σ (0) . (61)
In terms of the above matrix, we obtain the following energy expression in second order (see Eqs. (A6)-(A12) of the Appendix):

E (2) = σσ RTr f σO Ξ σσ OV U σ σ(1) V O . ( 62
)
In third order, we find (see Eqs. (A13)-(A22) of the Appendix):

E (3) = σσ RTr G σσ (1) V V σ U σ σ (1) V O f σ O U σ σ(1) OV † - 1 2 (f σO + f σ O ) G σσ (1) OO σ U σσ (1) OV † U σ σ (1) V O , (63)
and, finally, in fourth order (see Eqs. (A23)-(A36) of the Appendix):

E (4) = σσ RTr σ(0) O (f σ O -f σO ) U σσ (2) OO † U σ σ(2) OO + G σσ (2) V O f σ O U σ σ(2) OV † + σ(0) O - σ (0) V    2U σσ (1) V O σ U σ σ (2) OO f σ O U σ σ(1) OV † -f σO U σσ (2) OV † U σ σ(2) V O    + 1 2 (f σO + f σ O )    G σσ (2) OO U σ σ(2) OO -G σσ (1) OO σ U σ σ (1) OV † U σ σ(2) V O + U σ σ (2) OV † U σ σ(1) V O    + 1 2 G σσ (2) V V σ U σ σ (1) V O f σ O U σ σ(1) OV † + 2G σσ (1) V V σ U σ σ (2) V O f σ O U σ σ(1) OV † . (64)
Thus, in accordance with Wigner's 2n+1 rule, we are able to calculate E (2) and E [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] using only first order quantities whereas E [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] depends on second order quantities as well.

I. Density Variables of the SCDFT in Orders of Perturbation Theory

Following Refs. 60 and 72, we can write the perturbed density variables (N ≥ 1) of the SCDFT in terms of the density matrix. To write them in a compact way, it is useful to introduce the notation: (66) This gives rise to the perturbation expressions: (67) for the particle number density and

X µν (r) = χ µ (r) χ ν (r) , ( 65
) Y µν (r) = χ µ (r) ∇χ ν (r) -∇χ µ (r) χ ν (r) ,
n (N ) (r) = µν R P αα(N ) µν + P ββ(N ) µν X µν (r) ,
m (N ) x (r) = µν R P βα(N ) µν + P αβ(N ) µν X µν (r) , (68a)
m (N ) y (r) = µν I P βα(N ) µν -P αβ(N ) µν X µν (r) , (68b) m (N ) z (r) = µν R P αα(N ) µν -P ββ(N ) µν X µν (r) . (68c)
for the magnetization. The perturbed orbital-current density reads

j (N ) (r) = - 1 2 µν I P αα(N ) µν + P ββ(N ) µν Y µν (r) , (69)
while the perturbed spin-current densities are

J (N ) x (r) = - 1 2 µν I P βα(N ) µν + P αβ(N ) µν Y µν (r) , (70a)
J (N ) y (r) = 1 2 µν R P βα(N ) µν -P αβ(N ) µν Y µν (r) , (70b)
J (N ) z (r) = - 1 2 µν I P αα(N ) µν -P ββ(N ) µν Y µν (r) . (70c)
Substituting Eq. ( 46) in Eqs. ( 67) and (68c), we find that the first order perturbed particle-number density and z-component of the magnetization vanish:

n (1) (r) = 0 , (71a)
m (1) z (r) = 0 . ( 71b
)
From Eqs. ( 60)-( 64), we conclude that the perturbed n (N ) and m

(N ) z only affect the total energy in fourth order. On the other hand, all of the other perturbed density variables m

(N ) x , m (N ) y , j (N ) , J (N )
x , J (N ) y and J (N ) z have non-vanishing contributions for N = 1, and therefore affect the total energy in second order.

Eq. ( 71) or, equivalently, Eq.( 46), has consequences with regard to which electron-electron repulsion terms need to be evaluated in a calculation. For instance, substituting Eq. ( 46) in Eq. (15a), we find that the first order Coulomb matrix elements vanish:

C σσ (1) µν = 0 , (72)
The perturbed Coulomb potential, therefore, does not affect the total energy until fourth order. Furthermore, for collinear formulations of the xc potential, substituting Eqs. ( 23) and ( 71), in Eqs. ( 21) and ( 22), we obtain in first order:

V σσ (1) µν = 0 , (73)
and in second order:

V σσ (2) µν = δ σ,σ σ µ| ∂ 2 F xc n α(0) , n β(0) ∂n σ(0) ∂n σ (0) n σ (2) |ν . (74)
Therefore, the perturbed collinear xc potential also does not affect the total energy until fourth order.

III. STATEMENT OF THE PROCEDURE

Our CPKS calculation procedure is described in the following. In this context a quantity obtained at cycle k is denoted as, for example, as P (1) (k):

1. Solve the reference one-component SR KS-DFT problem, which yields, P (0) , c (0) and (0) .

2. Set k = 1.

3. Set P (1) (k) = 0 and H (1) (k) = u SO .

4. Start first order CPKS procedure:

(a) Calculate G (1) (k) using H (1) (k) and Eq. ( 37).

(b) Calculate U

(1)

OV (k) from Eq. (45a). (c) Calculate E (2) (k) and E [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] (k) from Eqs. (62) and (63).

(d) If k > 1, check for convergence on E (2) (k) and E (3) (k):
i. If convergence on the energy is reached, exit to step 5, otherwise continue to step 4e.

(e) Calculate P (1) (k + 1) from Eq. ( 35).

(f) Calculate H (1) (k + 1) from Eq. (14b), using subset of electron-electron terms needed through order three in the energy, as determined from Table I.

(g) Set k → k + 1 and go back to step 4a.

5. If occupied orbital energies are desired at first order (i.e. without

(2)

OO ), calculate (0)
O +

(1)

OO from Eq. (48a) and diagonalize the resulting matrix.

6. If virtual orbital energies are desired at first order (i.e. without

(2)

V V ), calculate (0)
V +

(1)

V V from Eq. (48b) and diagonalize the resulting matrix. 7. If E [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] and/or P (2) and/or (2) are desired, set k = 1. Otherwise, exit to step 12.

8. Set P (2) (k) = 0 and H (2) (k) = 0.

9. Start second order CPKS procedure:

(a) Calculate G (2) (k) using H (2) (k) and Eq. ( 37).

(b) Calculate U (2) (k) from Eqs. ( 49)-( 53).

(c) Calculate E (4) (k) from Eq. ( 64).

(d) If k > 1, check for convergence on E (4) (k):

i. If convergence on the energy is reached, exit to step 12, otherwise continue to step 9e.

(e) Calculate P (2) (k + 1) from Eq. ( 36).

(f) Calculate H (2) (k + 1) from Eq. (14c).

(g) Set k → k + 1 and go back to step 9a.

10. If occupied orbital energies are desired at second order calculate

(0) O + (1)
OO +

(2)

OO from Eq. ( 55) and diagonalize the resulting matrix.

11. If virtual orbital energies are desired at second order calculate (0)

V +

(1)

V V + (2)
V V from Eq. ( 56) and diagonalize the resulting matrix.

End

TABLE II: SOC contributions to the total energy (in E h ) as calculated by the CPKS method, using the PBE0 functional, for the homonuclear diatomic molecules of the halogen series. Our results were obtained using the fully relativistic shapeconsistent RECPs, with associated basis sets, of the Columbus group available at Ref. 73. Seven electrons were included in the valence space of each atom. The first row reports the reference 2c-SCF result, while subsequent rows report the difference through N th order between the CPKS and reference 2c-SCF values (e.g. ∆2 = E (2) -ESOC(2c-SCF))

F2 Cl2 Br2 I2 2c -7.03×10 -05 -2.94×10 -04 -7.13×10 -03 -3.82×10 -02
∆2 +7.59×10 -08 -2.81×10 -07 -1.18×10 -04 -1.13×10 -03 ∆3 -7.95×10 -10 -4.00×10 -08 -3.33×10 -05 -6.58×10 -04 ∆4 -4.25×10 -10 +5.80×10 -09 -2.59×10 -07 +2.37×10 -04

IV. COMPUTATIONAL DETAILS

Reference one-component SR KS-DFT calculations were performed with a developer's version of the Crys-tal17 code, [START_REF] Dovesi | Quantum-Mechanical Condensed Matter Simulations with CRYSTAL[END_REF] in which we also implemented our CPH-F/KS treatment of SOC. The CPHF/KS cycles were initiated after convergence of the reference 1c-SCF calculations within a criterion of 1×10 -12 Hartree a.u. (E h ) on the total energy. The same tolerance was used as a convergence criterion for the perturbed energies obtained from the first-and second-order CPHF/KS procedures.

As was the case in Part I, [START_REF] Desmarais | Perturbation Theory Treatment of Spin-Orbit Coupling. Part I: Double Perturbation Theory from a Single-Reference Initial Approximation[END_REF] our calculations were performed using RECPs of the energy-and shape-consistent kinds, with their associated basis sets, from the libraries available in Refs. 75 and 73. Experimental bond lengths of 1.42, 2.00, 2.28, 2.67 and 3.00 Å were used for F 2 , Cl 2 , Br 2 , I 2 and At 2 , respectively. Calculations were also performed on the halogen hydride molecules HBr, HI, HAt and HTs, using bond lengths of 1.41, 1.61, 1.72 and 1.72 Å , with the basis set of Gatti et al. 76 for the hydrogens.

The xc contribution was calculated by numerical quadrature using Gauss-Legendre radial and Lebedev angular point distributions, [START_REF] Towler | Density functional theory in periodic systems using local Gaussian basis sets[END_REF][START_REF] Lebedev | Quadratures on a sphere[END_REF][START_REF] Lebedev | Spherical quadrature formulas exact to orders[END_REF] with the quadrature weights proposed by Becke. [START_REF] Becke | A multicenter numerical integration scheme for polyatomic molecules[END_REF] We used a pruned grid consisting of 99 radial points and 1454 angular points (keyword XXLGRID in the Crystal17 manual). [START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF] The Xc-Fun library 82 was employed for taking the xc functional second derivatives required for the second-order CPKS procedure (i.e. for a calculation of the xc contribution to E [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] ).

Energies calculated by the CPHF/KS method were compared to those obtained from the 2c-SCF implementation in the Crystal code. [START_REF] Desmarais | Adiabatic connection in spin-current density functional theory[END_REF][START_REF] Desmarais | Spin-orbit coupling in periodic systems with broken time-reversal symmetry: Formal and computational aspects[END_REF][START_REF] Desmarais | Development of Tools for the Study of Heavy-Element Containing Periodic Systems in the CRYSTAL Code and their Application[END_REF][START_REF] Desmarais | Spin-orbit coupling from a two-component self-consistent approach. I. Generalized Hartree-Fock theory[END_REF][START_REF] Desmarais | Fundamental Role of Fock Exchange in Relativistic Density Functional Theory[END_REF] The CPKS calculations on the closed-shell (F 2 , Cl 2 , Br 2 , I 2 , At 2 , HBr, HI, HAt, HTs) and open-shell (F - 2 , Cl - 2 , Br - 2 , I - 2 , At - 2 ) systems were performed with the standard collinear formulation of the SVWN5, BHandH (i.e. SVWN5 with a 0.5 fraction of exact exchange), PBE and PBE0 xc functionals. [START_REF] Slater | A Simplification of the Hartree-Fock Method[END_REF][START_REF] Vosko | Accurate Spindependent Electron Liquid Correlation Energies for Local Spin Density Calculations: a Critical Analysis[END_REF][START_REF] Becke | A new mixing of Hartree-Fock and local density-functional theories[END_REF][START_REF] Perdew | Generalized Gradient Approximation Made Simple[END_REF][START_REF] Adamo | Toward Reliable Density Functional Methods without Adjustable Parameters: the PBE0 Model[END_REF]

V. COMPUTATIONAL RESULTS FOR THE HALOGEN SERIES OF DIATOMIC HOMONUCLEAR AND HYDRIDE MOLECULES

The CPKS approach corresponds to a perturbation theory treatment of a 2c-SCF calculation for SOC. Below we compare the two with regard to i) total energy ii) orbital energy levels and iii) spatial distribution of SCDFT density variables.

A. Convergence of Perturbation Series for the Total Energy

In this sub-section we compare total energies. Our calculations were performed on halogen diatomic molecules, because of the large contribution of SOC to their total energies as well as the availability of many sets of RECPs and associated valence basis sets. The tables provide results obtained with the PBE0 functional; similar results for the SVWN5, BHandH, and PBE functionals, as well as the HF method, are provided in the electronic supporting information (ESI). The first row of Table II reports the SOC contribution to the total energy (in E h ), as calculated using the reference 2c-SCF approach with the shape-consistent RECPs of Ref. 73 and 7 electrons in the valence space of each atom. Lower rows report the energy differences, through N th order ∆ N (N = 2, 3 or 4), for the CPKS approach w.r.t. the reference 2c-SCF. The ∆ N are found to be more than one order of magnitude smaller than the SOC contribution itself (in absolute value), indicating that even just through second-order the CPKS energy is an excellent approximation to the 2c-SCF value. Moreover, through fourth-order the ∆ N are monotonically decreasing (in absolute value). Similar results are obtained for the same systems using other functionals in Table S1 of the ESI.

Br - 2 I - 2 I - * 2 At - 2 2c -3.28×10 -03 -6.14×10 -03 -5.71×10 -03 -5.72×10 -02 ∆2 -3.29×10 -05 -2.56×10 -05 -2.57×10 -05 -3.09×10 -03 ∆3 +3.18×10 -06 -1.17×10 -05 -1.12×10 -05 -1.15×10 -03 ∆4 -4.61×10 -07 -1.71×10 -07 +6.99×10 -08 +8.95×10 -04
In Table III we report the results of similar calculations, which were performed instead with the energyconsistent RECPs and associated basis sets from Ref. 75. These potentials were available for the systems Br 2 , I 2 and At 2 , again including 7 electrons in the valence space of each atom. The reported ∆ N indicate that, again, the second-order CPKS values represent excellent approximations to the 2c-SCF energies, except for At 2 which requires the third-order correction to reduce the deviation from roughly 25 to 5 %. We can see from Table III that At 2 is significantly more challenging than the other lighter-element systems. In this case the perturbation series is no longer monotonically convergent at fourthorder. This suggests the desirability of using a thirdorder treatment in general, for the energy, especially considering the reduction in computational effort since it is only necessary to solve the first-order perturbation equation. Moreover a calculation through third-order in the energy requires only the exact-exchange contribution to the electron-electron potential, as may be recalled from Table I. These conclusions are unchanged when calculations are performed with other xc functionals, as can be seen from Table S2 of the ESI. Finally, we note that in the case of Br 2 , for some functionals (namely PBE, PBE0 and SVWN5), the ∆ N are not monotonically convergent at fourth-order (because

|∆ 4 | > |∆ 3 |), however the E (N ) are monotonic, as |E (2) | > |E (3) | > |E (4) |.
Tables IV and V provide results for systems related to those of Tables II andIII, by the addition of one electron to the lowest unoccupied molecular orbital. These results reenforce the conclusions that were drawn from tables II and III. In fact, the deviations from the 2c-SCF reference are now monotonically decreasing through fourth-order in all cases. Similar results are provided in Tables S3 and S4 of the ESI with other functionals and the HF method. III, but now for the halogen hydride molecules (instead of the homonuclear halogen diatoms). For the super-heavy element Ts, 25 electrons are treated in the valence space with the associated double-zeta valence basis set.

HBr HI HAt HTs 2c -1.59×10 -03 -3.11×10 -03 -2.80×10 -02 -7.23×10 -01

∆2 -1.57×10 -05 -1.29×10 -05 -8.70×10 -04 -1.01×10 -01 ∆3 +1.59×10 -06 -5.33×10 -06 -4.14×10 -04 -2.67×10 -02 ∆4 -2.05×10 -07 -1.38×10 -06 -3.52×10 -05 +1.05×10 -02

Table VI provides results that, again, use the energyconsistent potentials, this time for the halogen hydride molecules. The hydrides have a wider HOMO-LUMO gap than the homonuclear diatomics and, thus, might be considered more suitable for our perturbation treatment. Indeed, in comparing Tables III and VI we do see improved agreement in every common instance.

Results are also provided for Tennessine hydride (HTs, Ts being the super-heavy element No. 117). In that case, we used the ECP92MDFQ potential available from Ref. 75 with the associated double-zeta valence basis set. Our perturbation theory treatment appears well-behaved also for this p-block super-heavy element system. Similar results are provided in Table S5 of the ESI by using other xc functionals and the HF method.

In summary, the proposed CPKS treatment leads to a perturbation series through 3rd order in the total energy that converges monotonically to a value close to the 2c-SCF result for the neutral closed-shell, and singly negative open-shell systems considered here (i.e., F 2 to At 2 , F - 2 to At - 2 , HBr to HAt, as well as HTs). In one instance, At 2 , this monotonically convergent behavior breaks down in 4th order. For At 2 , monotonicity is restored by scaling the SOC operator (scaling factor ≤ 0.68). Simultaneously, the ratio of the second-order energy to the HOMO-LUMO gap (= 2.7 eV) is reduced to less than 0.43. By examining a plot of |E (2) |/E gap for the whole series from F 2 to At 2 , (see Figure 1) we see that this ratio increases monotonically with atomic number and, therefore, speculate that it can be used as a test for, at least, asymptotic convergence through 3rd order in the total energy.

Of course, other types of system may not exhibit the same convergence behavior as found in this initial investigation. In fact, our preliminary studies indicate that the SOC perturbation series does not converge in low order for the singly positive homonuclear halogen diatomics (even for F + 2 ). In that case, the significant multi-reference character of the ground state caused by quasi-degeneracies leads to failure of the CPKS procedure. For such cases (and also to improve the treatment of molecules like At 2 with a small bandgap), we have begun to extend our single-reference approach so as to include low-lying electronic states through quasi- degenerate perturbation theory (QDPT). [START_REF] Kirtman | Simultaneous calculation of several interacting electronic states by generalized Van Vleck perturbation theory[END_REF][START_REF] Kirtman | Variational form of Van Vleck degenerate perturbation theory with particular application to electronic structure problems[END_REF] FIG. 3: KS eigenvalue spectrum calculated with the energyconsistent RECPs and PBE0 functional for the specified hydrides. The eigenvalues are also provided in Tables S9-S11 of the ESI.

B. Convergence of Perturbation Series for Orbital Energies

We now provide the KS eigenvalue spectrum calculated by our perturbation treatment, at first and second order, for a representative set of the most challenging systems considered above. Energy levels for I 2 , I - 2 and At - 2 are given in Figure 2 and for the halogen hydrides HI, HAt and HTs in Figure 3. The orbital energies are plotted on the vertical axis as differences with respect to the HOMO, i.e. -HOMO , in atomic units. Panels (b) and (c) of Figure 2 show the spin splitting of the energy levels for the negatively charged open-shell ions I - 2 and At - 2 due to SOC. The plots show that the eigenvalue spectrum is nearly exactly reproduced, at just first order, except for the two most challenging cases represented by At - 2 and HTs. For both of these the second order treatment yields significant improvement, leading to close agreement with the reference 2c-SCF calculations, for the virtual as well as the occupied energy levels.

C. Convergence of Perturbation Series for SCDFT Density Variables

Finally, we test the ability of our perturbation theory approach to reproduce the density variables of SCDFT. Results are presented for the halogen homonuclear diatoms, rather than hydrides, because we have already shown in sections V A and V B that the agreement between the perturbation theory and 2c-SCF approaches is superior for the hydrides, which means that the diatoms represent a more challenging test.

Insight into the physical meaning of the density variables of SCDFT is provided by the corresponding continuity equations. [START_REF] Vignale | Current-and spin-densityfunctional theory for inhomogeneous electronic systems in strong magnetic fields[END_REF][START_REF] Bencheikh | Spin-orbit coupling in the spin-currentdensity-functional theory[END_REF] In the absence of SOC, but in the presence of an external magnetic field, these equations (Eqs. 6.8a and 6.8b of Ref. 54) show that the orbitalcurrent density j can be interpreted as the velocity field for n, while the ith Cartesian component of the spincurrent density J i is the velocity field for m i . In the presence of SOC, however, both the orbital-and spincurrent densities j and J x , J y , J z couple all variables n, m x , m y and m z .

A difference plot of the density variables as calculated from the 2c-SCF and the order-by-order CPKS methods provides a visual representation of the ability of the perturbation theory approach to reproduce the eight different spin blocks (see Eqs.( 67) -(70c)) of the SOC contribution to the density matrix. Figure 4 provides such a representation for the system I - 2 obtained using the PBE0 functional and the energy-consistent RECP. The columns of panels represent (from left to right) the density variables n, m, j, J x , J y and J z . The first row of panels (denoted as ∆2c) is the SOC contribution to the density variables, as calculated by taking the difference between the 2c-SCF and 1c-SCF distributions (the 1c-SCF distributions is non-zero only for n and m z ). Subsequent rows denote the differences between the first and second order perturbation theory and the 2c-SCF values. All values are in atomic units, with a common color for each density variable, and in which the intensity of the coloration correspond to associated values provided with the color bar on the right of individual columns.

FIG. 4: Spatial distribution of differences of the SCDFT density variables for the system I - 2 in the xz plane, with the molecular axis along z, as obtained using the PBE0 functional and the energy-consistent RECPs. The quantities in the first row are the differences ∆2c = 2c -1c where 2c denotes the distribution of the density variable calculated by 2c-SCF and 1c denotes a 1c-SCF calculation. Subsequent rows report ∆PT1 = PT1 -2c and ∆PT2 = PT2 -2c, where PT1 (PT2) denote the density variables obtained from a first-(second-)order perturbation theory treatment. From left to right the columns show the particle number density n, the magnetization m, the orbital-current density j and the three spin-current densities Jx, Jy and Jz. When vector fields are plotted, the length and direction of the arrows represents the projection in the xz plane. The color intensity represents the magnitude of the three dimensional vector. All quantities are plotted in atomic units.

Figure 4 shows that all three spin-current densities are well-reproduced by the first-order perturbation theory treatment. On the other hand, the first-order perturbed n and m z vanish (see Eqs. (71a), (71b); as a result there are important differences in first order for the particle number density, magnetization and orbital-current density (which couples strongly with the magnetization).

Thus, for open-shell systems with non-vanishing magnetization, a second order treatment is necessary. It is also sufficient as can be seen from the bottom row of Figure 4 for the case of I - 2 . Figure 5 provides a similar plot, this time for the system At - 2 . In this case, we again see that first-order perturbation theory is able to reproduce the spin-current densities of the reference 2c-SCF. However, in this more challenging case, visible differences remain in the distributions of n, m and j, even with a second-order treat-ment. Although the error in the total energy (cf. Table V) is quite small, it is much more significant in the distribution of the n, m and j density variables.

VI. COMPARISON OF COMPUTATIONAL REQUIREMENTS BETWEEN 2C-SCF AND CPKS

Although we have not yet developed a fully optimized CPHF/CPKS computer code, it is possible to compare the timing of this procedure with that of the alternative 2c-SCF approach analytically. Our analysis shows that in the case of large systems (as defined below) the number of floating point operations required for the former is more than an order of magnitude the lesser of the two. We start by analyzing the computational scaling with system size for the 2c-SCF approach. One cycle of the latter is composed of three computationally intensive steps: i) evaluation of two-electron integrals ii) transformation of Hamiltonian matrix from atomic orbital (AO) to molecular orbital (MO) basis iii) diagonalization of MO Hamiltonian matrix. [START_REF] Szabo | Modern quantum chemistry: introduction to advanced electronic structure theory[END_REF] In the following discussion of these three aspects, the terms "small system" and "large system" are utilized to describe systems for which cpu time is dominated by step i), or by steps ii) and iii), respectively.

In step i), the calculation of bielectronic integrals formally scales as the fourth power of the number of AOs, N AO , but may be reduced to linear scaling for large sys-tems by using standard approaches for integral screening, such as the Tolinteg strategy used in the Crystal code, [START_REF] Dovesi | CRYSTAL17 User's Manual[END_REF][START_REF] Pisani | Hartree-Fock ab initio treatment of crystalline systems[END_REF] Schwarz integral estimates, and others. [START_REF] Ochsenfeld | Linearscaling methods in quantum chemistry[END_REF] On the other hand, steps ii) and iii) scale according to the third power of N AO , although, with a smaller prefactor. Thus, for small systems (typically containing ten to a few dozen atoms or less) the cpu time is dominated by step i) while, for large systems, (typically containing tens to hundreds of atoms or more), it is dominated by steps ii) and iii). Note that the exact size thresholds are variable, depending, for instance, on the number of processors being used. [START_REF] Erba | Large-scale condensed matter DFT simulations: performance and capabilities of the CRYSTAL code[END_REF] The conversion in Step ii) from H to H = c † H c is done via two successive 3-index transformations. This process scales as (3/2) × 4(a + m) (2N AO ) [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] = 48(a + m)N 3 AO (a and m denote, respectively, floating point addition and multiplication. Here the multiplicative factor of 2 is for spin; the factor of 4 = 2 × 2(a + m) accounts for the fact that each matrix element is complex; and, finally, the factor of 3/2 = 1 + 1/2 (rather than 2) occurs because the second transformation results in an Hermitian matrix, thereby halving the number of independent elements.

The diagonalization in Step iii) can be achieved, for example, using the Jacobi algorithm [START_REF] Press | Numerical recipes 3rd edition: The art of scientific computing[END_REF] , which is the procedure currently used in our 2c-SCF code. At each Jacobi iteration, a "pivot" element (the largest off-diagonal element) is selected and a 2 × 2 rotation of the basis is performed. The updating of the Hamiltonian matrix requires a total of (4N AO )×(4a+4m) operations where the factor of 4N AO arises because 2 columns (or rows) containing 2N AO elements must be updated and the factor of 4(a + m) is due to the complex arithmetic involved. The eigenvectors (MO coefficients) are given by the product of Jacobi rotation matrices. As it was for the Hamiltonian, updating the eigenvector matrix involves another (4N AO ) × (4a + 4m) operations. After the update has been performed, a new pivot element is chosen until all of the (1/2)(2N AO ) 2 off-diagonal matrix elements have been considered. Thus, one Jacobi "sweep" consists of (1/2)(2N AO ) 2 rotations. The procedure is repeated for K sweeps until all the off-diagonal elements are smaller than a pre-set criterion. In our 2c-SCF calculations we have found that, typically, K >> 1 on the first SCF cycle, but << 1 on subsequent cycles. This gives a total scaling of the 2c-SCF diagonalization per cycle for large systems of 128K ×(a+m)×N [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] AO where K is an average over all cycles.

Next we turn to the CPKS treatment. For small systems some savings are obtained through third-order in the energy by the fact that only SOC and Fock exchange integrals need be calculated (no SR mono-electronic, Coulomb or exchange-correlation integrals are necessary, see Table I). This is particularly significant if non-hybrid functionals are used. On the other hand, no savings is expected for small systems if E [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] is needed.

For large systems the story is quite different. Diagonalization is avoided and the calculation of the first order coefficients scales as N occ × N virt × N AO , with N occ << N AO . Thus, there is just one significant timeconsuming step that scales as N 3 AO , namely the calculation of G σσ (N ) = c σσ(0) † H σσ (N ) c σ σ (0) with N = 1 (and N = 2 if E (4) is needed). This closely resembles step ii) of the 2c-SCF procedure. The only difference is that, here, the relevant orbital coefficients are the 1c-SCF, rather than 2c-SCF, coefficients. They are real rather than complex and, therefore, the factor of 4 in the 2c-SCF expression above is halved.

In furthering the analysis it is convenient to treat the spin blocks G σσ (N ) individually. Then, the factor (2N AO ) 3 in the 2c-SCF expression above reduces to N 3 AO for each spin-block. Within the diagonal spin-blocks of G (1) the real matrix elements vanish (cf. Eq.( 38)), so we must divide by another factor of 2, to obtain a total of (3/2)(a + m)(N AO ) 3 operations needed to calculate each Hermitian diagonal spin-block. Unlike the diagonal spin-blocks, the individual off-diagonal spin blocks are, in general, non-Hermitian, which means that the factor of 3/2 is replaced by 2 × 2 = 4. There are two diagonal spin-blocks and one-independent off-diagonal spin-block (the other is its Hermitian conjugate). Hence, the total number of operations is [2×(3/2)+4] = 7(a+m)(N AO ) 3 . For G (2) there is an extra multiplicative factor of 2 in the formula for the diagonal spin-blocks (since the real matrix elements do not vanish) and, then, the number of operations becomes 6 + 4 = 10(a + m)(N AO ) 3 .

Comparing the combination of steps ii) and iii) for 2c-SCF with those above for CPKS, we conclude that, through E [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] , the CPKS approach yields a saving factor for the number of floating point operations (N flops ) of (48 + 128K )/7 for large systems. From our experience we expect K to be typically greater than 1/6, which would lead to an order of magnitude savings in cpu time. If E [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] is needed, then the savings reduces to (48 + 128K )/(7 + 10).

We note that even further savings of cpu time may be achievable since: 1) no diagonalization is required for CPKS so that it is unnecessary to construct irreducible representations of the double-group in order to fully exploit symmetry; and 2) the number of iterations required for convergence in the CPKS method is likely to be less than for 2c-SCF. As regards memory usage, a saving factor of 2 may be expected vis-à-vis 2c-SCF. This is because the diagonalization step of the 2c-SCF necessitates simultaneously storing both the upper and lower triangular halves of the Hamiltonian and molecular orbital coefficient matrices. On the other hand, for the CPKS calculation, the Lower triangular elements may always be determined by (anti-)Hermiticity.

We expect to develop an optimized implementation of the CPKS procedure within the CRYSTAL code in the near future and will, then, present actual timings for comparison with 2c-SCF in conjunction with new applications.

VII. CONCLUSIONS

We have presented a non-canonical coupled perturbed Kohn Sham DFT/HF perturbation treatment of spinorbit coupling (SOC). Our treatment, based on an initial relativistic effective core potential approximation (RECP) accounts for spin-and orbital-current densities as well as the particle number density and magnetization. Perturbation expressions that satisfy the 2n + 1 rule are given through fourth-order. Working procedures for obtaining orbital energies that are correct through second order, along with magnetization and density variables, are also provided.

Tests on the halogen homonuclear diatomic and hydride molecules, including 6p and 7p elements, show that the proposed approach is capable of providing property values through second-or third-order that essentially match those from 2c-SCF calculations. The only exception was for At - 2 , in which case second-order differences remained in the particle-number density, magnetization and orbital-current density.

Our spin-current density functional perturbation theory treatment provides an efficient means of adding SOC effects to the SR KS-DFT/HF approximation. In particular, diagonalization in the two-component spinor basis is avoided, leading to saving factors on the number of re-quired floating point operations that may exceed an order of magnitude. We intend to adapt this treatment for solid state calculations in the near term. Applications to systems with small gaps or, in general, with strong multi-reference character would require the use of quasidegenerate perturbation theory (QDPT) as in ensemble DFT [START_REF] Englisch | Exact Density Functionals for Ground-State Energies. I. General Results[END_REF][START_REF] Englisch | Exact Density Functionals for Ground-State Energies II. Details and Remarks[END_REF][START_REF] Gross | Density-functional theory for ensembles of fractionally occupied states. I. Basic formalism[END_REF][100][101][102][103][104][105][106][107][108][109][110][111][112][113][114] , which we leave for further in the future.

Supporting Information

See Supporting Information at URL for a comparison a benchmark of the total energies calculated from our perturbation theory approach with the SVWN5, BHandH and PBE functionals, as well as the HF method. Also provided are tables containing the orbitals energies for Figures 2 and3.

We obtain the second and third lines of Eq. (A6) by combining the second and fourth lines of Eq. (A1c) while using the expression for the first order Hamiltonian matrix from Eq. (14b).

Let us now first consider the contribution in the first line of Eq. (A6). Using Eq. ( 36) we obtain: To get the last line of Eq. (A10), we have used Eq. (51).

Let us now consider the term in the second line of Eq. (A6). Proceeding as in Eqs. (A4) and (A5), we obtain: where, in the last passage, we have used Eq. (45a). If we now insert Eqs. (A10) and (A11) in Eq. (A6) the second-order energy reduces to:

E (2) =
1 2 σσ RTr u σσ SO P σ σ(1) . (A12)

Finally, substituting Eq. ( 35) and (61) in Eq. (A12), we obtain Eq. ( 62).

Third Order Energy

From Eq. (A1d) and Eqs. (14a)-(14c) we find that:

E (3) = σ
RTr H σσ(0) P σσ(3) + σσ RTr u σσ SO P σ σ(2) + σσ RTr V σσ (1) -aK σσ (1) P σ σ(2) . (A13)

The first term of the first line of Eq. (A13) is obtained by combining the first and fourth lines of Eq. (A1d) and using the relation (for N, M = 0, 1, 2, 3, 4):

σσ RTr K σσ (N ) P σ σ(M ) = σσ RTr K σσ (M ) P σ σ(N ) , (A14) as well as: σσ RTr V σσ (N ) P σ σ(M ) = σσ RTr V σσ (M ) P σ σ(N ) .

(A15)

The second line of Eq. (A13) is obtained in a similar way, by combining the second and third lines of Eq. (A1d) and using: σσ RTr C σσ (N ) P σ σ(M ) = σσ RTr C σσ (M ) P σ σ(N ) , (A16) Substituting Eqs. (14b) and (14c) in Eq. (A13) we obtain the simplified expression:

E (3) = σ
RTr H σσ(0) P σσ [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] + σσ RTr H σσ (1) P σ σ(2) , (A17)

Let us now consider the term in the first line of Eq. (A17). Using Eq. (A2), as well as Eqs. ( 37) and (39), we find: × f σ U σσ(3) † + U σσ [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF]

f σ + σ U σσ (2) f σ U σ σ(1) † + U σσ (1) f σ U σ σ(2) † . ( A18
)
Using the third order orthogonality condition of Eq. (31) in Eq. (A18) and the fact that the diagonal matrices σ(0)

and virt-virt blocks, as:

+ σ V σ U σσ (3) V O f σ O U σ σ(1) OV † + U σσ (1) V O f σ O U σ σ(3) OV † + U σσ (2) V O f σ O U σ σ(2) OV † + σ O σ U σσ (2) OO f σ O U σ σ(2) OO † . (A26)
Moving onto the term H σσ (1) P σ σ [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] in Eq. (A23) and using the expression for the third order density matrix from Eq. (A2), we find:

× f σ U σ σ(3) † + U σ σ(3) f σ + σ U σ σ (2) f σ × U σ σ(1) † + U σ σ (1) f σ U σ σ(2) † . (A27)
By taking advantage of the first order perturbation equation, Eq. ( 40), Eq. (A27) can be re-expressed as:

× σ f σ U σ σ(3) † -σ U σσ (1) f σ U σ σ(3) † + σσ (1) f σ U σ σ(3) † + U σσ (1) σ U σ σ(3) f σ -σ U σσ (1) U σ σ(3) f σ + σσ (1) U σ σ(3) f σ + G σσ (1) σ U σ σ (2) f σ U σ σ(1) † + U σ σ (1) f σ U σ σ(2) † . ( A28
)
We now write Eq. (A28) in terms of occ-occ occ-virt and virt-virt blocks and use Eq. (48a), to find:

f σO + G σσ (1) V O σ U σ σ (2) OO f σ O U σ σ(1) OV † + G σσ (1) OV σ U σ σ (1) V O f σ O U σ σ(2) OO † + G σσ (1) V V σ U σ σ (2) V O f σ O U σ σ(1) OV † + G σσ (1) V V σ U σ σ (1) V O f σ O U σ σ(2) OV † . ( A29
)
Then we substitute Eqs. (A26) and (A29) into Eq. (A23). In doing so, terms 1, 2, 5 and 6 of Eq. (A26) cancel with terms 1, 5, 4 and 2 of Eq. (A29), respectively, leading to:

E (4) = σσ RTr -σ O U σσ (2) OO † U σ σ(2) OO f σO -σ O U σσ (2) OV † U σ σ(2) V O f σO + σ V U σσ (2) V O f σ O U σ σ(2) OV † + σ O U σσ (2) OO f σ O U σ σ(2) OO † + G σσ (1) OO f σ O U σ σ(3) OO † + G σσ (1) OO U σ σ(3) OO f σO + G σσ (1) V O σ U σ σ (2) OO f σ O × U σ σ(1) OV † + G σσ (1) OV σ U σ σ (1) V O f σ O U σ σ(2) OO † + G σσ (1) V V σ U σ σ (2) V O f σ O U σ σ(1) OV † + G σσ (1) V V σ U σ σ (1) V O f σ O U σ σ(2) OV † + 1 2
σσ RTr H σσ (2) P σ σ(2) . (A30) Substituting Eq. (54a) for the occ-occ blocks of U [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] in

  FIG. 1: Ratio of |E (2) | over the HOMO-LUMO gap (plotted on a log scale) for the series F2 to At2 (black circles), for the calculations of Table II (F2 and Cl2), and Table III (Br2, I2 and At2), with the large-core Columbus and Stuttgart RECPs and the PBE0 functional. The red horizontal line denotes the value |E (2) |/Egap = 0.43, and the blue triangle and asterisk (also shown in close-up panel) are for At2 with the SOC operator scaled down by a factor of 0.68 and 0.69, respectively.

2c - 5 .

 5 96×10 -05 -2.37×10 -04 -5.65×10 -03 -3.11×10 -02 ∆2 +1.04×10 -07 +2.51×10 -07 +4.70×10 -06 +6.96×10 -04 ∆3 +3.13×10 -10 -9.13×10 -09 -8.12×10 -06 -6.10×10 -05 ∆4 -1.03×10 -10 -8.76×10 -10 -4.04×10 -07 +5.69×10 -07

FIG. 2 :

 2 FIG. 2: KS eigenvalue spectrum calculated with the energyconsistent RECPs and PBE0 functional for I2, I - 2 and At - 2 . The eigenvalues are provided also in Tables S6-S8 of the ESI.

FIG. 5 :

 5 FIG.5: Same as Figure4, but this time for At - 2 with the PBE0 functional and the energy-consistent RECPs.

σ

  RTr H σσ(0) P σσ(2) = = σ RTr c σσ(0) † H σσ(0) c σσ(0) f σ U σσ(2) † + U σσ(2) f σ + σ U σσ (1) f σ U σ σ(1) † . (A9)Since c σσ(0) † H σσ(0) c σσ(0) = σ(0) it follows from the Hermiticity of U σσ(2) OO (see Eq. (53)), as well as Eqs. (45a)-(45b), that: σ RTr H σσ(0) P σσ(2)

σσ

  RTr H σσ (1) P σ σ(1)

σ

  RTr H σσ(0) P σσ(3) = σ RTr σ(0)

σ

  RTr H σσ(0) P σσ[START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF]

σσ

  RTr H σσ (1) P σ σ(3) = σσ RTr G σσ (1)

σσ

  RTr H σσ (1) P σ σ(3) = σσ RTr U σσ (1)

σσ

  RTr H σσ (1) P σ σ[START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF]

  TableIsummarizes the minimum order of the energy affected by the various terms of the electron-electron potential, both for closed-shell and open-shell systems. As can be seen from the table, in most cases a calculation can be carried out through third order in the energy by including only the exact-exchange contribution to the electron-electron potential.

TABLE I :

 I Minimum order of the total energy affected by the various terms of the electron-electron potential. V is the SCDFT xc potential in an LDA or GGA treatment (either from a collinear or non-collinear formulation), Ĉ is the Coulomb potential and K is the exact-exchange potential.

		V		Ĉ	K
		Collinear Non-Col.
	Closed-Shell	4	4	4	2
	Open-Shell	4	2	4	2

TABLE III :

 III Same as TableII, but using instead the energyconsistent RECPs and associated basis sets of the Stuttgart group, available at Ref. 75, with 7 electrons in the valence space of each atom. The asterisk denotes a calculation on the I2 molecule using the larger valence basis set of Ref.91

TABLE IV :

 IV Same as table II, but now for the open-shell radicals, in which one electron has been added to the lowest unoccupied orbital.

TABLE V :

 V Same as table III, but now for the open-shell radicals in which one electron has been added to the lowest unoccupied molecular orbital.

TABLE VI :

 VI Same as table

* Electronic address: jacqueskontak.desmarais@unito.it † Electronic address: kirtman@chem.ucsb.edu 1 Chang, C.; Pelissier, M.; Durand, P. Regular twocomponent Pauli-like effective Hamiltonians in Dirac theory. Phys. Scr. 1986, 34, 394. 2 Buenker, R. J.; Chandra, P.; Hess, B. A. Matrix representation of the relativistic kinetic energy operator: Twocomponent variational procedure for the treatment of

Acknowledgements J.K.D. is grateful to the National Science and Engineering Research Council of the Government of Canada for a Postdoctoral fellowship application No. 545643.

Appendix A: Details on the Derivation of the Energy Expressions

Taking the N th derivative w.r.t. λ (N = 0, 1, . . . 4) of Eq. (57), and using Eqs. (14a)-(14c), (15a), (15b) and (72), to evaluate the result at λ = 0 we obtain:

E (1) = 1 2 σ RTr v + u AR + H σσ(0) P σσ(1) + 1 2 σ RTr V σσ(1) -aK σσ(1) P σσ(0) + σ RTr u σσ SO P σσ(0) , (A1b)

+ 1 2 σσ RTr V σσ (1) -aK σσ (1) P σ σ(1) + 1 2 σ RTr C σσ(2) + V σσ(2) -aK σσ(2) P σσ(0) + σσ RTr u σσ SO P σ σ(1) , (A1c)

σ RTr v + uAR + H σσ(0) P σσ [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] + 1 2 σσ RTr V σσ (1) -aK σσ (1) P σ σ(2) + 1 2 σσ RTr C σσ (2) + V σσ (2) -aK σσ (2) P σ σ(1) + 1 2 σ RTr C σσ(3) + V σσ [START_REF] Almlöf | A variational approach to relativistic effects in LCAO calculations[END_REF] -aK σσ(3) P σσ(0)

and

Eqs. (A1d) and (A1e) make use of the third and fourth order perturbed density matrices:

and

In deriving Eqs. (A1a)-(A1e) we have used H σσ (0) = P σσ (0) = 0 for σ = σ, as given by Eq. (14a).

1. Zeroth Order Energy Substituting Eqs. ( 58) and (34) in Eq. (A1a), we obtain Eq. (59).

First Order Energy

Using Eqs. ( 12), ( 46) and (15b), we find that the terms appearing in the second and third lines Eq. (A1b) must vanish since R V σσ(1) = 0, R K σσ(1) = 0, which follows directly from Eq. ( 46) and I P σσ(0) = 0. Then, substituting Eq. ( 35) in the first line of Eq. (A1b), we find:

Since the trace is invariant to a cyclic permutation of the matrices in Eq. (A4) we obtain:

(A5) using Eqs. (45a)-(45b), as well as Eq. (58).

Finally, given that Θ σσ is real Hermitian and U σσ(1) OV is imaginary anti-Hermitian, both terms in Eq. (A5) can be combined to yield Eq. ( 60).

The Second Order Energy

From Eqs. (A1c) and (14b), we find:

We obtain the first line of Eq. (A6) by combining the first and third lines of Eq. (A1c) and noticing that, as an example, for the exchange operator, from Eq. (15b):

and for the xc operator:

and f σ commute, we obtain:

Next, we turn to the second term of Eq. (A17). With the aid of Eqs. ( 36) and ( 37), we obtain:

Then, taking advantage of the first order perturbation equation (Eq. ( 40)) in Eq. (A20) we find:

Substitution of Eqs. (48a) and (51) in Eq. (A21), leads to:

Finally, using Eqs. (A19) and (A22) in Eq. (A17), terms one through four of Eq. (A19) cancel with terms 1, 5, 4, 2, respectively, of Eq. (A22), such that we obtain Eq. (63).

Fourth Order Energy

From Eq. (A1e), we find the following simplified expression by using Eqs. (14a)-(14c), as well as Eqs. (A14)-(A16):

The first term of the first line of Eq. ( A23) is obtained by combining the first and fifth lines of Eq. (A1e). The second term of the first line of Eq. (A23) is obtained by combining the second, fourth and sixth lines of Eq. (A1e). Finally, the second line of Eq. ( A23) is obtained directly from the third line of Eq. (A1e) using Eq. (14c). We begin by considering the term H σσ(0) P σσ [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF] in Eq. (A23). Using the expression in Eq. (A3) for the fourth order density matrix, we find:

Given that f σ and σ commute, substitution of the fourth-order orthogonality condition (Eq. ( 32)) in Eq. (A24) yields:

Eq. (A25) may be written in terms of occ-occ occ-virt Eq. (A30), we obtain:

Let us now consider the term H σσ (2) P σ σ(2) in the last line of Eq. (A31). Writing the second order density matrix as in Eq. ( 36), we obtain:

Then, using Eq. (A32) in Eq. (A31), yields:

We can combine terms seven and eight of Eq. (A33) by writing:

The ninth and tenth terms of Eq. (A33) can be combined similarly:

Finally, using Eqs. (A34) and (A35), as well as the Hermiticity/anti-Hermiticity of the different blocks of G (1) , G (2) , U (1) and U (2) , Eq. (A33) becomes:

+ 2G

σσ (1)

Terms one through six of Eq. (A36) are obtained, respectively, by combining terms one plus four, two plus three, eleven plus thirteen, twelve plus fourteen, fifteen, and five plus six of Eq. (A33). Finally Eq. (A36) can be simplified algebraically to Eq. ( 64) by factorization.