Hyperion: A New Computational Tool for Relativistic Ab Initio Hyperfine Coupling.

Herein we describe Hyperion, a new program for computing relativistic picture-change-corrected magnetic resonance parameters from scalar relativistic active space wave functions, with or without spin-orbit coupling (SOC) included a posteriori. Hyperion also includes a new orbital decomposition method for assisting active space selection for calculations of hyperfine coupling. For benchmarking purposes, we determine hyperfine coupling constants of selected alkali metal, transition metal, and lanthanide atoms, based on complete active space self-consistent field spin-orbit calculations in OpenMolcas. Our results are in excellent agreement with experimental data from atomic spectroscopy as well as theoretical predictions from four-component relativistic calculations.

Introduction

Magnetic resonance spectroscopy techniques, such as electron paramagnetic resonance (EPR) and nuclear magnetic resonance (NMR), are capable of providing very accurate information on the interactions between electron spins, nuclear spins, and external magnetic fields. The information available from EPR and NMR spectra is often encoded via iterative fitting of the experimental data to a model spin Hamiltonian in a set of effective parameters, each related to a specific type of coupling between magnetic entities. Quantities probing interactions between electron and nuclear spins, namely, EPR hyperfine coupling constants (HFCCs) and paramagnetic NMR shifts, depend strongly on unpaired electron (spin) density[1,2] and are central to the study of various phenomena encountered in molecular systems of interest. Examples include using HFCCs as proxies for covalency in actinide complexes,[3] quantifying spin decoherence in molecular qubits via theoretically determined hyperfine coupling (HFC) tensors,[4] and employing contact NMR shifts—which have an intrinsic dependence on HFC parameters[5]—to understand the solution structure and magnetic properties of paramagnetic MRI contrast agents.[6]

The steady evolution of electronic structure algorithms has made it computationally feasible to study heavy-element complexes fully ab initio. As the atomic number increases, relativistic effects become more important, to the extent where they cannot be regarded as mere perturbations of the Schrödinger picture; instead, a four-component Dirac formalism becomes necessary. The onset of the relativistic regime is especially important for magnetic interactions, which couple the electronic and positronic degrees of freedom of a Dirac spinor, therefore requiring an explicit description of the latter. Hence, theoretical frameworks developed for magnetic properties such as HFC must account for relativistic effects—both spin-independent (scalar relativistic, SR) and spin-dependent (e.g., spin–orbit coupling, SOC)—to ensure a wide range of applicability and to keep up with the latest experimental advances.

HFC between an unpaired electron spin S and a nuclear spin IN is most frequently modeled non-relativistically as the sum of anisotropic dipolar coupling, known as the spin-dipole (SD) mechanism, and isotropic Fermi coupling (FC)[7]where we use μ0, ge, μB, gN, μN, and IN to denote the vacuum permeability, electron g-factor, Bohr magneton, nuclear g-factor, nuclear magneton, and nuclear spin vector, respectively; rN = r – RN is the position vector of the unpaired electron with respect to the magnetic nucleus. The strength of FC is proportional to the spin population at the magnetic nucleus, ρNα–β, thus providing a convenient probe for spin delocalization in a molecule. Although the single-configurational view of electronic structure suggests that this term arises solely through s-type atomic orbital (AO) contributions to the singly occupied molecular orbital (SOMO), several factors complicate this interpretation.

Differences between the interaction of core spin-up and spin-down electrons with the unpaired spin give rise to spin polarization (SP),[2] which gives additional contributions to the isotropic HFC and cannot be described by a single electronic configuration. Thus, the quantitative interpretation of HFCCs requires quantum chemical techniques to model the electronic structure accurately and capture SP effects. For this purpose, there are two possible solutions: spin-unrestricted single-configurational methods, such as unrestricted Hartree–Fock (UHF) and unrestricted Kohn–Sham density functional theory (UKS-DFT), which offer low computational cost at the expense of broken spin symmetry, and more expensive spin-adapted multiconfigurational methods, which preserve spin symmetry and are designed to handle electron correlation.

Unrestricted DFT is currently the most widely used approach for determining HFCCs, due to the lower computational cost relative to wave-function-based (ab initio) algorithms. Hybrid functionals give good predictions for HFCCs of organic radicals and transition metal complexes;[8,9] however, the accuracy of these results is believed to be caused by fortuitous error cancellation.[10] Moreover, a recent study[11] shows that the best choice of functional for HFC is system-dependent. Aside from HFC-related shortcomings, the single-configurational framework underlying Kohn–Sham DFT is inappropriate for describing static correlation, an important feature of f-element complexes. Such systems require a multiconfigurational approach, usually in the form of active space wave function optimization techniques such as complete active space self-consistent field (CASSCF), in order to obtain meaningful predictions of energies and molecular properties.

Although far from black-box, active space algorithms are not only cheaper than fully correlated alternatives such as full configuration interaction or coupled cluster methods but also more flexible, as the active space is user-defined. This framework can be leveraged to obtain accurate theoretical HFCCs by including, in addition to static correlation, a selection of dynamical correlation effects (e.g., SP) that significantly influence HFC. The challenge, then, is developing a strategy for choosing computationally feasible active spaces that result in accurate theoretical HFCCs. Previous studies employing multireference configuration interaction (MR-CI) algorithms have already analyzed the convergence of HFCCs with respect to CI excitation level, orbital selection threshold, and basis set completeness.[2,12−19] However, most of the established trends are only applicable to the non-relativistic regime; in particular, the observations about AO contributions to HFCCs result directly from the delta distribution form of the FC term (eq , right), which is specific to the non-relativistic Schrödinger–Pauli framework.[20] Relativistic treatments of HFC hence require new work to establish updated guidelines for theoretical investigations.

