Efficient evaluation of the Breit operator in the Pauli spinor basis

Four-component Dirac–Hartree–Fock (DHF) is the starting point for accurate all-electron relativistic calculations.^1–4 The frequency-independent instantaneous two-electron interaction provides the most accurate description^5–8 of electron–electron interaction before going to a genuine relativistic quantum electrodynamics theory of many-electron systems.^9–14 Derived from quantum electrodynamics, the form of the instantaneous two-electron interaction operator, accurate to $O(c−2)$⁠, depends on the choice of gauge. In the Feynman gauge, it can be written as

where {i, j} are electron indices. Equation (3) is also known as the Gaunt term. The components of the α matrices are defined as

and σ consists of Pauli matrices,

In the Coulomb gauge, the instantaneous two-electron operator is

where the Breit term [Eq. (7)] includes the magnetic interaction and an additional gauge term. Note that the magnetic interaction of the Breit term in the Coulomb gauge is half of the Gaunt term in the Feynman gauge.

While the Dirac–Coulomb operator accounts for the leading two-electron energy correction in the relativistic regime, the frequency-independent Breit or Gaunt term cannot be ignored for accurate description of electronic structures of late-row and heavy elements.^5–8 The difference in the two-electron interaction in Feynman and Coulomb gauge arises from different truncation schemes in the instantaneous interaction approximation, leading to gauge-dependent results. Although both forms of two-electron interactions [$VeeF$ in Eq. (1) and $VeeC$ in Eq. (6)] are accurate to $O(c−2)$⁠, there has been work showing that the Breit operator in the Coulomb gauge is more accurate than the Gaunt operator in the Feynman gauge.^7,8 It was also suggested that the variational treatment of the frequency-independent Breit interaction is important for both correlation calculations and property evaluations.^15,16 DHF with the Coulomb and Gaunt interactions has been implemented in quantum chemistry software packages for atomic and molecular systems.^17–30 However, the gauge term in the Breit interaction is computationally more expensive in both integral evaluation and integral–density contraction.

In the 1980s, frequency-dependent and -independent Breit interactions were included perturbatively in a multi-configuration Dirac-Fock formalism with finite difference methods in atomic systems.^31–33 Later on, variational treatment of the frequency-independent Breit term in atomic systems with kinetically balanced Slater type basis³⁴ and Gaussian basis^35–37 was realized computationally.

Molecular calculation with the Breit interaction came much later. Since radial-angular separation is not applicable for molecular systems, an analytical integral evaluation algorithm for the gauge term in Gaussian type orbitals is required.^38–40 The first practical gauge integral evaluation using the McMurchie–Davidson algorithm was implemented and applied for molecular systems by Quiney and co-workers.^41,42 Recently, a gauge integral with the generalized Rys-quadrature method was developed for molecules with the density fitting technique.^28,43 Still, the high computational cost of the gauge term limits its application in molecular relativistic calculations.⁴⁴ 

Some techniques have been applied to reduce the computational cost of Dirac–Hartree–Fock calculations, such as the quaternion formalism,²³ resolution of the identity,^28,45 and quasi-four component approximations.^46,47 Recently, we introduced a Pauli matrix representation with an optimal spin- and component separation algorithm that results in a minimal floating-point count algorithm for integral–density contraction.³⁰ 

In this work, we generalize the integral–density contraction in the Pauli spinor representation in the previous work³⁰ to gauge term contraction. A restricted kinetically balanced (RKB) Gaussian basis is used in this work. Component- and spin-separations⁴⁸ are applied for the Fock matrix and density matrices for maximal optimization of the contraction scheme. We also derived and developed two new methods for gauge term integral evaluation in Gaussian type orbitals: one uses the generalized Obara–Saika recursion, and the other is the second derivative of Coulomb interaction with Rys quadrature. This work includes the gauge term variationally in a relativistic four-component Dirac–Hartree–Fock framework for molecular calculations.

All expressions in this work are in atomic units. The eigenfunction of the Dirac Hamiltonian ψ has a four-spinor form,

where ϕ^L and ϕ^S are spinors of the large and small components, expanded in spinor basis,

where τ ∈ {α, β} and N is the number of spatial basis functions. The large component basis is defined as

where χ_μ is the one component spatial basis function. The relationship between the large and small component spinor basis can be defined via the restricted kinetic balance (RKB) condition^49–51 to ensure the correct nonrelativistic limit of the positive energy states.^52,53 The RKB condition can be written as

