Bastien Mussard

Sandeep Sharma

One-step treatment of spin-orbit coupling and electron correlation in large active spaces

In this work we demonstrate that the heat bath configuration interaction (HCI) and its semistochastic extension can be used to treat relativistic effects and electron correlation on an equal footing in large active spaces to calculate the low energy spectrum of several systems including halogens group atoms (F, Cl, Br, I), coinage atoms (Cu, Au) and the Neptunyl(VI) dioxide radical. This work demonstrates that despite a significant increase in the size of the Hilbert space due to spin symmetry breaking by the spin-orbit coupling terms, HCI retains the ability to discard large parts of the low importance Hilbert space to deliver converged absolute and relative energies. For instance, by using just over 10 7 determinants we get converged excitation energies for Au atom in an active space containing (150o,25e) which has over 10 30 determinants. We also investigate the accuracy of five different two-component relativistic Hamiltonians in which different levels of approximations are made in deriving the one-electron and two-electrons Hamiltonians, ranging from Breit-Pauli (BP) to various flavors of exact two-component (X2C) theory. The relative accuracy of the different Hamiltonians are compared on systems that range in atomic number from first row atoms to actinides.

I. INTRODUCTION

Relativistic effects play an important role in a variety of photochemical processes, various magnetic spectroscopies and are responsible for unusual phase behaviors in transition metal oxides. In practice, the most accurate procedure for treating relativistic effects in molecular systems is through the solution of the Dirac-Coulomb-Breit (DCB) equation [START_REF] Barysz | Relativistic Methods for Chemists[END_REF][START_REF] Liu | Ideas of relativistic quantum chemistry[END_REF][START_REF] Fleig | Invited review: Relativistic wave-function based electron correlation methods[END_REF][START_REF] Saue | Relativistic Hamiltonians for Chemistry: A Primer[END_REF][START_REF] Kutzelnigg | Solved and unsolved problems in relativistic quantum chemistry[END_REF][START_REF] Pyykkö | Relativistic Effects in Chemistry: More Common Than You Thought[END_REF][START_REF] Pyykkö | The Physics behind Chemistry and the Periodic Table[END_REF] . Although recent developments have made it possible to perform self-consistent field [START_REF] Pernpointner | Parallelization of four-component calculations. I. Integral generation, SCF, and four-index transformation in the DiracFock package MOLFDIR[END_REF][START_REF] Yanai | A new computational scheme for the DiracHartreeFock method employing an efficient integral algorithm[END_REF][START_REF] Thyssen | A direct relativistic four-component multiconfiguration selfconsistent-field method for molecules[END_REF][START_REF] Kelley | Large-scale Dirac Fock Breit method using density fitting and 2-spinor basis functions Large-scale Dirac Fock Breit method using density fitting and 2-spinor basis functions[END_REF] , density functional theory [START_REF] Yanai | A new implementation of four-component relativistic density functional method for heavy-atom polyatomic systems[END_REF][START_REF] Saue | Four-component relativistic KohnSham theory[END_REF][START_REF] Belpassi | Recent advances and perspectives in fourcomponent Dirac Kohn Sham calculations[END_REF][START_REF] Mizukami | A dual-level approach to four-component relativistic density-functional theory[END_REF] , coupled cluster [START_REF] Nataraj | General implementation of the relativistic coupled-cluster method General implementation of the relativistic coupled-cluster method[END_REF][START_REF] Sørensen | Two-and four-component relativistic generalized-active-space coupled cluster method : Implementation and application to BiH Two-and four-component relativistic generalizedactive-space coupled cluster method : Implementation and application to BiH[END_REF] , explicitly correlated [START_REF] Bylicki | Relativistic Hylleraas configuration-interaction method projected into positive-energy space[END_REF][START_REF] Li | Relativistic explicit correlation : Coalescence conditions and practical suggestions Relativistic explicit correlation : Coalescence conditions[END_REF][START_REF] Ten-No | Communication : Explicitly correlated four-component relativistic secondorder Møller-Plesset perturbation theory Communication : Explicitly correlated four-component relativistic secondorder Møller-Plesset perturbation theory[END_REF] and multireference [START_REF] Almoukhalalati | Electron correlation within the relativistic no-pair approximation Electron correlation within the relativistic no-pair approximation[END_REF][START_REF] Bates | Fully relativistic complete active space self-consistent field for large molecules : Fully relativistic complete active space selfconsistent field for large molecules : Quasi-second-order minimax optimization[END_REF][START_REF] Fleig | The generalized active space concept for the relativistic treatment of electron correlation . II . Large-scale configuration interaction implementation based on relativistic 2-and 4-spinors and its application The generalized active space concept for the[END_REF][START_REF] Fleig | The generalized active space concept for the relativistic treatment of electron correlation . III . Large-scale configuration interaction and multiconfiguration self-consistent-field four-component methods with application to The generalized active space[END_REF][START_REF] Knecht | group general active space implementation with application to BiH Large-scale parallel configuration interaction . II . Two-and four-component double-group general active space implementation with application to BiH[END_REF][START_REF] Shiozaki | Relativistic Internally Contracted Multireference Electron Correlation Methods[END_REF] calculations on the DCB Hamiltonian, these methods remain problematic due to their high computational cost.

An alternative is to use two-component Hamiltonians which are capable of delivering quantitative accuracy for relativistic problems. They are derived by either eliminating the small component (ESC) or by using a unitary transformation called the Foldy-Wouthuysen (FW) transformation [START_REF] Foldy | On the Dirac Theory of Spin 1/2 Particles and Its Non-Relativistic Limit[END_REF] to decouple the large and the small components in the one-body Dirac or Dirac-Fock Hamiltonian. These two are of course related and FW transformation can be viewed as a normalized elimination of small component (NESC). This transformation is also applied to the two-electron integrals to generate a relativistically correct electron interaction. The resulting two-component Hamiltonian can be partitioned into the Schrödinger equation and additional relativistic spin-free and spin-dependent terms. The Schrödinger equation together with the spin-free terms are referred to as the scalar relativistic Hamiltonian and calculating the ground state energy of this Hamiltonian is usually no more expensive than that of the Schrödinger equation. The remaining spin-dependent terms, known as the spinorbit coupling (SOC), are usually treated using perturbation theory. However, for heavier elements these terms can become large and it is more appropriate to treat them on an equal footing with the spin-free terms and electron correlation. In this work we follow the latter approach.

The three most commonly used two-component Hamiltonians are the Douglass-Kroll-Hess [START_REF] Douglas | Quantum electrodynamical corrections to the fine structure of helium[END_REF][START_REF] Hess | Relativistic electronic-structure calculations employing a two-component no-pair formalism with external-field projection operators[END_REF][START_REF] Nakajima | The DouglasKrollHess Approach[END_REF][START_REF] Reiher | DouglasKrollHess Theory: a relativistic electrons-only theory for chemistry[END_REF] , the Barysz-Sadlej-Snijders [START_REF] Barysz | Infinite-order twocomponent theory for relativistic quantum chemistry[END_REF][START_REF] Barysz | Recovering four-component solutions by the inverse transformation of the infinite-order two-component wave functions[END_REF][START_REF] Iliaš | Theoretical study of PbO and the PbO anion[END_REF] and the exact two-component (X2C) Hamiltonians [START_REF] Dyall | Interfacing relativistic and nonrelativistic methods. I. Normalized elimination of the small component in the modified Dirac equation[END_REF][START_REF] Kutzelnigg | Quasirelativistic theory equivalent to fully relativistic theory[END_REF][START_REF] Liu | Quasirelativistic theory. II. Theory at matrix level[END_REF][START_REF] Iliaš | An infinite-order twocomponent relativistic Hamiltonian by a simple one-step transformation[END_REF] .

Both Douglass-Kroll-Hess and Barysz-Sadlej-Snijders Hamiltonians carry out the transformation analytically. The transformation operator contains a complicated function of the momentum operator and its integrals cannot be calculated analytically, instead the matrix representation of the momentum operator in a finite basis is calculated. In X2C, one forgoes the analytic transformation and the entire process is carried out algebraically using the matrix representation of operators, starting from the solution of the Dirac equation. The commonly used approximation, known as X2C-1e, consists of solving the spin-free non-interacting Dirac equation in one step (without electron-electron interaction a self-consistent procedure is not needed) and the transformation matrix is derived from this solution. The motivation for doing so comes from the fact that only a small fraction of the total cost of performing a correlated calculation is used in the solution of the onebody problem and thus solving it to exactly decouple the large and small component is relatively cheap.