Fully relativistic four-component approaches are unfeasible for all but the simplest systems,[21,22] unless combined with a low-cost electronic structure algorithm such as Dirac–Hartree–Fock (DHF)[23] or DFT.[24] As such, multiple strategies that decouple the upper (electronic) and lower (positronic) components of the Dirac Hamiltonian have been developed to lower the cost of relativistic calculations for application to real molecular systems. Use of a decoupling transformation is, in effect, a change in the reference frame of the wave function, and therefore, it must also be applied to the property operators; ignoring this second step leads to the so-called picture-change error (PCE).[25] Due to the picture-change correction, relativistic property operators are different from their non-relativistic counterparts and, in the case of relativistic HFC, the isotropic contribution is no longer proportional to the spin density at the nucleus.

Approximate (quasi-relativistic) decoupling techniques, such as the regular approximation (RA) and Douglas–Kroll–Hess (DKH), are now fairly widespread and provide excellent predictions for the energies of most relativistic systems.[26,27] These two-component approaches can be simplified further by disregarding spin-dependent effects, yielding a one-component scalar relativistic formalism. For magnetic resonance applications, the zeroth-order regular approximation (ZORA) is widely used in combination with DFT.[28−31] Despite providing a good description of valence properties,[32,33] ZORA affords large errors for core-dependent properties of heavy-element systems (e.g., core-hole X-ray excitations or absolute nuclear shielding tensors);[32,34] DKH is significantly more accurate in this respect. Indeed, HFCCs obtained using SR-DKH2 operators are in good agreement with experiment,[35,36] although Nguyen Lan et al.[37] show that a higher-order decoupling transformation (DKH3) is needed for a converged relativistic HFC picture. Quasi-relativistic HFCCs can also be derived via the infinite-order regular approximation (IORA), and when combined with a multiconfigurational wave function method, this approach produces highly accurate atomic hyperfine structure constants for alkali and coinage metals.[38]

Exact two-component decoupling schemes are also available, together with one-component SR variants. For the purpose of theoretical HFC, these offer two main improvements over quasi-relativistic theories: there is no uncertainty regarding the appropriate order of decoupling, and picture-change corrections are straightforward to implement via matrix multiplication. The exact-2-component (X2C)[39−41] and normalized elimination of the small component (NESC)[42] approaches are two popular, fully numerical, choices for theoretical studies of relativistic HFC. In such cases, it is worth keeping in mind that the form of the relativistic HFC operator changes on a case-by-case basis, since the picture-change transformation is only defined in matrix form.

Our goal was to devise a general methodology for determining relativistic, picture-change-corrected HFCCs for chemical systems of arbitrary size and complexity. With this in mind, we developed Hyperion, a Python-based program that computes SR-X2C-decoupled magnetic resonance parameters from complete active space (CAS) or restricted active space (RAS) wave functions, with or without spin–orbit coupling (SOC) added a posteriori (CASSCF-SO/RASSCF-SO); a method based on the same theoretical formalism has been developed in parallel by Autschbach and co-workers and implemented as part of OpenMolcas.[41] Although our code has to date been tested only with OpenMolcas, Hyperion is a stand-alone package that can be straightforwardly extended to allow inputs from other quantum chemical software. Herein, we demonstrate the use of Hyperion to obtain HFCCs of selected atoms and benchmark our results against experimental data from atomic spectroscopy and predictions from four-component calculations. As well as showing the performance of Hyperion, we demonstrate various strategies for tackling specific combinations of electron correlation, SR, and SOC effects, highlighting both merits and limitations of CASSCF-SO/RASSCF-SO methods. The comparatively small number of electrons and high symmetry of atoms mean that high levels of theory are achievable with relatively little computational cost, but as our goal is to extend our approach to molecules in the near future, we work with usual electronic structure approximations, such as contracted basis sets and RAS subspaces of limited size, which are the only feasible strategies for molecular calculations.

Theory

SR-X2C Magnetic Properties

Here, operators are denoted using the hat symbol (T̂) and vector quantities are shown in bold (A). One-component and two-component operators are represented using normal font (ĤSF-X2C,UU(1)), while bold-italic notation (SF-X2C,UU(1)) is used for their matrix representations in a scalar AO basis. We denote four-component operators as matrix operators (ĤmDE) to emphasize their 2 × 2 block structure. Finally, we employ double-struck notation () for matrix representations of four-component operators.

The starting point for any relativistic treatment of magnetic properties is the four-component Dirac equation under a scalar potential V and a vector potential Awhere σ is the 3-vector of Pauli spin matrices, p̂ is the electron linear momentum operator, me is the electron mass, and c is the speed of light in a vacuum; SI units are used throughout this work. We represent all four-component operators in block form, with each matrix element denoting a two-component operator. Note that, in two-component equations, σ = 2s, where s is the electron spin vector.

The Dirac wave function, also known as a four-component spinor, has an upper component ψU, sometimes referred to as the large component, as well as a lower component (or small component) ψL. In order for eq to approach the correct non-relativistic limit, the basis sets chosen for ψU and ψL must obey the restricted kinetic balance (RKB) condition,[43−46] namely, that, for an upper component basis set {ϕμ}, the lower component basis set is . Although using RKB alone is formally justified only in the absence of magnetic fields, this approach is reasonable, as the magnetic-field-dependent terms are treated as a perturbation herein (vide infra). By substituting the RKB condition into the four-component Hamiltonian, we arrive at the modified Dirac equation[47,48]withNotice that the lower component of eq has been replaced by the pseudo-upper component χU, which can be represented in the same basis as that used for ψU. Additionally, the modified Dirac Hamiltonian (eq ) has been separated into A-independent and A-dependent contributions; we henceforth treat the A-dependent term as a first-order perturbation, Ĥ(1). Using the Dirac relation, we rewrite Ŵ asand discard the σ-dependent term to obtain a spin-free four-component operator. By ignoring the spin-dependent contribution, this relativistic formalism becomes compatible with the one-component framework employed by most electronic structure packages, including OpenMolcas. The resulting spin-free (SF) modified Dirac operator iswhich shall serve as the zeroth-order (unperturbed) Hamiltonian within our theoretical framework.