where c is the speed of light and p = −i∇ is the linear momentum operator.

Introducing the RKB condition [Eq. (11)] into the Dirac equation results in the four-component Dirac formalism in matrix representation,⁵⁴ 

where T and S are the kinetic energy and overlap matrices, respectively. The solution of Eq. (12) consists of sets of positive/negative eigenvalues ({ϵ⁺}, {ϵ⁻}) with corresponding molecular orbital coefficients $CL+CS+$ and $CL−CS−$ for the positive/negative energy solutions.

The computation of integrals of operators in Eq. (6) is the dominant cost in building the Dirac–Coulomb–Breit Hartree–Fock matrix and is usually done in the atomic-orbital (AO) basis. To build the AO Dirac–Coulomb–Breit Hartree–Fock, one can scatter the four-component density matrix into 16 component- and spin-blocked density matrices,

where τ₁, τ₂, τ₃, τ₄ ∈ {α, β}. With the RKB condition [Eq. (11)], component- and spin-blocked integrals can be constructed from four-index scalar AO integrals. Then, one can build the Fock matrix by contracting AO spinor integrals with component- and spin-blocked density matrices, resulting in the contraction scheme in the two-spinor basis [see Eqs. (A1)–(A4) in  Appendix A].

In a recent study,³⁰ we introduced a minimal floating-point-operation (FLOP) count algorithm based on the Pauli quaternion formalism of maximal component- and spin-separation. In this approach, AO density matrices are further component- and spin-separated and directly contracted with scalar integrals with a minimal FLOP approach but at the cost of an increased arithmetic complexity.

In this work, we extend the Pauli quaternion formalism for building the gauge term in the Dirac–Coulomb–Breit Hartree–Fock equation. Note that the magnetic term [first term in Eq. (7)] is one-half of the Gaunt term, which was previously derived and implemented.³⁰ We refer readers to Ref. 30 for detailed implementation of the Gaunt term the with Pauli quaternion formalism.

In order to build the gauge term contribution in the Pauli component basis, both gauge integrals and density matrices need to be expressed in the Pauli matrix representation. The Pauli components of the density matrix can be transformed from the spin-blocked matrices,⁵⁴ 

for each component-block (LL, SS, LS, SL).

The separation of gauge integrals that contract with the density matrices in the Pauli component basis utilizes the Dirac identity and orthogonality of the quaternion. Mathematical derivations are lengthy and complex, and we only present the key equations and the final working expressions [Eqs. (A17)–(A30) in  Appendix A].

Collecting all Pauli components, the two-spinor Fock matrix can be constructed:

where X, Y ∈ {L, S}. The Dirac–Hartree–Fock equation can be solved self-consistently.

In this work, we derived and implemented two methods to evaluate the gauge integral with the operator form defined as

In  Appendix B, we show the working equations to evaluate the gauge integral using the Rys-quadrature approach in terms of second derivatives of the four-center two-electron Coulomb integral. In  Appendix C, we derived the recursion relations using the generalized Obara–Saika algorithm. Both methods are implemented and benchmarked against each other, and the resulting integrals only differ by less than 10⁻¹² a.u. Integral code optimization requires extensive work, and therefore, we will not provide timing comparisons between Rys-quadrature and Obara–Saika algorithms for the gauge term.

In this section, we focus on analyzing the memory requirement and computational cost for building the Dirac–Fock matrix from four-index electron-repulsion-integrals (ERIs). We use the in-core formalism where all four-index integrals are computed and stored in memory for the theoretical cost analysis. In this analysis, we used a Kramers-unrestricted condition where all four Pauli components are non-zero. The memory requirement and FLOP count for Dirac–Coulomb and Gaunt have been previously analyzed in Ref. 30. Note that the magnetic contribution in the Breit term has the same computational cost as the Gaunt term. Therefore, the previous cost analyses were directly used in this work with the addition of the gauge contribution.