To get quantitative accuracy one also needs to include the relativistically correct two-body terms. Such terms are derived from the transformation of the Coulomb and Gaunt operators, which gives rise to several spin-free and spin-dependent terms. In this work we will ignore the relativistic correction to the spin-free two-body terms and only some of the spin-dependent terms are included that appear at order α 2 . The effect of those two-body SOC terms is often treated approximately using the spin-orbit mean field (SOMF) approximation [START_REF] Heß | A mean-field spin-orbit method applicable to correlated wavefunctions[END_REF][START_REF] Marian | A new meanfield and ECP-based spin-orbit method. Applications to Pt and PtH[END_REF] , which replaces the two-body terms by an effective fock-like one-body term. Although the approximation sounds drastic, in practice it is known to be extremely accurate and is almost universally used. The SOMF approximation sim-plifies the correlated calculation considerably by reducing the memory and CPU requirement of storing and transforming the different sets of relativistic two-electron integrals.

Several algorithms and programs are now available for performing relativistic calculations with two-component Hamiltonians [START_REF] Buenker | Comparison of spin-orbit configuration interaction methods employing relativistic effective core potentials for the calculation of zero-field splittings of heavy atoms with a 2Po ground state[END_REF][START_REF] Berning | Spin-orbit matrix elements for internally contracted multireference configuration interaction wavefunctions[END_REF][START_REF] Vallet | A two-step uncontracted determinantal effective Hamiltonian-based SOCI method[END_REF][44][START_REF] Gordon | In Theory and Applications of Computational Chemistry[END_REF][START_REF] Kleinschmidt | CI: A multireference spin-orbit configuration interaction method for large molecules[END_REF][START_REF] Gauss | Calculation of Electronic g -Tensors using Coupled Cluster Theory[END_REF][START_REF] Neese | The ORCA program system[END_REF][START_REF] Mai | Perturbational treatment of spin-orbit coupling for generally applicable highlevel multi-reference methods[END_REF][START_REF] Werner | version 2015.1, a package of ab initio programs[END_REF][START_REF] Aquilante | Molcas 8: New capabilities for multiconfigurational quantum chemical calculations across the periodic table[END_REF] . Some of these methods include the spin-orbit coupling in the self-consistent field (SCF) calculation, due to which the orbital relaxation effects are fully included, but the resulting orbitals are complexvalued spinors [START_REF] Esser | Complex two-and four-index transformation over twocomponent kramers-degenerate spinors[END_REF][START_REF] Kim | Spin-orbit effects calculated by a configuration interaction method using determinants of two-component molecular spinors: test calculations on Rn and T1H[END_REF][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] Kim | Two-component multiconfigurational second-order perturbation theory with Kramers restricted complete active space self-consistent field reference function and spin-orbit relativistic effective core potential[END_REF][START_REF] Kim | Spin-orbit coupling and electron correlation in relativistic configuration interaction and coupled-cluster methods[END_REF] . In the intermediate approach the spin-orbit coupling terms are only introduced during the correlated calculation [START_REF] Kleinschmidt | CI: A multireference spin-orbit configuration interaction method for large molecules[END_REF][START_REF] Sjøvoll | A determinantal approach to spin-orbit configuration interaction[END_REF][START_REF] Yabushita | SpinOrbit Configuration Interaction Using the Graphical Unitary Group Approach and Relativistic Core Potential and SpinOrbit Operators[END_REF][START_REF] Fedorov | A study of the relative importance of one and two-electron contributions to spinorbit coupling[END_REF][START_REF] Ganyushin | A fully variational spin-orbit coupled complete active space self-consistent field approach: application to electron paramagnetic resonance g-tensors[END_REF] , and are treated on an equal footing with the dynamical electron correlation. A further approximation is possible which is usually called the quasi-degenerate perturbation theory or the state interaction approach, whereby one first calculates several electronic state wavefunctions of the spin-free Hamiltonian in different spin sectors. The matrix representation in the basis of these states of the complete two-component Hamiltonian including spin-orbit coupling is then diagonalized to obtain the spin-orbit coupled results [START_REF] Roos | Relativistic quantum chemistry: the multiconfigurational approach[END_REF][START_REF] Ganyushin | First-principles calculations of zero-field splitting parameters[END_REF][START_REF] Roemelt | Spin orbit coupling for molecular ab initio density matrix renormalization group calculations: Application to g-tensors[END_REF][START_REF] Sayfutyarova | A state interaction spin-orbit coupling density matrix renormalization group method[END_REF][START_REF] Knecht | A Nonorthogonal State-Interaction Approach for Matrix Product State Wave Functions[END_REF] . The accuracy of the method is limited by the number of states included in the state interaction approach, which is usually on the order of a few tens of states and in rare cases of heavy elements a few hundred [START_REF] Liu | A Stable Pentagonal Bipyramidal Dy(III) Single-Ion Magnet with a Record Magnetization Reversal Barrier over 1000 K[END_REF] . A related approach is the EOM(SOC)-CC [START_REF] Klein | Perturbative calculation of spin-orbit splittings using the equation-of-motion ionization-potential coupled-cluster ansatz[END_REF][START_REF] Cao | Combining the spin-separated exact two-component relativistic Hamiltonian with the equation-of-motion coupled-cluster method for the treatment of spin-orbit splittings of light and heavy elements[END_REF] where the spin-orbit coupling terms are included only during the equation of motion part of the calculation and thus are effectively treated perturbatively.

In this work we perform multireference calculations in large active spaces where the SOC is treated on an equal footing with the electron correlation. Heat-bath configuration interaction (HCI) algorithm [START_REF] Holmes | Heat-bath Configuration Interaction: An efficient selected CI algorithm inspired by heat-bath sampling[END_REF][START_REF] Sharma | Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory[END_REF] is the multireference method which is used to calculate the zero-field splittings using real-valued orbitals obtained from a nonrelativistic SCF calculation. HCI is an efficient variant of the general class of methods in which a selected configuration interaction is performed which is followed by perturbation theory. As we will show in the results section, HCI is able to systematically discard large parts of the Hilbert space to deliver accurate excitation energies at a much reduced cost relative to a full configuration interaction (FCI) calculation. To the best of our knowledge current calculations and the density matrix renormalization group (DMRG) theory calculations performed by Sayfutyarova et al. [START_REF] Sayfutyarova | A state interaction spin-orbit coupling density matrix renormalization group method[END_REF] and Knecht et al. [START_REF] Knecht | Communication: Four-component density matrix renormalization group[END_REF] are the only ones in literature where relativistic calculations are performed on an active space larger than 16 electrons in 16 orbitals. However, unlike the DMRG calculations of Sayfutyarova, here we treat SOC non-perturbatively on an equal footing with the electron correlation.

The rest of the article is organized as follows. In Section II, we present the derivation and the working equations of the Breit-Pauli and X2C Hamiltonians used in this work. In Section III, we describe the Heat-bath Configuration Interaction algorithm and the extensions to it that allow us to treat the SOC terms. In Section IV and V, we give the computational details and results respectively.

II. THEORY

In this section, we present the working equations of the two-component Hamiltonians used in this work: the Breit-Pauli (BP) and the exact two-component (X2C) Hamiltonians. To give some context, we start with the matrix formulation of the one-electron four-component Dirac Hamiltonian

V ¯ne c σ σ σ • p c σ σ σ • p ¯V ¯ne -2mc 2 Ψ = EΨ, (1)
where V ¯ne is the electron-nucleus potential, p ¯is the momentum and σ σ σ = {σ σ σ x , σ σ σ y , σ σ σ z } is the set of Pauli matrices. The eigenvectors Ψ of the Dirac Hamiltonian are four-component bispinors which contain the "large" Ψ L and "small" Ψ S components that are themselves twocomponent wavefunctions:

Ψ = Ψ L Ψ S . (2)
Although the Dirac Hamiltonian is fully Lorentz invariant, introducing the Coulomb interaction between the electrons breaks this invariance and one would need to go to quantum electrodynamics (QED) to obtain a fully Lorentz invariant theory of interacting electrons, but it is not obvious that a fully Lorentz invariant many-electron Hamiltonian can be derived from QED. In practice, the Dirac-Coulomb-Breit Hamiltonian with the non-retarded electron-electron interaction

Vee (i, j) = 1 r ij - α α α i • α α α j r ij + α α α i • α α α j 2r ij - (α α α i • r ¯ij )(α α α j • r ¯ij ) 2r 3 ij (3)
is used (in brackets are the Coulomb, Gaunt and Gauge terms, while α α α i are the Dirac matrices) and is expected to be sufficiently accurate for chemical applications.

In principle, it is possible to obtain exact twocomponent Hamiltonians by solving the DCB Hamiltonian and using the many-body solutions to decouple the small and the large components exactly. This is of course impractical because the full solution to the DCB Hamiltonian would be required. Instead, one tries to decouple the large and small components at the one-body level, which can be done approximately analytically (as is done to derive the Breit-Pauli Hamiltonian) or exactly algebraically (as is done to derive the exact two-component Hamiltonian). This class of procedures is tractable because one needs to solve at most the one-body problem. The transformation obtained this way is then applied to the two-body Coulomb, Gaunt and Breit terms to derive two-component two-body SOC Hamiltonians.

A. Breit-Pauli Hamiltonian

The Breit-Pauli (BP) Hamiltonian used in this work is a two-component Hamiltonian obtained through an analytic FW transformation performed on the fourcomponent one-body Dirac Hamiltonian of Eq. ( 1) with two-body terms obtained by contributions from the transformed Coulomb and Gaunt terms (first and second terms of Eq. ( 3)). The final expression of the Breit-Pauli Hamiltonian [START_REF] Dyall | Introduction to Relativistic Quantum[END_REF][START_REF] Reiher | Relativistic Quantum Chemistry[END_REF] is

H ¯BP = H ¯BP SF + H ¯BP SO (1) + H ¯BP MF (1) , (4)
where the three terms are defined in Eqs. ( 5), ( 6) and (15).

One-body part

The one-body part of the BP Hamiltonian, known as the Pauli Hamiltonian, contains the kinetic energy, electron-nuclear interaction, one-body spin-orbit coupling, mass-velocity and Darwin terms. In this work we will ignore the divergent mass-velocity and Darwin terms and are thus left with H ¯Pauli = H ¯BP SF + H ¯BP SO [START_REF] Barysz | Relativistic Methods for Chemists[END_REF] , where the spin-free part

H ¯BP SF = T ¯+ V ¯ne (5)
contains the kinetic energy T ¯and the electron nuclear attraction V ¯ne . The one-body spin-orbit coupling term is

H ¯BP SO (1) = -i α 2 4 iA p ¯i Z A r iA × p ¯i • σ σ σ i , (6)
where, α is the fine structure constant, p ¯i is the momentum of electron i, Z a is the atomic number of nucleus A, and r iA the distance between electron i and nucleus A.

In second quantization the one-body SOC term can be written as

H ¯BP SO (1) = -i α 2 4 iA pq µνλ λ µν p µ | Z A r iA |q ν σλ pq (7) = -i α 2 4 pqλ BP S λ pq σλ pq , (8)
where λ, µ, ν can be (x, y, z), p and q are orbitals, p µ = ∂p ∂µ are nuclear derivatives of the orbitals, and λ µν is the Levi-Civita symbol. Note that the three σ operators are

σx pq = â † pα âqβ + â † pβ âqα σy pq = -i(â † pα âqβ -â † pβ âqα ) (9)
σz pq = â † pα âqα -â † pβ âqβ .
We can infer from the definition of BP S λ pq that the onebody spin-orbit coupling integrals are anti-hermitian.

Two-body part

The application of the FW transformation to the Coulomb and Gaunt operators spawns several spin-free and spin-dependent terms, out of which in this work we will only include the spin-same-orbit and spin-other-orbit interactions, i.e. respectively the first and second term in

H ¯BP SO (2) =i α 2 4 i =j p ¯i 1 r ij × p ¯i • (σ σ σ i + 2σ σ σ j ). ( 10
)
In second quantization, this takes the form

H ¯BP SO (2) =i α 2 4 pqrsλ 2J λ rspq + J λ pqrs Dλ pqrs , (11)
where the two-electron integrals are

J λ pqrs = µν λ µν p µ r|q ν s , (12)
and the operator is

Dλ pqrs = σλ pq Êrs -δ qr σλ ps . (13)
For real orbitals one can utilize the 4-fold symmetry

J λ pqrs = -J λ qprs = -J λ qpsr = J λ pqsr (14)
to reduce the memory cost of storing the two electron integrals.

Note that the two-body spin-orbit coupling term is of opposite sign to the one-body term of Eq. ( 8) and acts as a screening potential. Several other terms including the spin-spin dipole interactions and Fermi contact interaction are ignored. For most applications such simplifications are justified, however some noteworthy exceptions exist such as the oxygen molecule [START_REF] Langhoff | Ab initio evaluation of the fine structure of the oxygen molecule[END_REF][START_REF] Vahtras | Ab initio calculations of zero-field splitting parameters[END_REF] , where the spinspin dipole interactions are of the same order of magnitude as the spin-orbit coupling terms.

Although the two-body terms can be easily incorporated in a calculation, the cost of storing and transforming integrals for large molecules become too expensive. This motivates the use of the spin-orbit mean field approximation [START_REF] Heß | A mean-field spin-orbit method applicable to correlated wavefunctions[END_REF][START_REF] Marian | A new meanfield and ECP-based spin-orbit method. Applications to Pt and PtH[END_REF][START_REF] Li | On the spin separation of algebraic two-component relativistic Hamiltonians: Molecular properties[END_REF] , where one writes an effective onebody operator by integrating out the spins of two out of the four orbitals in the two-electron integrals of Eq. (12). Such an integration leads to the effective one-body term [START_REF] Heß | A mean-field spin-orbit method applicable to correlated wavefunctions[END_REF] written in second quantization as

H ¯BP MF (1) = i α 2 4 pq BP Sλ pq σλ ps , (15)
where

BP Sλ pq = rs γ rs J λ pqrs - 3 2 J λ prsq - 3 2 J λ sqpr ( 16
)
and γ is the one-body density matrix. Although it is not used here, the cost of calculating the two-body terms can be further substantially reduced by using the one-center approximation which makes use of the local nature of the spin-orbit interaction and introduces only a small error 77 .

B. X2C Hamiltonian

The derivation of the X2C Hamiltonian begins by using the restricted kinetic balance 78 that replaces the small component (ψ S ) with the pseudo-large component (φ L )

ψ S = σ σ σ • p 2c φ L , (17)
which gives rise to the modified four-component oneelectron Dirac equation

V ¯ne T T ¯W ¯-T ¯ ψ L φ L = S ¯0 0 ¯α2 4 T ¯ ψ L φ L E, (18)
where

W ¯= α 2 4 (σ σ σ •p ¯)V ne (σ σ σ •p ¯)
and S ¯is the non-relativistic metric. The restricted kinetic balance transformation ensures that the large and the pseudo-large components have the same symmetry properties and can be expanded with a common basis set. Also, the two components become equal to each other in the non-relativistic limit.

The only spin-dependent term in the modified Dirac equation of Eq. ( 18) is found in W ¯, which can be exactly partitioned into a spin-free and spin-dependent term

W ¯= α 2 4 p ¯• V ne p ¯+ i α 2 4 σ σ σ • p ¯Vne × p ¯ =W ¯SF + W ¯SD . ( 19
)
This allows us to separate the spin-free (H ¯4c SF ) and spindependent (H ¯4c SD ) terms in the four-component modified Dirac equation