Scalar relativistic exact-2-component (SR-X2C) theory[49−52] can be employed to determine a unitary transformation matrix (eq ) that decouples the upper and pseudo-upper components by block-diagonalizing the matrix representation of ĤSF-mDE(0), herein denoted as . The electronic problem can then be solved by diagonalization of the upper–upper matrix, SF-X2C,UU(0), without any reference to the other two-component blocks.The SF-X2C-decoupled perturbation is obtained in an analogous manner, by applying the same unitary transformation to . For simplicity, we employ operator notation throughout the remainder of this section; however, it should be noted that, within the Hyperion implementation, the decoupling transformation of Ĥ(1) and indeed all subsequent steps in the calculation of expectation values are carried out in matrix form. As this work pertains to electronic properties, only the upper–upper block of the transformed operatoris of interest, and we henceforth drop the UU subscript from Hamiltonian operator notation.

The interaction between an unpaired electron and a nuclear spin, i.e., HFC, can be derived from eq by substituting in the vector potential AN induced by a point-like magnetic dipole NAlternatively, the nuclear magnetic dipole can be modeled as a Gaussian distribution (eqs and 13)—this is also implemented in Hyperion; however, results herein employ the point nucleus expression for simplicity.The electronic HFC perturbation operator within the SF-X2C framework is thereforewhere we distinguish between a spin-dependent contribution, ĤN,SF-X2CFC+SD, and an imaginary, spin-independent contribution, ĤN,SF-X2CPSO. The former reduces to the sum of the FC and SD operators in the non-relativistic limit, while the latter represents the interaction between the electronic orbital angular momentum and the nuclear spin and is known as the paramagnetic spin–orbit coupling (PSO) term.

It is worth highlighting that eq does not include a delta-function, that would imply only contributions at the nucleus and characteristic of traditional interpretations of FC; as shown by Kutzelnigg,[20] this hallmark only arises in the non-relativistic limit of a two-component framework. As a result, relativistic spin-dependent HFC cannot be interpreted as a combination of classical dipolar coupling (SD term) and an FC term sampling the spin density at nucleus N.[7]

An important aspect that bears discussion is the spin–orbit contribution to HFC. Herein, SOC effects are included in the zeroth-order Hamiltonian, and as such, the spin–orbit HFC contribution is a first-order property,[34] computed as the expectation of the PSO operator, eq . An alternative approach involves modeling SOC as a perturbation and mapping the HFC operator onto a true spin effective Hamiltonian,;[53] in this case, the spin–orbit HFC term only appears in the second-order response, and the second-order HFC energy includes both the PSO interaction and SOC effects on the wave function.[54] Therefore, even though both approaches yield spin–orbit contributions to the HFC energy, the two are not directly comparable. The perturbative SOC formalism is the de facto choice in non-relativistic HFC studies;[9,31] meanwhile, relativistic four- and two-component methods already account for SOC in the wave function; hence, the HFC energy includes a first-order PSO term.[55] Within a one-component formalism such as that employed by OpenMolcas and Hyperion, SOC can be modeled using either approach. Nevertheless, as we are targeting a variety of chemical systems, including atoms and molecules where strong SOC significantly alters the manifold of electronic states, we do not implement the perturbative SOC method within Hyperion. Instead, effects of SOC are included a posteriori in the zeroth order wave function via the restricted active space state interaction (RASSI) formalism[56] before calculation of the hyperfine coupling.

Spin Hamiltonian Parametrizations

The cornerstone of EPR theory is the spin Hamiltonian, an effective operator describing all of the relevant magnetic interactions as couplings in a (2S + 1)-dimensional space, where S is the spin quantum number of the system. The HFC spin Hamiltonian iswhere the HFC tensor N quantifies the strength and anisotropy of the interaction. Quantum mechanical determinations of EPR parameters are based on mapping ab initio operators onto the spin Hamiltonian and differentiating to obtain the tensor components of N:where . In a spin-free/scalar relativistic framework, provided the electronic system is spin-only (no SOC), the tensor components are computed aswhere we have defined the vector operator ĥN asNote that in the case of spin-free electronic structure approaches, such as spin-free CASSCF, the spin projection expectation value ⟨Ŝ⟩ is exactly equal to S.

This spin-only parametrization is not appropriate, however, for systems exhibiting non-negligible SOC, as S is no longer a good quantum number. Instead, the size of the model space is chosen as the number of low-lying electronic states in the SO-coupled energy spectrum, resulting in an effective multiplet, where is the pseudospin quantum number. Replacing S with in eq yields the pseudospin Hamiltonian, from which the elements of symmetrized tensor NNT are computed via the method of Chibotaru[57]where |I⟩ and |J⟩ denote eigenstates in the pseudospin manifold. A similar approach is widely used to derive g-tensors for strongly SO-coupled systems.[58] The eigenvectors and eigenvalues of the symmetrized tensors correspond to the principal axes and the squared principal values of the original tensors; as a result, the pseudospin parametrization yields unsigned HFCCs. Although the present work is focused on HFCC magnitude, we note in passing that sign information becomes important if the parameters determined theoretically are employed in further computational work, such as simulations of EPR spectra.

Finally, Hyperion is also capable of determining SF-X2C g-values, using a similar methodology to that described above for HFCCs. The key expressions employed in g-value calculations can be easily obtained by replacing eqs and 17 with the vector potential induced by an external magnetic field and with the electron-Zeeman spin Hamiltonian, respectively. This additional capability will be assessed as part of a future study on molecules.

HFC Orbital Decomposition