The gauge term is built from 16 four-index scalar (or one-component) integrals shown in Eqs. (A11)–(A14). In the two-spinor RKB approach, these 16 types of gauge ERIs are assembled into 64 two-spinor ERIs, and the integral–density contraction is carried out in the two-spinor form [see Eqs. (A1)–(A4)]. In contrast, the Pauli component approach does not assemble the 16 gauge ERIs into a large dimension, although several auxiliary scalar integrals [e.g., (σ · σ), (σ × σ), and (σσ)_xz + (σσ)_zx] can be constructed to reduce the cost of contraction. Instead, we scatter the Dirac–Fock and density matrices into the Pauli representation so that all gauge ERIs remain in the scalar form and the integral–density contraction is done in Pauli quaternion representation (see Sec. II B). Tables I and II list the estimated computational costs in terms of memory requirement and FLOP count for forming the Dirac–Fock matrix from integrals.

Table I shows that the Pauli matrix formalism has a much smaller memory storage requirement compared to the two-spinor matrix form. This is mainly because in the Pauli matrix basis, all integrals are stored in one-component real-valued scalar form. In contrast, in the two-spinor basis, scalar integrals are assembled into two-component complex-valued tensors. As a result, the in-core integral storage requirement is much bigger in the two-spinor basis than in the Pauli matrix basis.

Table II gives the FLOP count for each type of Coulomb J and exchange K contractions with the density. For both the Dirac–Coulomb and the Breit terms, we see that the Pauli matrix formalism requires a much smaller number of J contractions than the two-spinor approach. However, the number of K matrix formations is much bigger in the Pauli basis than in the two-spinor basis. Since all integrals and density matrices in the two-spinor basis are complex-valued, each integral–density contraction requires 8N⁴ FLOPs. In contrast, each integral–density contraction needs 4N⁴ FLOPs in the Pauli basis because integrals are real-valued. As a result, the computational cost in terms of FLOP count is much more efficient in the Pauli basis.

The purpose of the data presented here is to display the magnitude of the effect of the Dirac–Coulomb–Breit Hamiltonian in atomic systems. In particular, DHF energies are compared to one another in this section. Integral computation is performed with the LIBCINT integral library.⁵⁵ All calculations in this section are performed with a development version of the Chronus Quantum software package.⁵⁶ The speed of light utilized in this section is 137.035 999 679 94 a.u.

Ground state DHF energies were computed utilizing the Dirac–Coulomb and Dirac–Coulomb–Breit Hamiltonians. All atoms are defined to be neutral in these calculations. A core Hamiltonian guess is utilized, such that no spin multiplicity need to be defined for the wavefunction guess. The self-consistent-field optimization proceeds in the Kramers unrestricted framework in which the lowest N_e number of positive energy orbitals is singly occupied. In order to analyze the importance of relativistic two-electron corrections in Dirac Hartree–Fock, we also compare the results obtained using the (LL|LL) approximation (previously referred to as the Bare–Coulomb³⁰ approximation), i.e., only (LL|LL) integrals are used without contributions from the small component. The four-component Hamiltonians presented here differ only in their treatment of the two-electron interaction. This allows for the isolation of the effect of various two-electron corrections, particularly the terms of the Dirac–Coulomb–Breit Hamiltonian presented here.

The ground state energies of selected atomic species from across the Periodic Table were computed and relativistic two-electron contributions are compared in Table III. To obtain the ground state energies, we have carried out separate self-consistent-field (SCF) calculations with each of the three Hamiltonians and the (LL|LL) approximation. Relativistic two-electron corrections are defined as the differences between different SCF energies (see the caption of Table III). The uncontracted ANO-RCC basis^58–59 set was utilized for all cases. The results are listed in Table III and plotted in Fig. 1. As the atomic number Z increases, the Dirac–Coulomb (ΔE^DC) and Breit (ΔE^Breit) corrections increase as expected, as shown in Fig. 1(a). Compared to ΔE^DC, the Breit correction becomes less important as the atomic number increases, as indicated by the decreasing ratio of $ΔEBreitΔEDC$ [see Fig. 1(b)].

We plot in Fig. 1(c) the magnetic and gauge contributions (ΔE^Mag and ΔE^Gauge) in the Dirac–Coulomb–Breit Hamiltonian and their ratio $ΔEGaugeΔEMag$ in Fig. 1(d). Both the magnetic and gauge terms grow exponentially with respect to the atomic number and both contributions are equally important in the Breit Hamiltonian. Figure 1(d) shows that the ratio $ΔEGaugeΔEMag$ approaches unity toward light elements and ∼0.78 toward heavy elements. Recall that the Breit operator is the Gaunt term in the Feynman gauge, and in the Coulomb gauge, consists of half of the Gaunt term (i.e., magnetic contribution) and the gauge term. Therefore, another way to understand Fig. 1(d) is that the energetic difference between the Feynman gauge and the Coulomb gauge is approaching ∼0% for light elements and ∼20% for heavy elements.