H ¯4c SF = V ¯ne T T ¯W ¯SF -T ¯ (20) H ¯4c SD = 0 ¯0 0 ¯W ¯SD . ( 21
)
An FW transformation can be derived by solving only the spin-free part of the one-body Hamiltonian (Eq. 20) to obtain the X2C-1 one-body SOC terms [START_REF] Li | On the spin separation of algebraic two-component relativistic Hamiltonians: Molecular properties[END_REF][START_REF] Li | On the spin separation of algebraic two-component relativistic Hamiltonians[END_REF][START_REF] Li | Combining spin-adapted open-shell TD-DFT with spinorbit coupling[END_REF] . Another possibility is to obtain a complex-valued FW transformation by solving the total one-body Hamiltonian (containing the total W ¯of Eq. 19) to obtain the X2C-N one-body SOC terms [START_REF] Cheng | Perturbative treatment of spin-orbit coupling within spin-free exact two-component theory[END_REF] . In both cases the two-body part is obtained by applying the X2C-1 FW transformation to the Coulomb and Gaunt operators. The final expression of the resulting exact two-component Hamiltonians used in this work is thus:

H ¯X2C-x = H ¯X2C SF + H ¯X2C-x SO (1) + H ¯X2C MF (1) , (22)
where "x" distingishes between the X2C-1 and X2C-N one-body terms. The H ¯X2C SF is defined in Eq. ( 28), H ¯X2C-1 SO [START_REF] Barysz | Relativistic Methods for Chemists[END_REF] is defined in Eq. ( 29), H ¯X2C-N SO [START_REF] Barysz | Relativistic Methods for Chemists[END_REF] is defined in Eq. (31) and the H ¯X2C MF [START_REF] Barysz | Relativistic Methods for Chemists[END_REF] is defined in Eq. (32).

One-body part: X2C-1 scheme

As stated above, in the X2C-1 scheme one uses only the spin-free four-component Hamiltonian H ¯4c

SF to solve the modified Dirac equation of Eq. ( 18) in a single step. The data of the obtained four-component positive energy solutions Ψ + is used to apply an FW transformation to block-diagonalize H ¯4c SF and retain its large component block. The same FW transformation is subsequently used on H ¯4c SD and on the Coulomb and Gaunt terms to obtain the complete X2C-1 Hamiltonian.

Without going into the detailed derivation which can be found for example in Ref. 82, it can be shown that the following FW transformation matrix block-diagonalizes the spin-free modified Dirac equation

U ¯= 1 ¯X X ¯1 ¯ R ¯0 0 ¯R ¯ , ( 23
)
where the first matrix does the decoupling and the second renormalizes the states to bring them to the Schrödinger metric. Only X ¯and R ¯are required to determine the large component of the block-diagonalized modified Dirac Hamiltonian. The matrix X ¯is defined through the relation

φ L + = X ¯ψL + ( 24
)
between the large component ψ L + and the pseudo large component φ L + of the positive energy bispinors

Ψ + = ψ L + φ L + ( 25
)
and the matrix R ¯is defined as

R ¯= (S ¯-1 S ¯)-1/2 (26) S ¯= S ¯+ α 2 2 X ¯ †T ¯X ¯. ( 27
)
Using this unitary transformation, one can show that the spin-free two-component X2C Hamiltonian H ¯X2C SF is given as

H ¯X2C SF = R ¯ †(V ¯ne + T ¯X ¯+ X ¯ †T ¯+ X ¯ †(W ¯SF -T ¯)X ¯)R ¯. ( 28
)
The same FW transformation is subsequently applied to the spin-dependent four-component Hamiltonian H ¯4c SD to obtain the one-body spin-orbit coupling operator H ¯X2C-1

SO (1) = -i α 2 4 Ai (X ¯R ¯) † p ¯i Z ai r Ai × p ¯i • σ σ σ i (X ¯R ¯), (29)
which involves integrals

X2C S λ pq = rs (X ¯R ¯) † pr BP S λ pq (X ¯R ¯)sq , (30)
that are equal to the integrals found in the Breit-Pauli case transformed by the regularizing matrices X ¯and R ¯.

2. One-body part: X2C-N scheme

An alternative scheme is to solve the complete modified Dirac equation found in Eq. ( 18), including both the spin-free and spin-dependent parts, to calculate the corresponding complex-valued FW transformation. The expressions for calculating the matrices X ¯and R ¯are formally equivalent to Eqs. ( 24) and ( 26) respectively, with the difference that these are now 2n × 2n dimensional matrices. The X ¯and R ¯matrices can be used in Eqs. ( 28) and ( 29) to obtain the complete two-component one-body hamiltonian containing both the spin-free and spin-dependent terms

H ¯X2C-N = H ¯X2C-N SF + H ¯X2C-N SO (1) . ( 31
)
It is worth pointing out that the H ¯X2C-N SF contains contributions from high powers of W ¯SD and is not the same as the genuine spin-free Hamiltonian H ¯X2C SF of Eq. ( 28). However, as expression in Eq. ( 22) indicates, we nevertheless use H ¯X2C

SF in the X2C-N scheme.

Two-body part

The two-body operator and its mean field approximation for the X2C theory are derived in great detail in Ref. 79 and 76 and we merely give a brief outline of the derivation.

The X2C two-body operators are obtained by performing the same two sequential transformations (the restricted kinetic balance transformation and the spinfree unitary transformation) described in the "One-body part: X2C-1 scheme" section. As a result, one obtains several terms, the first of which is the modified spin-free Coulomb term. This term is customarily replaced by the bare Coulomb term and the difference is ignored. Out of the remaining spin-dependent terms only the spin-sameorbital and spin-other-orbital are retained in our work. These are similar to their Breit-Pauli counterparts multiplied by the regularizing matrices X ¯and R ¯. The mean field approximation to the two-body part is derived by performing the X2C transformation on the spin-orbit contribution of the two-electron integrals in the modified Dirac-Fock-Coulomb-Breit operator [START_REF] Li | On the spin separation of algebraic two-component relativistic Hamiltonians: Molecular properties[END_REF] . The resulting X2C mean field operator is

H ¯X2C MF (1) = i α 2 4 pq X2C Sλ pq σλ ps , (32)
where the one-body integrals are

X2C Sλ pq = R pr g ¯LL,λ + g ¯LS,λ + g ¯SL,λ + g ¯SS,λ rs R sq . (33)
The matrices g ¯are g LL,λ mn = -

pq 2K λ pmqn γ SS pq g LS,λ mn = - pqr K λ mpqr + K λ pmqr γ LS pq X rn = g SL,λ nm (34)
g SS,λ mn = -

pqrs 2 K λ rspq + K λ rsqp -K λ rpsq γ LL pq X † mp X qn ,
the matrices X ¯and R ¯are formally defined in Eq. ( 24) and ( 26) and the matrices γ γ γ are

γ γ γ LL = R ¯ †γ γ γ SA R ¯(35) γ γ γ LS = γ γ γ LL X ¯= γ γ γ SL † (36) γ γ γ SS = X ¯ †γ γ γX ¯. ( 37
)
where γ γ γ SA = 1 2 (γ γ γ α + γ γ γ β ) is the spin-averaged density matrix. Finally, the two-electron integrals appearing in Eq. ( 34) are defined as

K λ pqrs = µν λ µν p µ r ν |qs . (38)
The mean field X2C Hamiltonian exactly reduces to the mean field BP Hamiltonian when R ¯= X ¯= I ¯.

III. HEAT-BATH CONFIGURATION INTERACTION