One aspect of theoretical determination of HFCCs using multiconfigurational methods is choosing the most appropriate active space. This is largely a task of trial and error, for which it is hard to define general rules. Hence, we have developed an orbital decomposition method to assess the involvement of particular MOs in the spin-dependent HFCC of a spin-free state ΨSS. To do so, we use the second quantization formalism to represent the expectation value of each vector operator component ĥN,FC+SD as a sum over pairs of MOs {μν}where ↠and â denote spin–orbital creation and annihilation operators, respectively, and Pα–β is a spin density matrix element. We have used ÔN, to denote the l = z components of a rank-2 tensor operator defined asNote that ÔN represents the spatial part of hN,SF-X2CFC+SD, as derived via eqs and 20.

We consider each term from the summation in eq separately, combine the k = x, y, z components of ÔN, via the vector norm, and divide the result by the spin projection to obtain a two-dimensional symmetric matrix with elementsThe matrix contains the same information as the spin-dependent HFCC, albeit in a modified form wherein MO degrees of freedom are not integrated out. We note that the factor of ⟨Ŝ⟩–1, which reduces to S–1 for spin-free states, mimics the expression for HFC operators in the spin-only formalism. Within this representation, the diagonal elements can be interpreted as individual orbital contributions, while the coupling between two different orbitals is quantified by .

Inspired by density matrix renormalization group (DMRG) entanglement diagrams,[59] we designed a similar pictorial representation of orbital involvement in HFC (Figure ), using the matrix elements of . Individual orbital contributions are shown as markers on the circumference of a circle, while pairwise contributions are shown as chords, color-coded by order of magnitude.

Figure 1

HFC orbital decomposition diagrams obtained from relativistic determinations of HFCCs using Hyperion. (a) 133Cs atom, RASCI(45,52) using CASSCF(11,15)-optimized orbitals; (b) 197Au atom, RASCI(45,54) using RASSCF(19,26)-optimized orbitals. Lines across the diagram correspond to off-diagonal matrix elements of . Individual orbital contributions, as indicated by circular markers, correspond to diagonal matrix elements .

It is worth emphasizing that this model is based on a CI-type wave function that is not relaxed after applying the first-order HFC perturbation and that our orbital decomposition diagrams only show a static picture of HFC. To further understand this, we adopt the direct/indirect terminology used by Engels[14] to describe the influence of specific configurations in a CI wave function. The HFC orbital decomposition matrix for a given electronic state is computed via its spin density matrix, which is completely determined by the CI expansion. The direct influence of an (active) orbital or pair of orbitals, as indicated by , is then proportional to the total CI contribution from all configurations where the orbitals of interest are singly occupied. However, the spin density associated with one orbital can vary between different active space selections, due to the indirect effect of all of the other correlated orbitals. Therefore, the active space selection should not be based solely on the orbital decomposition diagram. HFC calculations should be carried out both with and without a specific orbital in the CAS/RAS to determine the magnitude of its indirect effect.

Computational Details

All electronic structure calculations use a local version of OpenMolcas[60] v19.11 adapted to print the X2C decoupling matrices, while the spatial HFC integrals are evaluated analytically with Libcint.[61] Full ANO-RCC basis sets[62−66] are used throughout to ensure sufficient flexibility and to maintain consistency with previous relativistic HFC studies.[39,41,67] The electronic wave function is optimized using either the complete active space self-consistent field (CASSCF) approach[68] or a restricted active space (RAS) approach[69] with (RASSCF) or without orbital optimization (RASCI). For RASSCF/RASCI calculations, the maximum allowed number of RAS1 holes and the maximum allowed number of RAS3 electrons are both set to 2. Spherical symmetry in the optimized MOs is enforced via the ATOM keyword, which prevents rotations among orbitals belonging to different representations of the full rotation group. Where necessary, SOC is added a posteriori using the RASSI approach,[56] which employs a Pauli spin–orbit Hamiltonian,[70] represented using atomic mean field integrals (AMFIs).[71] To keep RAS calculations tractable, the RAS2 subspace is restricted to the partially occupied shells. As such, for 2S1/2 systems, only the singly occupied s orbital is included in RAS2, while, for transition metals or lanthanides with partially filled d or f shells, respectively, RAS2 is made up of the valence nd or nf and (n + 1)s orbitals as required.

Results and Discussion

Spin Polarization vs Electron Correlation

We take this opportunity to discuss a conceptual aspect that is often encountered in computational HFC literature: the overlap between spin polarization and electron correlation. Both effects have physical interpretations related to the interaction between electrons in a many-electron system; however, a distinction arises in electronic structure theory due to the use of single-configurational (also known as mean-field) references with respect to which perturbations such as SP, static correlation, and dynamic correlation are described.

Unrestricted Hartree–Fock (UHF) theory is the simplest wave-function-based electronic structure framework that accounts for SP, albeit not accurately, due to spin contamination. The UHF wave function can be expressed as a perturbation expansion from a restricted open-shell Hartree–Fock (ROHF) reference, with the first-order term comprising only single excitations.[72] It can be deduced that SP arises mainly from singly excited configurations; a spin-adapted CI ansatz with singles (S–CI) should therefore be able to provide a more accurate, non-spin-contaminated description of SP.[2]

Electron correlation—the instantaneous interaction of electrons, which is not captured by a mean-field ansatz—is also routinely described by a CI wave function, which is most accurate when all possible configurations obtained by distributing all electrons among the available orbitals are included (full CI, FCI). The correlation energy is defined as the difference between the FCI and the reference Hartree–Fock energies. For open-shell systems represented in a spin-adapted framework (i.e., with a ROHF reference), it follows that electron correlation includes an SP contribution,[72] whose magnitude directly relates to the overall weighting of single excitations in the CI expansion. Note, however, that, if a UHF reference is used instead, the boundary between SP and “other” correlation becomes unclear. This is fortunately not a concern herein, as we employ spin-adapted wave function methods exclusively.