In this section, we analyze the contribution of the Breit Hamiltonian to the ground state energy and its effect on the electronic structure characteristics of selected heavy-element compounds from the actinyl oxycation series. Linear geometries and experimental bond lengths of 1.76, 1.75, and 1.74 Å were utilized for uranyl, neptunyl, and plutonyl, respectively.⁶⁰ 

Table IV gives the energetic decomposition of different relativistic two-electron interaction terms in the total energy of uranyl, neptunyl, and plutonyl. The uncontracted ANO-RCC-VDZP basis set was used for both the oxygen and actinide atoms. The main trends seen across the actinyl series reflect the growing nuclear charge of the actinide atom and mirror the atomic periodic trends seen in Sec. III B. Although the relative contribution of the Breit term is small (⁠$∼0.1%$ of the total energy), its effect on electronic structure characteristics of the molecular bonding could be significant.

Density of states plots were generated using the occupied Dirac–Coulomb and Dirac–Coulomb–Breit molecular orbitals for the three actinyls presented in Table IV. These plots are presented in Fig. 2 and show a total large component density as well as a large component contribution from each angular momentum orbital type, i.e., s, p, d, and f characters. Orbital character is summed over the three atoms in the molecule.

Figures 2(a) and 2(b) show that the highest occupied orbitals of uranyl consist of six 5f-2p and six 6d-2p bonding orbitals with a very small contribution of s character. The occupied nonbonding 5f electron molecular orbitals in neptunyl and plutonyl are the next highest energy orbitals, appearing lower in energy than the 12 bonding orbitals.

Comparing the DOS plots for this series, the energies of bonding molecular orbitals in neptunyl are more spread out than for those of uranyl and plutonyl. Comparing the DOS plots computed using the Dirac–Coulomb and Dirac–Coulomb–Breit Hamiltonians, we see that the Breit term does not significantly alter the character of the orbitals or the general ordering of the molecular orbitals presented here, but it does alter their relative energetics. For example, a distinct sharpening of the highest energy peak in each of these systems can be observed in the Dirac–Coulomb–Breit case. This is a result of narrowing of the frontier orbital energy gap near HOMO. However, for deeper valence orbitals, there is no noticeable difference in the DOS plots computed using Dirac–Coulomb and Dirac–Coulomb–Breit Hamiltonians.

In order to investigate the effect of the Breit operator on energetic observables, the dissociation energies of uranyl, neptunyl, and plutonyl into 2 × O²⁻ and U(VI), Np(VI), and Pu(VI), respectively, were calculated with both Dirac–Coulomb and Dirac–Coulomb–Breit. To compute the SCF energies of atomic An(VI) cations, a core Hamiltonian guess is utilized, such that no spin multiplicity need to be defined for the wavefunction guess. The SCF optimization proceeds in the Kramers unrestricted framework in which the lowest N_e number of positive energy orbitals is singly occupied. These dissociation energies are given in Table V. Note that this study is meant to analyze the importance of the Breit contribution to bond dissociation energy at the Dirac–Hartree–Fock level. Accurate bond dissociation energies will require calculations using electron correlation methods.

Dissociation energies, and therefore bonding energies, can be seen to be higher when calculated with Dirac–Coulomb–Breit than with the Dirac–Coulomb Hamiltonian in these compounds, suggesting that the Breit Hamiltonian stabilizes the chemical bonds between actinide and oxygen. This can be understood from the DOS plots in Fig. 2. The Breit Hamiltonian narrows the energy gap among frontier orbitals near HOMO. As a result, orbitals lower than HOMO will contribute more to the molecular bonding with Dirac–Coulomb–Breit than with the Dirac–Coulomb Hamiltonian. The difference between the calculated dissociation energies grows larger with an increasing atomic number of the actinide atom.

In this work, we applied Pauli spinor representation in the restricted–kinetic–balance condition for four-component Dirac–Hartree–Fock calculations with the Coulomb–Breit two-electron interaction. Along with the density–integral contraction formalism, we introduced two different methods for calculating frequency-independent gauge term integrals with the Gaussian basis. One method utilizes the generalized Obara–Saika algorithm, and the other is the second derivative of Coulomb interaction with Rys quadrature.