Heat-bath Configuration Interaction (HCI) is a recently-developed 69 variant of the general class of methods that perform a selected configuration interaction calculation followed by perturbation theory (SCI-PT) [START_REF] Liu | A Stable Pentagonal Bipyramidal Dy(III) Single-Ion Magnet with a Record Magnetization Reversal Barrier over 1000 K[END_REF][START_REF] Ivanic | Identification of deadwood in configuration spaces through general direct configuration interaction[END_REF][START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF][START_REF] Buenker | Individualized configuration selection in CI calculations with subsequent energy extrapolation[END_REF][START_REF] Evangelisti | Convergence of an improved CIPSI algorithm[END_REF][START_REF] Harrison | Approximating full configuration interaction with selected configuration interaction and perturbation theory[END_REF][START_REF] Steiner | The efficient treatment of higher excitations in CI calculations. A comparison of exact and approximate results[END_REF][START_REF] Wenzel | Multireference Basis-Set Reduction[END_REF][START_REF] Neese | A spectroscopy oriented configuration interaction procedure[END_REF][START_REF] Abrams | Important configurations in configuration interaction and coupledcluster wave functions[END_REF][START_REF] Bytautas | A priori identification of configurational deadwood[END_REF][START_REF] Knowles | Compressive sampling in configuration interaction wavefunctions[END_REF][START_REF] Schriber | Communication: An adaptive configuration interaction approach for strongly correlated electrons with tunable accuracy[END_REF][START_REF] Caffarel | Using CIPSI Nodes in Diffusion Monte Carlo[END_REF][START_REF] Garniron | Hybrid stochastic-deterministic calculation of the second-order perturbative contribution of multireference perturbation theory[END_REF] . HCI, similar to other SCI-PT methods, consists of a variational step and a perturbative step. In the variational step, a set of important determinants is iteratively identified and the Hamiltonian is diagonalized in the space of these determinants. In the perturbative step, the energy obtained during the variational step is corrected using Epstein-Nesbet perturbation theory. For details of the algorithm we refer the reader to our recent publications [START_REF] Holmes | Heat-bath Configuration Interaction: An efficient selected CI algorithm inspired by heat-bath sampling[END_REF][START_REF] Sharma | Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory[END_REF][START_REF] Holmes | Excited states using semistochastic heat-bath configuration interaction[END_REF][START_REF] Smith | Cheap and near exact CASSCF with large active spaces[END_REF] , here we briefly describe the key aspects of the method relative to other SCI-PT approaches.

A. The heat-bath criterion

At each iteration of the variational stage, the multireference ground state wavefunction

|Ψ = Di∈V c i |D i (39)
with energy E 0 is calculated by diagonalizing the Hamiltonian in the current space V of important determinants |D i to obtain their coefficient c i . The space V is augmented with a set of new determinants |D a that satisfy the HCI criterion max

Di∈V |H ai c i | > 1 , (40)
where H ai = D a | Ĥ|D i and 1 is a user-defined parameter. The HCI criterion is different from the one used in CIPSI [START_REF] Huron | Iterative perturbation calculations of ground and excited state energies from multiconfigurational zeroth-order wavefunctions[END_REF][START_REF] Evangelisti | Convergence of an improved CIPSI algorithm[END_REF] , which is based on the contribution of a determinant D a to the perturbative correction to the wavefunction

|Di ∈V H ai c i E 0 -E a > 1 . (41)
Although the HCI criterion is in principle suboptimal at picking out the most important determinants, in practice the difference is minimal. Moreover, this slight difference is more than made up for by the significant advantage of the HCI criterion, which is that one generates only the determinants |D a that satisfy the criterion in Eq. ( 40), and no resources are spent on generating determinants that would have been discarded. The same criterion can also be used to speed up perturbation theory step.

B. Stochastic perturbation theory

The Epstein-Nesbet perturbative correction is evaluated as

E 2 ( 2 ) = |Da 1 E 0 -E a   ( 2) |Di ∈V H ai c i   2 , ( 42
)
where the symbol ( 2) designates a "screened sum" in which terms smaller in magnitude than a user-defined parameter 2 are discarded and where the first sum is over the determinants |D a that meet this criterion. The exact perturbation correction is recovered in the limit of 2 → 0, and an 2 much smaller than 1 is needed in order to obtain a good approximation to the perturbative correction.

Although the screened sum allows one to discard a large fraction of determinants, in most cases this is not sufficient to eliminate the memory bottleneck of having to store all the determinants |D a and their perturbative contributions in memory. To overcome this memory bottleneck, we use semistochastic perturbation theory [START_REF] Sharma | Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory[END_REF] to estimate the perturbative correction. In the semistochastic perturbation theory, an initial deterministic calculation is performed with a relatively loose parameter d 2 to obtain an approximate perturbative correction E D 2 ( d 2 ). The error in this calculated energy is then corrected stochastically by performing several iterations in which a stochastic perturbative correction is evaluated using the loose parameter d 2 and a much tighter parameter 2 . The near-exact perturbative correction with the tight parameter 2 can then be calculated as

E 2 ( 2 ) = E D 2 ( d 2 ) + E S 2 ( 2 ) -E S 2 ( d 2 ) . (43)
The key to the success of the semistochastic perturbation theory is that both the loose and tight stochastic perturbative corrections, E S 2 ( d 2 ) and E S 2 ( 2 ), are calculated with the exact same set of sampled determinants. This correlated sampling significantly reduces the variance in the estimated perturbative corrections, thereby reducing the stochastic error.

C. Spin-Orbit coupling and excited states

With the introduction of the spin-orbit coupling, the Hamiltonian does not commute with the Ŝz operator and consequently Ŝz is no longer a good quantum number. The SOC terms introduce non-zero matrix elements between determinants that contain different numbers of α and β electrons, and thus break the degeneracy between the 2S + 1 states of a spin multiplet. The energy splitting can be measured experimentally to determine the zero-field splitting (ZFS), and these energy differences are usually a small fraction of the absolute energies of the molecule. Thus two complications need to be addressed with the HCI algorithm: forming the wavefunction consisting of determinants with different numbers of α and β electrons, and the special care needed to calculate the small energy differences with sufficient accuracy.

To address the first complication, let us recall that the variational step of the HCI algorithm includes three operations, identifying the important determinants to add to the space V, finding the non-zero Hamiltonian matrix elements between all the determinants of V and diagonalizing the Hamiltonian. Each of these three operations are modified due to the addition of the SOC terms. In identifying the important determinants, one now also has to include those that can be generated with an α → β or a β → α excitation. On the other hand, the advantageous search protocol for doubly excited determinants that meet the HCI criterion does not change because the SOC terms only contain one-body operators. The Hamiltonian is stored in the sparse storage format and most of the terms in it are constructed using the same algorithm as before [START_REF] Sharma | Semistochastic Heat-bath Configuration Interaction method: selected configuration interaction with semistochastic perturbation theory[END_REF] , which avoids having to perform the expensive double loop over all determinants in the space V. Additionally, matrix elements that break Ŝz symmetry are identified for each determinant by looping over all possible α-occupied β-unoccupied and β-occupied αunoccupied orbital pairs to generate new determinants and doing a binary search over V to check if the determinant is present. Finally, the generalization of the Davidson algorithm to diagonalize the complex-valued Hamiltonian poses no serious problem.

Calculating the small energy splitting accurately requires that the energy of each state is calculated with a similar accuracy. The standard procedure for doing so is the state-average algorithm and was recently introduced in the context of HCI by modifying the HCI criterion of Eq. ( 40) [START_REF] Holmes | Excited states using semistochastic heat-bath configuration interaction[END_REF] . In this work we introduce yet another different criterion that is more suitable for targeting degenerate states. The updated HCI criterion is

max Di∈V |H ai ci | > 1 ci = n |c n i | 2 ( 44
)
where ci is the square root of the sum of the squares of the coefficients c n i of determinant |D i in the n th state. This ensures that unitary transformations between degenerate states do not change the determinants that are added to V.

IV. COMPUTATIONAL DETAILS

Several different combinations of one-body scalar relativistic, two-body scalar relativistic, one-body SOC and two-body mean field SOC terms can be used. Here we perform calculations with 5 different Hamiltonians where the two-body scalar relativistic term is the Coulomb interaction while the other three terms are specified below:

• "bp-bp" uses the scalar relativistic, one-body and mean field two-body parts of the Breit-Pauli Hamiltonian. This involves the integrals in Eqs. ( 8) and ( 16).

• "x2cn-bp" uses the X2C-1 scalar relativistic Hamiltonian (Eq. ( 28)), the one-body part of the X2C Hamiltonian derived with the X2C-N scheme and the two-body part of the Breit-Pauli Hamiltonian. This involves terms in Eqs. ( 31) and ( 16).

• "x2cn-x2c" uses the X2C-1 scalar relativistic Hamiltonian (Eq. ( 28)), the one-body part of the X2C Hamiltonian derived with the X2C-N scheme and the two-body part of the X2C Hamiltonian. This involves the terms in Eqs. ( 31) and ( 33).