In a comprehensive review of the spin-polarization model for HFC, Chipman[2] refers to second- and higher-order correlation effects/true correlation, thus dividing spin-polarizing single excitations from other excited configurations sampled by FCI. While this split view helps balance computational efficiency with accuracy in multiconfigurational calculations, it has also led to a prevalence in the HFC literature of single-excitation-only results, with no discussion of higher-order effects.[41,67] This is likely a safe approximation for systems such as simple organic radicals;[73] however, there is no reason to assume a priori that the SP contribution to HFC is more important than contributions from higher-order excitations. In fact, early work shows that high-level CI approximations are required even for second period atoms to achieve quantitative agreement with experiment.[16] For the atomic systems studied here, we attempt to include most correlation afforded by a chosen active space using either CASSCF or, where this is unfeasible, RASSCF/RASCI with single and double excitations.

Alkali Metals

We first apply our methodology to study the hyperfine structure of atoms in the alkali series. This is a popular test set for theoretical HFC computations, requiring a good description of electron correlation, as well as an appropriate treatment of relativistic effects for the heavier elements.[22,74,75]Table shows predicted (unsigned) HFCCs for the 2S1/2 ground state of each alkali atom, determined via active space electronic structure methods. We note that our CASSCF(1,1) results are very close to the four-component Hartree–Fock data reported by Talukdar,[22] suggesting that, at least for the direct contribution to HFC, relativistic effects are correctly accounted for with the X2C decoupling in Hyperion. In order to include correlation contributions to the spin density, the active space must be expanded beyond the SOMO; the CAS size is limited to a maximum of 18 orbitals, and we turn to RAS methods to explore electron correlation effects in a larger orbital space. Note that, due to the computational scaling of RAS algorithms, we restrict the RAS2 subspace to the SOMO. HFCCs converge the fastest when the CASSCF algorithm is used, with CASSCF(11,15) results within approximately 10% of experimental HFC data for alkali Na–Fr. This contrasts with RASSCF and RASCI results, which display slower convergence with active space size and only reach CASSCF accuracy when the majority of MOs are correlated.

Table 1

Unsigned Isotropic HFCCs in MHz Computed for the Ground State of Alkali Atomsa

atom	wave function	RAS1	RAS2	RAS3	|A|
7Li	4c-HF				288.2
	CASSCF(1,1)		1 × s		282.3
	CASSCF(3,10)		4 × s, 2 × p		367.9
	CASSCF(3,15)		4 × s, 2 × p, 1 × d		367.9
	CASSCF(3,14)		5 × s, 3 × p		395.7
	RASSCF(3,29)	1 × s	1 × s	6 × s, 7 × p	392.8
	Experimental				401.7
23Na	4c-HF				633.4
	CASSCF(1,1)		1 × s		639.2
	CASSCF(9,13)		2 × s, 2 × p, 1 × d		724.6
	CASSCF(9,14)		3 × s, 2 × p, 1 × d		794.2
	CASSCF(11,15)		4 × s, 2 × p, 1 × d		801.9
	RASSCF(11,35)	2 × s, 1 × p	1 × s	6 × s, 6 × p, 1 × d	820.9
	RASSCF(11,40)	2 × s, 1 × p	1 × s	6 × s, 6 × p, 2 × d	839.6
	Experimental				885.8
39K	4c-HF				151.0
	CASSCF(1,1)		1 × s		152.9
	CASSCF(9,13)		2 × s, 2 × p, 1 × d		184.0
	CASSCF(9,14)		3 × s, 2 × p, 1 × d		197.1
	CASSCF(11,15)		4 × s, 2 × p, 1 × d		203.3
	CASSCF(11,16)		5 × s, 2 × p, 1 × d		210.7
	RASSCF(19,39)	3 × s, 2 × p	1 × s	6 × s, 6 × p, 1 × d	203.6
	Experimental				230.8
85Rb	4c-HF				666.9
	CASSCF(1,1)		1 × s		691.9
	CASSCF(9,13)		2 × s, 2 × p, 1 × d		831.2
	CASSCF(11,15)		4 × s, 2 × p, 1 × d		906.5
	RASSCF(19,26)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d, 1 × f	859.1
	RASCI(37,54)	4 × s, 3 × p, 1 × d	1 × s	5 × s, 6 × p, 1 × d, 1 × f	908.8
	Experimental				1011.9
133Cs	4c-HF				1495.5
	CASSCF(1,1)		1 × s		1539.1
	CASSCF(9,13)		2 × s, 2 × p, 1 × d		1863.9
	CASSCF(9,14)		3 × s, 2 × p, 1 × d		2024.6
	CASSCF(11,15)		4 × s, 2 × p, 1 × d		2041.3
	RASCI(35,40)	5 × s, 4 × p	1 × s	5 × s, 4 × p, 1 × d	2083.2
	RASCI(35,47)	5 × s, 4 × p	1 × s	5 × s, 4 × p, 1 × d, 1 × f	2095.7
	Experimental				2298.1
223Fr	4c-HF				5518.0
	CASSCF(1,1)		1 × s		5240.1
	CASSCF(9,13)		2 × s, 2 × p, 1 × d		6277.0
	CASSCF(11,15)		4 × s, 2 × p, 1 × d		6685.7
	RASCI(43,57)	6 × s, 5 × p	1 × s	5 × s, 6 × p, 1 × d, 1 × f	6940.0
	Experimental				7654.0

The RAS1, RAS2, and RAS3 columns indicate the number of atomic shells—separated by angular momentum—included in each subspace. Experimental HFCCs and four-component Hartree-Fock (4c-HF) HFCCs are reproduced from ref (22).