With this maximally component- and spin-separated formalism, the computational cost of density–integral contraction for gauge term in the Breit interaction is minimized. It is shown that this Pauli spinor representation is equivalent to the two-spinor formalism but with a much lower computational cost.

We applied this method for calculating atomic energies in order to analyze the periodic trend of the Breit term. It was shown that the Breit contribution becomes relatively less important compared to the Dirac–Coulomb term as the atomic number increases. The analysis also suggests that the difference between the Breit terms in the Feynman gauge and the Coulomb gauge approaches ∼20% toward the heavy elements.

We also studied the effect of the Breit term on the electronic structure characteristics of actinyl oxycations. The analysis shows that the Breit term modifies the relative energies of frontier molecular orbitals near HOMO. As a result, it leads to a slightly stronger chemical bond.

Future development will be focused on combining this Pauli spinor formalism with atomic balance^52,61 and minimal-FLOP Cholesky integral decomposition⁶² to further reduce the computational cost of four-component Dirac–Hartree–Fock.

X.L. acknowledges support from the U.S. Department of Energy, Office of Science, Basic Energy Sciences, in the Heavy-Element Chemistry program (Grant No. DE-SC0021100) for the development of electronic structure methods within the four-component Dirac framework.

The authors have no conflicts to disclose.

Shichao Sun: Data curation (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (equal); Writing – review & editing (equal). Jordan Ehrman: Data curation (equal); Writing – original draft (equal); Writing – review & editing (equal). Qiming Sun: Conceptualization (supporting); Software (supporting); Writing – original draft (supporting); Writing – review & editing (supporting). Xiaosong Li: Conceptualization (equal); Software (equal); Writing – original draft (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The gauge term contribution to the Fock matrix in the two-spinor basis reads

where subscript “3” denotes that the integral is $1r123$ in contrast to $1r12$ and superscript g denotes the gauge term contribution.

In this work, we perform component- and spin-separation for gauge integrals. We use the Dirac identity,¹ $(σ⋅μ)(σ⋅ν)=Iμ⋅ν+iσ⋅μ×ν$⁠, and introduce ξ_τ as the spin representation,

We have the following identities:

otherwise zero.

The gauge integral in Eq. (A1) can be written as

Similarly, the integral in the first term of Eq. (A3) is

The integral in the second term of Eq. (A3) is

and the exchange integral of Eq. (A3) is

and the exchange integral of Eq. (A2),

where ∇_κ is the nuclear coordinate derivative of the atomic basis χ_κ. I(2) and σ(2) are the identity and Pauli matrix for electron 2.

For simplicity and brevity, we introduce a simplified notation for two-spinor integrals,

The “^-” over ν, κ means derivatives are taken of these orbitals. To expand the spinor basis integral $(μτ1Lσντ2S|κτ3Sσλτ4L)3$ as linear combination of second derivatives of spatial basis integrals as shown in Eq. (A10), we introduce the following short notations for the second derivative of spatial basis gauge integrals:

where r₁₂r₁₂ in Eqs. (A11)–(A14) is the operator evaluated in the integral. We move it outside of the ()₃ for brevity. We give the definition for integrals with derivatives on the second and third centers.

Other types of integrals where the derivatives are taken on the second and fourth centers, or the first and fourth centers, can be derived from Eqs. (A11)–(A14) using symmetry. For the V^g,LL, V^g,LS, and V^g,SS terms, we need gauge spinor integrals of the type (LS|SL), (SL|LS), (LS|LS). By using the symmetry of the spinor integrals, we can express the gauge integrals of the type (SL|LS), (LS|LS) by Pauli components of (LS|SL). Comparing Eq. (A9) with Eq. (A6), we have

Similarly, comparing Eq. (A7) with Eq. (A6), we have

The discussion above illustrates the mathematical expressions and notations needed to derive the gauge contribution to the Fock matrix in the component- and spin-separated form. We skip the detailed derivations and only present final working expressions. Note that the gauge contribution to the Fock matrix has a prefactor of $1(2c)2$⁠. For brevity, we omit the prefactor in all equations.

The gauge term contribution to the Pauli components of the LL block of the Fock matrix is an exchange type,