• "x2c1-bp" uses the X2C-1 scalar relativistic Hamiltonian (Eq. ( 28)), the one-body part of the X2C Hamiltonian derived with the X2C-1 scheme and the two-body part of the Breit-Pauli Hamiltonian. This involves the integrals in Eqs. ( 30), ( 16).

• "x2c1-x2c" uses the X2C-1 scalar relativistic Hamiltonian (Eq. ( 28)), the one-an two-body parts of the X2C Hamiltonian derived with the X2C-1 scheme. This involves the integrals in Eqs. (30), and (33).

In all calculations we begin by performing stateaverage HCISCF [START_REF] Smith | Cheap and near exact CASSCF with large active spaces[END_REF] (which is a CASSCF-like method where FCI is replaced by HCI) calculations targeting the relevant states with an equal weight. These HCISCF calculations are performed with a scalar relativistic Hamiltonian and SOC terms are not included at this stage. The optimized active space and orbitals are then used to perform a one-step SHCI calculation where both the scalar relativistic and SOC terms are included. Thus the spinorbit coupling and electron correlation are treated on an equal footing during the correlated calculation. In some calculations the initial HCISCF calculation and the subsequent one-step SHCI calculation are performed with a different active space. This is reasonable because the aim of the HCISCF calculations is to generate orbitals that are equally weighted towards all the targeted states; the relevant SOC calculation is the one-step SHCI calculation.

All calculations were performed with a combination of the PySCF 99 and Dice codes. The BP and X2C integrals and the routines for performing HCISCF calculations are implemented in PySCF and the SHCI calculations are performed using the Dice code. The SOC integrals are implemented using the automated code generator of the libcint library 100 , which provides an efficient implementation of the standard integrals and their nuclear derivatives.

In the following it is important to note that SHCI has stochastic error associated with it, which is sometime shown in parentheses, but is in any case much smaller than the targeted zero-field splitting.

V. RESULTS

We present results of calculations done on the halogen group atoms (F, Cl, Br, I), the coinage metals group atoms (Cu, Au) as well as on the Neptunyl(VI) dioxide radical, NpO 2+ 2 , using the different Hamiltonians described in Section IV.

A. Halogens

We begin by performing a state-average HCISCF calculation where we target the three doubly-degenerate ground states corresponding to the singly occupied p x , p y or p z orbitals using a valence active space of 4 electrons in 7 orbitals with the ANO basis set 102 . We then perform a one-step SHCI calculation including SOC terms to calculate the six lowest energy states. We find a 4-fold degenerate 2 P 3/2 ground state and a 2-fold degenerate 2 P 1/2 excited state, and their energy difference corresponds to the zero-field splitting measured experimentally and tabulated in the NIST reference table of Ref. 101 (displayed in the last column of Table I).

It is important to note that the Hilbert space corresponding to the larger active spaces in Table I can be enormous, for example the number of determinants that contribute to the ground and excited state of the Iodine TABLE I. Zero-field splitting results (cm -1 ) for the atoms of the halogen group using SHCI with the ANO basis and using several active spaces and spin-orbit coupling methods. The stochastic errors on the zero-field splitting are systematically lower than 1 cm -1 and the numbers in brackets are our estimates of the error in ZFS relative to the fully converged results calculated using the extrapolation scheme (see text). The results are compared to previously reported data from Ref. 61 atom with an active space of (121o,17e) is greater than 10 24 (this number is 242 C 17 , which assumes that no symmetry other than the particle number symmetry is used). With this kind of enormous active space, the individual energies of the targeted states are not fully converged, however, the zero-field splitting is converged to better than 1% of its FCI value except in the case of Br where it is accurate to 3%. The estimated error in the energies is calculated using the extrapolation technique described in Ref. 97, which uses the fact that near convergence the SHCI energies are nearly perfectly linear with respect to the SHCI perturbative correction. Such an extrapolation FIG. 1. Convergence of the total energies (Ha+2604) of the 2 P 1/2 and 2 P 3/2 states (blue and red dots in the upper part of the graph) and the zero-field splitting (cm -1 +3685, black dots in the lower part of the graph) with the SHCI perturbative correction. The dotted lines are linear fits through these points which are used to extrapolate the energies to vanishing perturbation energies. There is a small noise (less than 1%) associated with the ZFS which can be traced to a combination of stochastic noise and uncertainty of the SHCI energies. is shown in Figure 1, where this trend seems to hold true in presence of the SOC terms as well. Such an extrapolation is used to estimate the error in the ZFS energies calculated using SHCI and is given in the bracket.

From Table I it is possible to compare the performance of the various SOC Hamiltonians. We find that the "x2c1-x2c" scheme delivers the most accurate ZFS energies for most of the atoms. Interestingly, the simple Breit-Pauli results are extremely accurate for lower atomic weight species including Florine and Chlorine, however the results progressively start to deteriorate, significantly over-estimating the ZFS with the increase in the molecular weight of the atom. This is to be expected because the Breit-Pauli SOC Hamiltonian is unbounded and thus over-estimates the energy splittings of a heavy atom in a variational calculation [START_REF] Dyall | Interfacing relativistic and nonrelativistic methods. I. Normalized elimination of the small component in the modified Dirac equation[END_REF] .

The disagreement between the experimental and "x2c1-x2c" ZFS energies can be attributed to three causes, basis set incompleteness error, neglect of corecorrelation and error in the relativistic Hamiltonian. For the Fluorine atom the core is fully correlated and it is surprising to note that the agreement with experiment becomes worse as we increase the active space to include the core electrons. This has to be attributed to basis-set incompleteness and deficiencies of the relativistic treatment, which include neglect of various two-body scalar relativistic and spin-orbit coupling terms and the use of the SOMF approximation. At the moment it is not possible to disentangle the two because the current implementation of SHCI in Dice only works with the one-body SOC terms, however, in future publications we intend to address this point in more detail. Similar errors can be seen with other Halogens as well. In the cases of Bromine and Iodine, the improvement in agreement with the experiment when one increases the active space electrons from 7 to 17 suggests that including inner core electrons in the active space might further improve the agreement with experiments.

B. Coinage metals

We have performed SHCI calculations on the more challenging coinage elements, Cu and Au. In all these calculations we again first perform a scalar relativistic state-average HCISCF calculation using Roos's ANO basis set 105 with the 11 valence electrons in 11 orbitals (including the valence s and d orbitals along with the virtual d orbitals to account for the double-d shell effect) to simultaneously optimize with equal weight the 6 lowest lying states, 5 of which are degenerate, where the five 3d and the 4s orbital are singly occupied. The obtained optimized orbitals are then used to perform a onestep SHCI calculation with the "x2c1-x2c" SOC Hamiltonian to obtain a 2-fold degenerate 2 S 1/2 ground state, a 6-fold degenerate 2 D 5/2 state and a 4-fold degenerate 2 D 3/2 state. Two such calculations were performed, one in which the same (11o,11e) active space was used and another in which all the virtual orbitals were included in the active space in addition to semi-core electrons. We report the result of these calculations in Table II which also contains the reference data from Ref. 104 and results of previous theoretical calculations.