The balance of doubly occupied and virtual (unoccupied) orbitals in the CAS/RAS has a crucial influence on HFCC accuracy; it is not sufficient to augment a minimal active space with core orbitals, as appropriate virtual orbitals are needed to correlate them. A similar observation was made in a recent coupled-cluster study[22] that emphasizes the need for high-energy unoccupied orbitals to correlate inner-core electrons. From our results, we deduce that radial correlation, introduced via virtual shells with the same angular momentum as the core–shells,[76] has the most significant effect on HFCCs. Compare, for example, the HFCCs obtained from CASSCF(9,13) and CASSCF(9,14); the virtual s shell included in the latter leads to a 6–8% improvement in accuracy. For calculations that include orbital optimization (i.e., CASSCF or RASSCF), a good basic principle for active space selection is to include one radially correlating virtual shell for each doubly occupied shell in the active space. However, the exponential scaling of CASSCF severely restricts this strategy, and as such, the largest CAS selections reported herein include only one virtual shell of each angular momentum.

It is also interesting to analyze the influence of polarization basis functions. The inclusion of one virtual d shell in CASSCF affects the HFCCs of all alkali by 4–5% (except Li), suggesting an angular correlation[76] effect related to the presence of core p electrons. Extrapolating, we postulate that, for core–shells with the highest angular momentum l, virtual shells with an angular momentum of at least l + 1 are needed to capture angular correlation. Note that, although previous theoretical work highlights the slow convergence of correlation energy with maximum angular momentum[77]—indicating that an accurate description of correlation likely requires much higher angular momenta—similar effects on HFCCs have only been explored for light atoms.[13,16] Nevertheless, testing this hypothesis here is unfeasible given the computational limitations of CASSCF and RASSCF algorithms; we therefore limit our approach to include one virtual d shell for atoms Na–Fr and additionally one virtual f shell for Rb–Fr. Nonetheless, we come within 8–12% total relative error for these atoms.

Spin-dependent HFCCs are known to be particularly sensitive to correlation effects from inner-core electrons; to investigate this, we performed RASSCF and RASCI calculations with RAS1 subspaces spanning most of the core region. All RASCI calculations are performed with CASSCF(11,15)-optimized orbitals for consistency. Orbital decomposition analysis of RASCI results (Figure a) reveals a trend of decreasing HFC contribution with increasing angular momentum. Most s and p orbitals are strongly coupled by the HFC operator and make significant (>10 MHz) direct contributions to the HFCC. The influence of (virtual) polarization functions (4–5% increase in HFCC accuracy) is not reflected by the relatively insignificant contributions (<10–3 MHz) shown in the orbital decomposition diagram. Therefore, such functions have a predominantly indirect effect, displacing spin density from orbitals with significant HFC contributions. Occupied d and f shells show similarly small orbital decomposition contributions, and additionally, our calculations suggest their inclusion leads to insignificant variations in the computed HFCC. As such, we conclude that core d and f orbitals of alkali atoms can be safely left out of the RAS1 subspace for the purpose of HFCC determinations. We note that the error with respect to experimental HFCCs plateaus around 10% for RASCI (Figure ), which could be a consequence of the restricted excitation level and/or the basis set size. These results are nevertheless very encouraging, considering the highly correlated nature of atomic systems, as well as the fact that our calculations explore a relatively limited parameter space (contracted basis sets, RAS1/RAS3 subspaces limited to no more than 35 orbitals).

Figure 2

Evolution of HFCC absolute error with active space size for alkali atoms. Absolute errors are computed with respect to experimental data from ref (22).

Lastly, the issue of orbital optimization in RAS calculations bears discussion. The super-CI algorithm employed by OpenMolcas “folds” single excitations involving inactive and secondary orbitals into the active MO coefficients,[68] thus capturing more correlation than a CI-only calculation with the same parameters. Theoretically, orbital optimization should allow for accurate HFCCs to be obtained with relatively little computational expense, since the active space needs only be large enough to capture the significant correlation contributions to HFC; indeed, this proves true for CASSCF; however, for RASSCF, the results are unpredictable. In the case of light alkali (Li–K), it was possible to run RASSCF calculations involving the entire core s and p manifold, together with polarization functions, which produced HFCCs in excellent agreement with experiment. On the other hand, this approach is prohibitive for heavier atoms—even if memory is not a limitation, orbital optimization is slow and prone to convergence issues. We therefore only employ RASCI for large active space calculations on Rb–Fr. RASSCF calculations involving fewer active orbitals were also attempted; however, the resulting HFCCs were on average less accurate than those from RASCI. It appears that the excitation level restriction severely hinders the efficiency of orbital optimization; stochastic CASSCF[78] and DMRG-CASSCF methods[79] are likely to provide more accurate HFCCs; however, due to additional uncertainty introduced by approximate CI solvers, such algorithms are not explored herein.

Coinage Metals

The treatment of HFC in coinage metals is very similar to that of the alkali metals, given they both have single-reference 2S1/2 ground states. A notable point of difference is the valence d shell which, unlike the core-like d shells in alkali atoms, imparts a significant contribution to the HFCC. Correlating the entire valence region is therefore not feasible with CASSCF, and in this case, RASSCF is the only option. We find that the benefits of a larger active space outweigh the limitations due to excitation level restrictions, leading to a net improvement in HFCC accuracy from CASSCF to RASSCF (Table ). This contrasts with the behavior observed for the alkali metals, indicating that it is perhaps the nature of the correlated orbitals (valence vs core), rather than the CI approximation, that makes the biggest difference to theoretical HFCCs.

Table 2

Unsigned Isotropic HFCCs in MHz Computed for the Ground States of Cu, Ag, and Aua

atom	wave function	RAS1	RAS2	RAS3	|A|
63Cu	CASSCF(11,6)		1 × s, 1 × d		4137.3
	CASSCF(19,14)		3 × s, 2 × p, 1 × d		4604.9
	RASSCF(19,19)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d	4611.3
	RASSCF(19,26)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d, 1 × f	4831.4
	RASSCF(19,35)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d, 1 × f, 1 × g	4813.7
	RASSCF(29,45)	3 × s, 2 × p, 1 × d	1 × s	3 × s, 2 × p, 1 × d, 1 × f, 1 × g	5004.8
	RASCI(29,50)	3 × s, 2 × p, 1 × d	1 × s	5 × s, 6 × p, 1 × d, 1 × f	5036.5
	Experimental				5866.9