It is interesting to note (see Figure 2) that although the ZFS splitting of the 2 D state into the 2 D 5/2 and 2 D 3/2 states, calculated by SHCI and other methods, are quite accurate, the excitation energies of these states relative to the 2 S ground state are very different. The agreement between the SHCI calculations and experiments is quite good with the maximum error being of just 0.03 eV, as opposed to other methods including the large active space DMRG-SISO 64 calculations where the error can be as large as 0.12 eV. This difference in accuracy can be attributed to the use of the ANO-TZ basis set and the use of perturbative treatment of the SOC terms in the DMRG-SISO calculations. We note that the accuracy of CASPT2-SO 103 is significantly better than CASSCF-SO 103 , which points to the fact that dynamical and core correlations are important. The ability to offer insight on magnetic and chemical properties of lanthanide and actinide containing complex 106 . We begin by performing a state-average HCISCF calculation with a (4o,1e) active space in which the energies of the doubly degenerate 2 Φ and 2 ∆ states are optimized using the X2C scalar relativistic Hamiltonian and we find an energy splitting equal to 2108 cm -1 . The energy difference between these two states is indicative of the strength of the crystal field splitting since in the absence of the axial oxygen atoms and relativistic effects, these states would be exactly degenerate. Next we perform the onestep SHCI calculations with two different active spaces, (4o,1e) and (143o,17e), to obtain the 4 doubly degenerate states which are shown in Table III. It is interesting to note that with the use of the same active space and geometry, the ZFS calculated here using the X2C Hamiltonian is significantly different from the one calculated by Gendron et al. 106 , where the DKH Hamiltonian along with atomic mean field spin-orbit coupling as implemented in Molcas was used. The ZFS splitting changes very little relative to this result even when slightly different geometry was used by Knecht et al. [START_REF] Knecht | A Nonorthogonal State-Interaction Approach for Matrix Product State Wave Functions[END_REF] . The difference between our results and those of Gendron and Knecht can be attributed to the difference in the Hamiltonian used to treat the relativistic effect. For the Np atom, it will certainly be interesting to perform a full four-component calculation with the DCB Hamiltonian to calculate the ZFS and compare the accuracy of the various approximate twocomponent Hamiltonians. Such an MRCI calculation was performed by Knecht, however only the calculated g-tensors was reported, and was in good agreement with the DKH results calculated there. It should be pointed out that even though the energy differences are quite different, we find that the ground state energy calculated with the (4o,1e) active space consists of 89% 2 Φ state and 11% 2 ∆ state, in agreement with the results of Knecht 2 , relative to the SOC ground-state energy. On the left are spin-free and spin-orbit coupling states calculated using SHCI and the "x2c1-x2c" method, and an active space of (4o,1e) (black) and (143o,17e) (red), and on the right are reference data from Ref.

VI. CONCLUSION

In this work we have extended the SHCI algorithm to treat the relativistic Hamiltonians containing spin-orbit coupling terms, by allowing one-body integrals and configuration interaction coefficients to be complex-valued numbers. The agreement with experimental results for several atomic species has been shown to be quite good and superior to previously calculated results. This is because not only are we treating large active spaces but we are also treating the electron correlation and relativistic effects on an equal footing. The main shortcoming of the current methodology is that we cannot effectively include core correlation and treat larger systems because we do not have the ability to add dynamical correlation. Work in this direction is under way and we are looking for ways to combine SHCI with the recently-developed multireference linearized coupled cluster theory [108][109][110] , which is formulated as a perturbation theory and uses Fink's partitioning 111,112 of the Hamiltonian. We are also looking to extend the methodology so that the relativistic effects can already be included at the self-consistent field cycle to obtain complex-valued spinors. This will reduce the burden on the correlated calculations because a large part of the relativistic effects will be treated at the SCF level.

FIG. 2 .

 2 FIG.2. Errors (eV) in the individual 2 D 5/2 and 2 D 3/2 states of the coinage atoms (bars), as well as of the ZFS between the two (black line) of different methods previously reported and of SHCI calculations from this work. Although the ZFS is well reproduced by all methods, SHCI performs well in its description of the individual states.

  et al.. This indicates that the g-tensors calculated using SHCI would most likely have agreed well with those of DKH and four-component MRCI calculations performed by Knecht et al..

FIG. 3 .

 3 FIG. 3. Energies of the electronic states of NpO 2+2 , relative to the SOC ground-state energy. On the left are spin-free and spin-orbit coupling states calculated using SHCI and the "x2c1-x2c" method, and an active space of (4o,1e) (black) and (143o,17e) (red), and on the right are reference data from Ref. 106.

  FIG. 3. Energies of the electronic states of NpO 2+2 , relative to the SOC ground-state energy. On the left are spin-free and spin-orbit coupling states calculated using SHCI and the "x2c1-x2c" method, and an active space of (4o,1e) (black) and (143o,17e) (red), and on the right are reference data from Ref. 106.

  , 68 and 76 and from experimental data from Ref. 101.

	Method		Active Spaces	Ref.
	F	(4o,7e) (87o,7e)
	bp-bp	405	399		403 a
	x2c1-bp	404	397		398 b
	x2cn-bp	404	398		421 c
	x2c1-x2c 405	399		404 d
	x2cn-x2c 405	399
	Cl	(4o,7e) (95o,7e) (100o,17e)
	bp-bp	834	805	873	823 a
	x2c1-bp	825	804	867	876 b
	x2cn-bp	822	802	858	907 c
	x2c1-x2c 827	806	866(1)	882 d
	x2cn-x2c 825	804	865
	Br	(4o,7e) (95o,7e) (100o,17e)
	bp-bp	3680	3655 3797	3404 a
	x2c1-bp 3407	3382 3432	3649 b
	x2cn-bp 3373	3346 3396	3723 c
	x2c1-x2c 3428	3403 3454(93)	3685 d
	x2cn-x2c 3394	3366 3415
	I	(4o,7e) (116o,7e) (121o,17e)
	bp-bp	8816	9346 9752	6961 a
	x2c1-bp 6951	7199 7453	7755 b
	x2cn-bp 6816	7049 7246	7752 c
	x2c1-x2c 7021	7277 7487(62)	7603 d
	x2cn-x2c 6886	7127 7379

a Ref. 61 b EOM-CCSD(SOC) from Ref. 68 c X2Cmmf-FSCCSD from Ref. 76 d Exp. from Ref. 101

TABLE II .

 II Zero-field splitting results (eV) for the atoms of the coinage metal group using SHCI with the ANO basis and the "x2c1-x2c" spin-orbit coupling method. We report also some previous theoretical work, from Ref. 103 and Ref.64.

	States		Active Spaces		Ref. 104
	Cu	CASSCF-SO 103 CASPT2-SO 103 DMRG-SISO 64	SHCI
		(11o,11e)	(11o,11e)	(45o,19e)	(172o,19e)
	2 D 5/2	1.49	1.43	1.31	1.39	1.39
	2 D 3/2	1.75	1.69	1.57	1.67	1.64
	2 D 5/2 -2 D 3/2	0.26	0.26	0.26	0.28	0.25
	Au	CASSCF-SO 103 CASPT2-SO 103 DMRG-SISO 64	SHCI
		(11o,11e)	(11o,11e)	(57o,43e)	(150o,25e)
	2 D 5/2	1.71	0.97	1.02	1.15	1.14
	2 D 3/2	3.22	2.49	2.55	2.64	2.66
	2 D 5/2 -2 D 3/2	1.51	1.52	1.53	1.49	1.52

TABLE III .

 III Relative energies of the electronic states of NpO 2+ 2 calculated with "x2c1-x2c" SOC method for two different active spaces. We have also shown the results obtained previously using the atomic mean field spin-orbit Hamiltonian as implemented in Molcas and the SOC terms treated perturbatively. , the Neptunyl(VI) dioxide radical 65,106,107 using the ANO basis set truncated to the triple zeta level. The molecule has C ∞v point group symmetry and the Np-O bond is of length 1.70 Å, taken from the work of Gendron et al.

	States	Energies		Ref. 106
		(4o,1e)	(143o,17e)
	2 Φ 5/2	0	0	0
	2 ∆ 3/2	4527	3857	3011
	2 Φ 7/2	8568	8675	8092
	2 ∆ 5/2	10659	10077	9192
	found in single molecular magnets is important because
	experiments on those compounds prove difficult due to
	their high toxicity. Here, we perform relativistic calcula-
	tions using the SHCI algorithm to calculate the ground
	and low lying excited states of NpO 2+ 2

ACKNOWLEDGMENTS

The research was supported through the startup package from the University of Colorado, Boulder. We would also like to thank several useful discussions and insightful comments from Zhendong Li, Qiming Sun and Juergen Gauss.


Schwilgué's calculating machines A collection of articles

Denis Roegel

In 2002, while doing research on the history of the astronomical clock built by Jean-Baptiste Schwilgué in the cathedral of Strasbourg between 1838 and 1843, I became aware of the existence of several calculating devices designed by Schwilgué at about the same time. Several of these mechanisms had been in the collections of the Strasbourg museums, but when I asked the curators about them in 2004, they did not know where they were or even if they did still exist.