107Ag	CASSCF(11,6)		1 × s, 1 × d		1325.1
	CASSCF(19,14)		3 × s, 2 × p, 1 × d		1455.1
	RASSCF(19,19)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d	1481.1
	RASSCF(19,26)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d, 1 × f	1594.6
	RASCI(37,54)	4 × s, 3 × p, 1 × d	1 × s	5 × s, 6 × p, 1 × d, 1 × f	1620.6
	RASCI(47,59)	4 × s, 3 × p, 2 × d	1 × s	5 × s, 6 × p, 1 × d, 1 × f	1619.4
	Experimental				1712.5
197Au	CASSCF(11,6)		1 × s, 1 × d		2362.1
	CASSCF(19,14)		3 × s, 2 × p, 1 × d		2536.0
	RASSCF(19,19)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d	2608.3
	RASSCF(19,26)	1 × s, 1 × p, 1 × d	1 × s	1 × s, 1 × p, 1 × d, 1 × f	2733.6
	RASCI(45,54)	5 × s, 4 × p, 1 × d	1 × s	4 × s, 5 × p, 1 × d, 1 × f	2704.6
	Experimental				3049.7

The RAS1, RAS2, and RAS3 columns indicate the number of atomic shells—separated by angular momentum—included in each subspace. Experimental HFCCs are reproduced from ref (75).

All RASCI calculations reported herein were performed using a RASSCF(19,26)-optimized orbital space for all coinage metals. Surprisingly, correlating most of the s and p manifold changes the initial RASSCF HFCCs very little; the largest variations are observed in calculations that correlate additional polarization shells. We observe improvements of 2–6% upon including one f shell in the orbital optimization step (cf. the RASCI results in Tables and S2). Despite the significant challenges associated with these elements, Hyperion is still able to achieve relative errors of less than 15%.

Some notable differences between the coinage metals and the alkali are highlighted by the orbital decomposition analysis. Figure shows orbital decomposition diagrams obtained from RASCI calculations on 133Cs and 197Au atoms; while both systems have a 6s1 ground configuration, their core configurations set them apart, which is reflected by the observed HFC. Compared to 133Cs, 197Au exhibits stronger coupling between d orbitals, as well as more significant contributions from the f polarization shell (however, the latter could be a consequence of including this shell in the orbital optimization step). Nevertheless, in both cases, the largest contributions to HFC are concentrated around orbitals 5s–7p, with additional non-negligible couplings involving the most diffuse s and p functions.

Groups VI-B (Cr) and VIII-B (Fe)

The hyperfine structure of transition metal (TM) atoms with partially filled d shells is by far the most challenging to model due to a number of competing factors. The orbital angular momentum couples to the nuclear spin through the PSO mechanism, and the resulting HFCC contribution is similar in magnitude to the spin-dependent FC+SD contribution (Figure S3). Both SOC and spin density must therefore be modeled accurately. On the one hand, the RASSI approach requires a sufficient number of spin-adapted states in order to represent the SOC states accurately, where the number of optimized roots corresponds to the lowest-energy Russell–Saunders (LS) terms (Table ). Additional LS terms were included for 101Ru, 183W, and 189Os to obtain a converged ordering of SO energies at the minimal CASSCF-SO level.

Table 3

Number of Spin-Adapted Roots Optimized for Each Spin S and Corresponding LS Terms

atom	number of roots	LS terms
53Cr	1 (S = 3), 6 (S = 2)	7S, 5D, 5S
95 Mo	1 (S = 3), 6 (S = 2)	7S, 5D, 5S
183W	1 (S = 3), 5 (S = 2), 3 (S = 1)	7S, 5D, 3P
57Fe	12 (S = 2)	5D, 5F
101Ru	7 (S = 2), 7 (S = 1)	5F, 3F
189Os	12 (S = 2), 10 (S = 1)	5D, 5F, 3P, 3F

On the other hand, the electronic states of TM atoms exhibit significant mixing between nd(n + 1)s2 and nd(n + 1)s1 configurations; these have competing influences on the form of the valence s orbital, as optimized singly occupied s functions are usually more radially expanded than doubly occupied s functions.[76] The dependence of the HFC operator amplifies such differences, leading to computed HFCCs that are very sensitive to variations in the CI expansion. This proves particularly problematic when CASSCF/RASSCF orbitals are averaged over states dominated by different configurations, such as the low-lying 5D (3d64s2) and 5F (3d74s1) terms of 57Fe.[82,83] Note that CASSCF(8,6)-SO HFCCs for 101Ru are in much better agreement with experiment than 57Fe HFCCs (Figure ); this is because the 101Ru 5F term is sufficiently energetically separated and can be modeled without requiring additional quintet roots in the spin-adapted CASSCF step. Similar state-averaging effects are observed in the quintet levels of 53Cr and 95Mo.

Figure 3

Unsigned HFCCs computed by Hyperion for selected energy levels of Cr group and Fe group atoms. Experimental HFCCs are reproduced from ref (80) (53Cr), (81) (95Mo, 183W), (82) and (83) (57Fe), (84) (101Ru), and (85) (189Os).