In the meantime, it occurred to me in 2003 that one of the mechanisms, a small adding machine with keys, was in fact predating what was commonly thought to be the oldest key-driven adding machine, namely Parmelee's machine patented in 1850. Schwilgué had patented his machine in 1844, six years before Parmelee. But it was only in 2006 that I found out that one of these machines was located in Zurich. I examined it in 2007, and at the same time a very similar machine, but bearing the name of Schilt, was dismantled for me by David Todd in Washington. These researches led to the publication of the first description of Schwilgué's machine (other than the patent), as well as to the first recognition of what I then thought was the earliest key-driven calculating machine.

The article was published in 2008, but in 2009 one Strasbourg museum located a number of items in a long forgotten cave, and among these items were several of Schwilgué's machines, in particular three small adding machines. I examined these machines shortly afterwards, but a summary of all known such machines was only published in 2015, owing in particular to my work on the LOCOMAT project (locomat.loria.fr). By that time, however, I knew that Schwilgué's machine was not the earliest key-driven adding machine, since Silvio Hénin had published his discoveries about Torchi's 1834 machine in 2010. Nevertheless, Schwilgué's machine still is, to my knowledge, the earlist extant key-driven adding machine which can be dated with certainty.

In 2003, I also became aware of the existence of another peculiar adding machine. This machine has keys too, but only three, and it could nevertheless add any figure between 1 and 9. This machine is undated, but can be ascribed to Schwilgué without any doubt. It was probably designed at around the same time as the 9-keys adding machine, or perhaps slightly before. It is therefore possible that this machine, of which only one copy is known, is in fact the oldest extant key-driven adding machine, although it will probably never be possible to prove it. A description of this machine, together with other machines with reduced sets of keys, was published in 2016.

But was Torchi the first to imagine or build a key-driven adding machine? It turns out that he was not. In 2003 or 2004, I was browsing White's A new century of inventions published in 1822, and I copied several of its illustrations, then layed them aside. It was only in 2015 that I had a closer look at my copies, and to my surprise, I realized what I had overlooked, namely that White had presented a key-driven adding machine in his catalogue of inventions. This machine may never have been built, but it was conceived, and White should currently be thought as the first to have imagined such a machine, a conclusion which was published in 2016.

Apart from his two small adding machines, Schwilgué has also built a much larger and specialized adding machine, to serve as an aid in the positioning of a gear-cutting machine. I had also known about this machine since 2002, and it was one of the items which surfaced in 2009 in Strasbourg. The analysis of this machine in the Summer of 2009 eventually led to the two articles published in 2015 and 2016.

Finally, also in 1844, Schwilgué patented a mechanical counter, which, although apparently very mundane, was a precursor of a whole range of modern mechanical counters, themselves important features of 19th and 20th century mechanical calculating machines. A description of Schwilgué's counter and other related devices was published in 2017.

Schwilgué's venture in calculating devices was brief. He seems to have been working on them mainly in the 1840s, but during this relatively short time, he was able to break serious ground, by building probably the first key-driven adding machine which was actually sold for a number of years, by building a seminal mechanical counter, and by experimenting with other more specialized devices.

In this collection of articles, I hope to give a comprehensive view of this venture.

Nancy, September 21, 2017.

Les machines à calculer de Schwilgué

Un recueil d'articles

Denis Roegel

En 2002, alors que j'effectuais des recherches sur l'histoire de l'horloge astronomique de la cathédrale de Strasbourg construite par Jean-Baptiste Schwilgué entre 1838 et 1843, j'ai appris l'existence de plusieurs mécanismes de calcul conçus par Schwilgué vers la même époque. Plusieurs de ces mécanismes se trouvaient autrefois dans les collections des musées de Strasbourg, mais lorsque j'ai posé la question aux conservateurs vers 2004, aucun d'entre eux ne savait où ils se trouvaient, ni même s'ils existaient encore.

Entre-temps, j'ai réalisé en 2003 que l'un des mécanismes, une petite additionneuse à touches, était en fait antérieure à celle qui était habituellement considérée comme la plus ancienne additionneuse à touches, à savoir la machine brevetée par Parmelee en 1850. Schwilgué avait breveté sa machine en 1844, six ans avant Parmelee. Mais ce n'est qu'en 2006 que j'ai trouvé que l'une de ces machines se trouvait à Zurich. Je l'ai examinée en 2007, et à la même époque une machine très similaire, quoique portant le nom de Schilt, a été démontée pour moi par David Todd à Washington. Ces recherches ont conduit à la publication de la première description de la machine de Schwilgué (autre que le brevet), et aussi à la première reconnaissance de ce que je croyais alors être la plus ancienne machine à calculer à touches.

L'article a été publié en 2008, mais en 2009 l'un des musées de Strasbourg a localisé un certain nombre de pièces dans une cave oubliée depuis longtemps et parmi ces pièces se trouvaient plusieurs machines de Schwilgué, en particulier trois petites additionneuses. J'ai examiné ces machines peu après, mais un panorama de ces machines n'a été publié qu'en 2015, du fait de mon travail sur le projet LOCOMAT (locomat.loria.fr). À ce moment, toutefois, je savais que la machine de Schwilgué n'était pas la plus ancienne additionneuse à touches, car Silvio Hénin avait publié ses découvertes sur la machine de Torchi de 1834 en 2010. La machine de Schwilgué est tout de même, à ma connaissance, la plus ancienne additionneuse à touches existante et qui puisse être datée avec certitude.

En 2003, j'ai aussi appris l'existence d'une autre curieuse additionneuse. Cette machine a aussi des touches, mais seulement trois, et elle pouvait tout de même additionner n'importe quel chiffre entre 1 et 9. Cette machine n'est pas datée, mais peut être attribuée à Schwilgué sans le moindre doute. Elle a probablement été conçue à la même époque que l'additionneuse à neuf touches, ou peut-être légèrement plus tôt. Il est par conséquent possible que cette machine, dont on ne connaît qu'un seul exemplaire, soit en fait la plus ancienne additionneuse à touches existante, quoique cela ne pourra probablement jamais être prouvé. Une description de cette machine, ainsi que d'autres machines avec des ensembles réduits de touches, a été publiée en 2016.

Mais est-ce que Torchi a été le premier a imaginer ou construire une additionneuse à touches ? Il s'avère que non. En 2003 ou 2004, j'avais parcouru le livre A new century of inventions publié par White en 1822, et copié quelques unes de ses illustrations, puis je les avais mises de côté. Ce n'est qu'en 2015 que j'ai à nouveau regardé ces copies de plus près et, à ma grande surprise, j'ai réalisé ce qui m'avait échappé, à savoir que White avait présenté une additionneuse à touches dans son catalogue d'inventions. Cette machine n'a peut-être jamais été construite, mais elle a été imaginée, et White devrait aujourd'hui être considéré comme le premier à avoir conçu une telle machine, une conclusion qui a été publiée en 2016.

Outre ses deux petites additionneuses, Schwilgué a aussi construit une machine bien plus grande et spécialisée, dans le but de faciliter le positionnement d'une machine à tailler les engrenages. J'avais aussi appris l'existence de cette machine en 2002, et c'était l'une des pièces qui ont fait surface en 2009 à Strasbourg. L'analyse de cette machine au cours de l'été 2009 a par la suite conduit aux deux articles publiés en 2015 et 2016.

Finalement, aussi en 1844, Schwilgué a breveté un compteur mécanique qui, quoique d'apparence banale, a été un précurseur de tout un éventail de compteurs mécaniques modernes, eux-mêmes d'importants éléments des machines à calculer mécaniques des 19 e et 20 e siècle. Une description du compteur de Schwilgué et de mécanismes apparentés a été publiée en 2017.

La contribution de Schwilgué en matière de mécanismes de calcul a été brève. Il semble y avoir travaillé principalement dans les années 1840, mais pendant cette période relativement courte, il a été en mesure de produire des résultats conséquents, en construisant probablement la première additionneuse à touches qui a été vendue pendant un certain nombre d'années, en construisant un compteur mécanique fondateur, et en expérimentant avec d'autres appareils plus spécialisés.

Dans ce recueil d'articles, j'espère donner une vue d'ensemble de cette contribution.

Nancy, 21 septembre 2017.

Content Sommaire

The articles listed in the table of contents should be included here at the indicated pages, starting at page 7, and keeping the original odd pages on the present odd pages.