Although relatively small in magnitude, the spin-dependent part of TM HFCCs depends on SP and higher-order correlation, similar to the spin-only HFCCs of alkali and coinage metals. Unlike the 2S1/2 systems, however, the expanded RAS2 subspace precludes large RAS optimizations. Therefore, calculations that correlate the entire core region were only feasible with RASCI for 3d TMs; however, this still gave very good results (Figure ). RASCI-SO HFCCs computed for 57Fe are overall improved compared to CASSCF(8,6)-SO; meanwhile, RASCI-SO HFCCs for all 53Cr energy levels except 5S2 are slightly worse than those obtained from CASSCF-SO. We note that the RASCI-SO energies determined for 53Cr do not match the ordering observed experimentally and that the wrong ground state is predicted, indicating inaccuracies in the electronic wave function that are reflected in the calculated HFCCs. RASSCF-SO calculations correlating the valence shells were performed for all six atoms, yielding HFCCs that are overall less accurate compared to CASSCF-SO for the 3d and 5d TMs. The improvement with RASSCF observed for 95Mo and 101Ru can be justified by their well separated ground terms, which prevents state contamination during the RASSCF optimization.

Overall, it appears that minimal CASSCF-SO provides the most balanced model for the hyperfine structure of TMs, with errors around 25%. Approaches that include more correlation effects can theoretically improve the accuracy of the FC+SD term; in practice, however, the state-averaging formalism employed in CASSCF-SO, combined with the number of optimized spin-free states—which are likely characterized by different dynamical correlation effects—worsen the quality of the SO-coupled wave functions,[86,87] which propagates to the computed HFCCs.

Lanthanides

The hyperfine structure of lanthanide (Ln) atoms is dominated by the PSO term, with previous work suggesting minor contributions from core polarization.[88] As the 6s orbital and the 4f manifold are energetically well separated, CI effects between the two are negligible and the HFC response is expected to arise predominantly from the 4f shell. Hence, we compute the HFCCs of multiple levels in 4f6s2 Ln atoms using a minimal CASSCF(N,7)-SO, with the number of optimized roots (Table ) selected using a similar approach to section .

Table 4

CASSCF(N,7)-SO-Optimized LS Terms for Each Ln Atom

atom	LS terms
141Pr	4I
143Nd	5I
147Pm	6H
147Sm	7F, 5D
151Eu	8S
159Tb	6H, 4I
161Dy	5I
165Ho	4I
167Er	3H, 3F, 1G
169Tm	2F

Predicted HFCCs (Table ) are in remarkable agreement with experiment overall.[89] The poorest agreement is observed for 151Eu and 165Ho; the former has a spin-only septet ground state, which only exhibits spin-dependent HFC (no PSO contribution), and hence requires a more sophisticated treatment of correlation. Inaccuracies in the latter could also be due to missing CI effects, as indicated in previous work,[89] but increasing the active space with CASSCF or RASSCF methods had no appreciable impact on the calculated HFCCs here.

Table 5

Unsigned Isotropic HFCCs in MHz Determined for Ln Atomsa

atom	level	|AFC+SD|	|APSO|	|Atot|	|Aexpt|
141Pr	4I9/2	28.6	933.9	962.5	926.2
	4I11/2	13.9	759.4	773.4	730.4
	4I13/2	0.4	654.8	655.1	613.2
	4I15/2	13.3	587.0	573.8	541.6
143Nd	5I4	6.2	207.5	201.3	195.7
	5I5	3.8	163.1	159.3	153.7
	5I6	1.7	137.7	135.9	130.6
	5I7	0.3	121.8	122.1	117.6
	5I8	2.4	111.2	113.6	110.5
147Pm	6H7/2	42.0	472.3	430.4	447.1
147Sm	7F1	34.8	70.4	35.6	33.5
	7F2	28.6	70.4	41.8	41.2
	7F3	20.4	70.4	50.0	50.2
	7F4	11.2	70.4	59.2	59.7
	7F5	1.7	70.4	68.7	69.1
	7F6	7.7	70.4	78.1	78.4
151Eu	8S7/2	1.8	0.0	1.8	20.1
159Tb	6H15/2	25.2	720.0	745.2	673.8
	6H13/2	2.8	773.4	770.7	682.9
	6H11/2	34.0	856.8	822.9	729.0
161Dy	5I8	2.7	124.9	127.6	116.2
	5I7	0.3	136.8	137.1	126.8
	5I6	1.9	154.6	152.7	139.6
	5I5	4.3	183.2	178.9	162.0
	5I4	6.9	233.2	226.2	205.2
165Ho	4I15/2	28.6	1253.7	1225.1	800.6
	4I13/2	0.7	1398.4	1399.0	937.2
	4I11/2	29.8	1621.9	1651.7	1035.1
	4I9/2	61.5	1994.5	2056.0	1137.7
167Er	3H6	9.2	136.2	127.0	120.5
	3F4	16.8	147.1	130.4	121.8
	3H5	4.7	158.0	162.7	159.5
	3H4	5.2	168.4	173.6	173.4
	3F3	4.6	149.8	145.2	143.5
	3F2	14.4	217.9	203.5	167.1
169Tm	2F7/2	48.4	433.0	384.7	374.1
	2F5/2	115.8	577.4	693.2	704.6

Experimental HFCCs are reproduced from ref (89).

Conclusions

We have presented a new computational package, Hyperion, that enables the evaluation of relativistic picture-change-corrected magnetic resonance parameters from CASSCF-SO and RASSCF-SO wave functions, along with a new orbital decomposition method to assist in choosing appropriate active spaces for HFCC calculations. We used this code to study the hyperfine structure of alkali metal, transition metal, and lanthanide atoms in order to understand the range of applicability afforded by this approach. Our best predicted HFCCs are within 10% accuracy for alkali, 15% for coinage metals, and 20% for lanthanides (although the vast majority of Ln HFCCs deviate less than 10% from the experimental value). The hyperfine structure of group VI-B (Cr) and group VIII-B (Fe) transition metals proved to be the most challenging to model; however, despite the larger percentage errors, we obtain theoretical HFCCs that closely follow experimentally observed trends. Based on these results, we devised a number of guidelines for modeling HFC in systems exhibiting important correlation effects, strong SOC, or a combination of both. In future work, these guidelines will be refined by using Hyperion to study HFC in molecular systems.