The fully correlated frequency-independent Dirac–Coulomb–Breit Hamiltonian provides the most accurate description of electron–electron interaction before going to a genuine relativistic quantum electrodynamics theory of many-electron systems. In this work, we introduce a correlated Dirac–Coulomb–Breit multiconfigurational self-consistent-field method within the frameworks of complete active space and density matrix renormalization group. In this approach, the Dirac–Coulomb–Breit Hamiltonian is included variationally in both the mean-field and correlated electron treatment. We also analyze the importance of the Breit operator in electron correlation and the rotation between the positive- and negative-orbital space in the no-virtual-pair approximation. Atomic fine-structure splittings and lanthanide contraction in diatomic fluorides are used as benchmark studies to understand the contribution from the Breit correlation.

## I. INTRODUCTION

In many-body electronic structure theory, the frequency-independent Dirac–Coulomb–Breit operator defined in the Coulomb gauge provides the most accurate description^{1–5} of electron–electron interaction before going to a genuine relativistic quantum electrodynamics theory.^{6–11} The Dirac–Coulomb–Breit operator can be written as

where {*i*, *j*} are electron indices. The components of the ** α** matrices are defined as

and ** σ** consists of Pauli matrices,

The Breit operator [Eq. (3)] includes the magnetic interaction and a gauge term, which is the most computationally intensive. A lower cost operator that still carries the magnetic interaction is the Dirac–Coulomb–Gaunt Hamiltonian defined in Feynman gauge,

It is generally well understood that the Breit operator in the Coulomb gauge is more accurate than the Gaunt operator in the Feynman gauge^{3,4} and that the variational treatment of the frequency-independent Breit interaction is important for both correlation calculations and property evaluations.^{12,13}

Although the Dirac–Coulomb and Dirac–Coulomb–Gaunt Hamiltonians have been used in correlated many-body calculations,^{14–16} due to the large computational cost of the gauge term in the Breit operator, the application of the Dirac–Coulomb–Breit Hamiltonian mostly remains in the mean-field theory or used as perturbative treatment in the correlated methods. Frequency-independent Breit interactions have been implemented with finite difference methods,^{17–19} Slater type basis,^{20} and Gaussian basis.^{5,21–29} Recently, we introduced a Pauli matrix quaternion representation with an optimal spin- and component-separation algorithm that results in a minimal floating-point count algorithm for building the Dirac–Coulomb–Breit Hamiltonian.^{5,29} The reduced computational cost of the four-component density-integral contraction allows for the development of fully correlated Dirac–Coulomb–Breit many-body methods.

In this work, we introduce the fully correlated Dirac–Coulomb–Breit many-body theory in the variational four-component multiconfigurational self-consistent-field (MCSCF) framework, including the complete active space self-consistent field (CASSCF) and density-matrix renormalization group SCF (DMRGSCF) methods. To the best of our knowledge, this is the first time the Dirac–Coulomb–Breit operator has been applied in the context of DMRG. We will place a special emphasis on the importance of positive-negative-energy orbital rotation in the MCSCF procedure and the contribution of the Breit correlation beyond the Dirac–Coulomb operator.

## II. METHOD

In this work, we develop correlated Dirac–Coulomb–Breit four-component multiconfigurational self-consistent-field methods (DCB-MCSCF) in the frameworks of the complete active space self-consistent field (CASSCF) and density matrix renormalization group self-consistent field (DMRGSCF). Since the fundamental expressions of four-component CASSCF and DMGRSCF are similar to exact-two-component (X2C) implementations, we recommend the readers to Refs. 30 and 31 for more information on Kramers-unrestricted CASSCF and DMRG, respectively. In this section, we focus on methodological developments that are unique to the DCB-MCSCF theory.

In this section, the following notation will be used: *i*, *j*, *k*, *l* label inactive four-spinors; *t*, *u*, *v*, *w* label active four-spinors; *a*, *b*, *c*, *d* label virtual four-spinors; *p*, *q*, *r*, *s* label general four-spinors; and *I*, *J*, *K* label Slater determinants. Indices with positive and negative superscript (e.g., *i*^{+}, *j*^{+} and *a*^{−}, *b*^{−}) refer to positive-energy and negative-energy four-spinors, respectively. The following discussions assume that all quantities are complex-valued unless stated otherwise.

### A. Four-component multiconfigurational wave function in no-virtual-pair approximation

The Dirac equation is cast in a finite Gaussian type basis. The four-spinor molecular orbitals (MO) are expanded in two-spinor basis as

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

where *χ*_{μ} are spatial basis functions. The relationship between the large and small component two-spinor basis can be defined via the restricted kinetic balance (RKB) condition, which ensures the correct nonrelativistic limit of the positive energy states,^{32–37}

In the RKB condition, the Dirac–Hartree–Fock Hamiltonian in matrix form can be efficiently built and solved in the Pauli matrix quaternion representation (see Refs. 29 and 5 for Dirac–Coulomb–Gaunt and Dirac–Coulomb–Breit Dirac–Hartree–Fock, respectively).

The solution of the four-component Dirac–Hartree–Fock equation consists of sets of positive and negative eigenvalues {*ϵ*^{+}}, {*ϵ*^{−}} with corresponding molecular orbital coefficients $(CL+CS+)T$ for the positive and $(CL\u2212CS\u2212)T$ for the negative energy solutions. These are mean-field spin-coupled, relativistically corrected molecular four-spinors that will be used for subsequent correlated calculations. We employ the conventional no-virtual-pair approximation (NVPA) where only positive-energy orbitals $(CL+CS+)T$ are considered in the multiconfigurational expansion.^{38–40}

#### 1. DCB-CASSCF

For the NVPA CAS wave function, |Ψ^{+}⟩ is described as a linear combination or a configuration interaction (CI) expansion of Slater determinants, |*K*^{+}⟩, constructed from a subset of the orthonormal positive-energy four-spinors,

where *N*_{det} is the total number of determinants in the expansion.

The correlated DCB-CASSCF energy is written as

where *ϵ*^{c} is the Dirac–Hartree–Fock (e.g., using the Dirac–Coulomb–Breit Hamiltonian) energy of the inactive orbitals (the “core energy”) and is defined as

where $VeeDCB$ is the Dirac–Coulomb–Breit operator in Eq. (1). $H\u0302CAS$ is the active space Hamiltonian and is defined as

where *h*^{c} is the core Dirac–Fock matrix and is given by

and $E\u0302p+q+$ is an operator that excites from four-spinor orbital *q*^{+} to *p*^{+} and can be written using the usual creation and annihilation operators $E\u0302p+q+=ap+\u2020aq+$. The excitation list is generated using Handy’s string based approach^{41} adapted to the four-component orbital space.

In NVPA DCB-CASSCF, only positive-energy four-spinors are used in the CI expansion and energy evaluation, and all corresponding two-electron integrals in both mean-field and correlation treatments are evaluated with the Dirac–Coulomb–Breit operator.

#### 2. DCB-DMRGSCF

The DMRG wave function can be written as

where $CI1I2\u2026IN$ are complex-valued configuration-interaction coefficients, *N* is the number of active four-spinors, and $I$ is an occupation number vector comprised of Kramers-unrestricted four-spinor orbitals. Upon refactoring to the matrix-product state (MPS) ansatz, the no-pair wave function used in DMRGSCF is

where {*I*_{p}} are occupations and {*a*_{i}} are the virtual indices (bond dimension) that are capped by the parameter *m*_{max} = 2^{N/2}, commonly referred to as *number of renormalized block states*.^{42} Variational relativistic DMRG methods are an emerging area of research.^{31,43–45}

We enforce the NVPA on the Kramers-unrestricted MPS wave function by placing constraints on occupations [*I*_{p} in Eq. (16)]. We restrict negative-energy four-spinors (*p*^{−}) to be unoccupied (vac) but allow positive-energy four-spinors (*p*^{+}) to be unoccupied or occupied (occ),^{43,44}

### B. Positive-negative-energy orbital rotation

The four-component multiconfigurational wave function is optimized by making the energy stationary with respect to variations of *both* the configuration parameters and the four-spinor coefficients. For the complex configuration interaction (CI) coefficients, this is achieved via solution of the CASCI eigenvalue equation.

The complex four-spinor orbitals are optimized via a unitary transformation of the existing orbitals,

where the **R** matrix may be written in terms of an anti-Hermitian orbital-rotation matrix **X**,

The CASSCF wavefunction is optimized when

There are three orbital optimization strategies investigated in this work:

No orbital rotation: In this scheme, the four-component multiconfigurational methods become DCB-CASCI and DCB-DMRG using positive- and negative-energy four-spinors from the solution of Dirac–Coulomb–Breit Hartree–Fock without further orbital optimization.

Rotation within positive-energy orbitals only: In this scheme, the orbital rotation matrix [Eq. (21)] is restricted to positive-energy space only, i.e., $Xp+q+$. The negative-energy four-spinors are from the solution of Dirac–Coulomb–Breit Hartree–Fock without further orbital optimization.

Rotation in the full space including positive- and negative-energy orbitals: In this scheme, the orbital rotation matrix [Eq. (21)] also considers the rotation between positive- and negative-orbital space in CASSCF and DMRGSCF energy minimization, i.e., $Xp+q\xb1$.

For Dirac–Hartree–Fock, we use rotation in the full space including both positive- and negative-energy orbitals and enforce NVPA by restricting the density construction to only positive-energy spinors (our single determinant has no negative-energy spinor occupation).

### C. AO-direct integral transformation of Dirac-Coulomb–Breit Hamiltonian

The correlated Dirac–Coulomb–Breit many-body theory requires two-electron repulsion integrals (ERIs) in the form of $(pq|VeeDCB|rs)$, where {*p*, *q*, *r*, *s*} are general four-component orbitals and $VeeDCB$ is defined in Eq. (1). In correlated relativistic treatment, the integral transformation from the atomic orbital to four-spinor basis is the most computational expensive step. Particularly, in the correlated DCB Hamiltonian, the gauge integral transformation is the most complicated and computational dominant procedure.

The contraction scheme in the Pauli quaternion representation^{5,29} recently developed by the team allows for fast transformation directly, without going through an in-core assembly of spinor integrals, from quantities computed in the atomic orbital basis to ERIs in four-spinor basis. For example, the AO-direct transformation of the gauge term starts from four basic types of scalar gauge integrals in the atomic orbital basis (see Ref. 5 for detailed derivations),

where the subscript “3” denotes that the integral is $1r123$ in contrast to $1r12$. The contraction between scalar AO integrals with density [$O(N6)$ method] or four-spinor coefficients [$O(N5)$ method] can be easily carried out in the Pauli quaternion representation. This approach results in an efficient transformation procedure that leads to the AO-direct-enabled correlated Dirac–Coulomb–Breit many-body theory. Although the technical advances are non-trivial to implement, the theoretical foundation is well illustrated in Refs. 5 and 29.

## III. RESULTS AND DISCUSSION

The Dirac–Coulomb–Breit-CASSCF calculations were performed in the Chronus Quantum^{46} software package. All calculations in this work are all electron and use a Kramers-unrestricted spinor basis in the no-virtual-pair approximation. The Dirac–Coulomb–Breit-DMRGSCF implementation is built off our previous modular interface between QCMaquis^{47–49} and Chronus Quantum^{46} computational packages.^{31} The integrals for the four-component Hamiltonians and orbital optimizations were computed in Chronus Quantum, and the matrix-product state optimizations were performed in QCMaquis. For more information on four-component integral evaluation, see Refs. 5 and 29. We use a “state-specific” approach for optimization of the matrix-product-state wave function^{50} and a “state-averaged” approach for orbital optimization.^{51} The matrix product state utilized no additional reordering of orbitals. Unlike CASSCF, the active–active orbital rotations are non-redundant for DMRGSCF (with unconverged *m*); however, we did not include them in this study.^{52,53} Orbital optimizations were performed with an approximated quasi-second-order Newton–Raphson method.^{54}

### A. Importance of positive-negative-energy orbital rotation

#### 1. Core electrons

In the no-virtual-pair approximation (NVPA), only positive energy molecular orbitals are used in the correlated wave function treatment, such as the configuration interaction expansion. However, positive- and negative-energy orbitals are not completely disentangled from each other due to the electron correlation. Therefore, the orbital rotation in the correlated relativistic method should be carefully examined. In this section, three strategies are examined within NVPA: no orbital rotation (e.g., CASCI), orbital rotation within the positive energy space (denoted as CASSCF^{+}), and orbital rotation in the full four-component space (denoted as CASSCF^{±}). In the full orbital rotation, including both positive and negative energy spaces, the negative energy orbitals are treated as virtuals.

We first examine the importance of orbital rotation for core-electrons. The ground state energy of Rn^{84+} was computed using three Hamiltonians, including Dirac–Coulomb (DC), Dirac–Coulomb–Gaunt (DCG), and Dirac–Coulomb–Breit (DCB), in the correlated CAS(2,10) treatment. The CAS(2,10) consisted of the *ns*, (*n* + 1)*s*, and (*n* + 1)*p*. Uncontracted ANO-RCC-MB (255 basis functions including up to *f* function) is used in this study. In this basis set, the tightest atomic *s*-orbital has an exponent of 53 906 998. Table I compares the correlation energy computed using the three different orbital rotation strategies in correlated CAS calculations. The correlation energy, *E*_{c}, is defined as the energy difference between four-component correlated calculation and the four-component Hartree–Fock mean-field theory with the same relativistic operator.

. | E_{c} (eV)
. |
---|---|

E(DCB-CASSCF^{±})−E(DCB-HF) | −3.099 |

E(DCB-CASSCF^{+})−E(DCB-HF) | −3.404 |

E(DCB-CASCI)−E(DCB-HF) | −0.155 |

E(DCG-CASSCF^{±})−E(DCG-HF) | −3.804 |

E(DCG-CASSCF^{+})−E(DCG-HF) | −4.237 |

E(DCG-CASCI)−E(DCG-HF) | −0.181 |

E(DC-CASSCF^{±})−E(DC-HF) | −1.141 |

E(DC-CASSCF^{+})−E(DC-HF) | −1.225 |

E(DC-CASCI)−E(DC-HF) | −0.059 |

. | E_{c} (eV)
. |
---|---|

E(DCB-CASSCF^{±})−E(DCB-HF) | −3.099 |

E(DCB-CASSCF^{+})−E(DCB-HF) | −3.404 |

E(DCB-CASCI)−E(DCB-HF) | −0.155 |

E(DCG-CASSCF^{±})−E(DCG-HF) | −3.804 |

E(DCG-CASSCF^{+})−E(DCG-HF) | −4.237 |

E(DCG-CASCI)−E(DCG-HF) | −0.181 |

E(DC-CASSCF^{±})−E(DC-HF) | −1.141 |

E(DC-CASSCF^{+})−E(DC-HF) | −1.225 |

E(DC-CASCI)−E(DC-HF) | −0.059 |

Table I suggests that the orbital rotation plays a very important role in recovering the correlation energy, exemplified by the large difference between CASCI and CASSCF^{+}/CASSCF^{±} calculations. When comparing the partial orbital rotation within the positive energy space in CASSCF^{+} and the full rotation in CASSCF^{±}, a ∼10% over-estimation of the correlation energy was observed for DCB-CASSCF^{+}, DCG-CASSCF^{+}, and DC-CASSCF^{+}. This analysis suggests that orbital rotation in the full four-component space including both positive and negative energy orbitals is essential in accurate correlated treatment of core electrons and cannot be neglected at the correlated level.^{15}

#### 2. Valence electrons

We now focus on the importance of positive-negative orbital rotation on valence electrons. We report in Table II correlation energy of CASCI, CASSCF^{+}, and CASSCF^{±} for alkaline earth metals of Be to Ra with an active space consisting of the outermost *ns* electrons and the *ns*, *np*, (*n* + 1)*s*, and (*n* + 1)*p* orbitals [CAS(2,16)] for Be–Ca and the *ns*, *np*, (*n* − 1)*d*, (*n* + 1)*s*, (*n* + 1)*p*, and *nd* orbitals for Sr–Ra [CAS(2,36)]. It is obvious that CASCI with no orbital optimization, compared to CASSCF (both “^{+}” and “^{±}”), misses a large part of the electron correlation. In contrast to the case for core electrons, doing the full positive-negative orbital rotation in CASSCF^{±} for valence electrons does not significantly change the correlation energy compared to CASSCF^{+}. This observation suggests that positive-negative orbital rotation with CASSCF becomes less important for valence electrons further away from the nucleus and can be treated at the Dirac–Hartree–Fock level. Comparing correlation energies computed using three different relativistic Hamiltonians will be discussed in Sec. III B.

. | Be . | Mg . | Ca . | Sr . | Ba . | Ra . |
---|---|---|---|---|---|---|

E(DCB-CASSCF^{±})−E(DCB-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.560 |

E(DCB-CASSCF^{+})−E(DCB-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.560 |

E(DCB-CASCI)−E(DCB-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

E(DCG-CASSCF^{±})−E(DCG-HF) | −1.194 | −0.854 | −0.748 | −0.683 | −0.635 | −0.560 |

E(DCG-CASSCF^{+})−E(DCG-HF) | −1.194 | −0.854 | −0.748 | −0.683 | −0.635 | −0.560 |

E(DCG-CASCI)−E(DCG-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

E(DC-CASSCF^{±})−E(DC-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.561 |

E(DC-CASSCF^{+})−E(DC-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.561 |

E(DC-CASCI)−E(DC-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

. | Be . | Mg . | Ca . | Sr . | Ba . | Ra . |
---|---|---|---|---|---|---|

E(DCB-CASSCF^{±})−E(DCB-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.560 |

E(DCB-CASSCF^{+})−E(DCB-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.560 |

E(DCB-CASCI)−E(DCB-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

E(DCG-CASSCF^{±})−E(DCG-HF) | −1.194 | −0.854 | −0.748 | −0.683 | −0.635 | −0.560 |

E(DCG-CASSCF^{+})−E(DCG-HF) | −1.194 | −0.854 | −0.748 | −0.683 | −0.635 | −0.560 |

E(DCG-CASCI)−E(DCG-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

E(DC-CASSCF^{±})−E(DC-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.561 |

E(DC-CASSCF^{+})−E(DC-HF) | −1.192 | −0.854 | −0.748 | −0.683 | −0.636 | −0.561 |

E(DC-CASCI)−E(DC-HF) | −0.319 | −0.332 | −0.665 | −0.559 | −0.458 | −0.262 |

### B. Importance of Breit operator in correlation energy

Tables I and II also compare the correlation energy computed using the three different relativistic Hamiltonians in correlated calculations. Figure 1 shows the trend of CAS(2,10) correlation energy in He-like two-electron systems as the nuclear charge increases. Larger active spaces, including CAS(2,18) and CAS(2,46), were also tested (see the supplementary material). For CASSCF^{+} and CASSCF^{±}, large active spaces do not change the results, suggesting that the benchmark tests shown in Fig. 1 have reached the full MCSCF limit. This is not the case for CASCI calculations, which need a much larger active space to reach the full CI limit.

From Table I, it is obvious that the correlation energy for Rn^{84+} arising from the Breit operator [−3.099 − (−1.141) = −1.958 eV] is much bigger than that from the Dirac–Coulomb term (−1.141 eV). This observation suggests that the Breit correlation becomes the dominant contribution for deep core electrons.^{7} Compared to the Dirac–Coulomb–Breit result (DCB-CASSCF^{±}), the Dirac–Coulomb–Gaunt Hamiltonian over-estimates the correlation energy by ∼25%, whereas the Dirac–Coulomb term underestimates it by ∼63%. These behaviors can be understood from the nature of the corresponding Hamiltonians. The Dirac–Coulomb Hamiltonian lacks the magnetic and gauge interactions, which play a significant role in the electron–electron interaction of deep core electrons. The Gaunt term is twice the magnetic interaction in the Breit Hamiltonian, which has an additional gauge term. As discussed in Ref. 5, the difference between the Gaunt and Breit Hamiltonians results in ∼20% difference in relativistic mean-field energy correction.^{5} Table I suggests that this difference is also manifested in the correlation energy.

Figure 1 shows that while the difference in correlation energy between DC, DCG, and DCB is relatively small for core electrons in light elements, it quickly becomes significant as the nuclear charge increases (see the supplementary material for all computed results, including CASCI). For Ne^{8+}, the DC-CASSCF^{±} and DCG-CASSCF^{±} differ by 31 meV, and DCG-CASSCF^{±} and DCB-CASSCF^{±} correlation energies differ by 7 meV. They quickly increase to 101 meV between DC-CASSCF^{±} and DCG-CASSCF^{±} and 22 meV between DCG-CASSCF^{±} and DCB-CASSCF^{±} for Ar^{16+}— noticeable differences among three different relativistic operators. In addition, the difference between CASSCF^{±} and CASSCF^{+} orbital rotation schemes also becomes noticeable starting at Ar^{16+}.

In contrast, for valence electrons, we find no strong dependence of electron correlation on the Breit operator, as shown in Table II. In addition, the CASSCF (both “^{+}” and “^{±}”) correlation energy for the valence electrons decreases as atomic number increases, exhibiting an opposite trend compared to that for the core electrons. The CASCI results do not display the same trend as CASSCF due to a lack of optimization of the virtual orbitals.

### C. Atomic fine-structure splitting

In Table III, we report atomic fine-structure splitting for Ga, In, Tl, La, and Te using state-average Dirac–Coulomb, Dirac–Coulomb–Gaunt, and Dirac–Coulomb–Breit 4C-CASSCF^{±} and 4C-DMRGSCF^{±}. For the atomic systems, an uncontracted ANO-RCC-VTZP basis set was used.^{56,57} The wave functions used a valence active space and averaged over the total number of *m*_{J} states involved in the fine-structure splitting. For Ga, In, and Tl, we used a SA(6)-CAS(3,8), for La, we used a SA(10)-CAS(3,12), and for Te, we used a SA(9)-CAS(6,8). SA(*N*) denotes a CASSCF obital optimization averaged over *N* number of states. Convergence thresholds of 10^{−7} au in energy and 5 × 10^{−4} au in orbital rotation gradient were used for 4C-CASSCF^{±} and 4C-DMRGSCF^{±}.

. | . | 4C-CASSCF^{±}
. | 4C-DMRGSCF^{±}
. | Expt.^{55}
. |
---|---|---|---|---|

Ga $(P1/22\u2192P3/22)$ | DC | 0.087 | 0.087 | 0.102 |

DCG | 0.085 | 0.085 | ||

DCB | 0.085 | 0.085 | ||

In $(P1/22\u2192P3/22)$ | DC | 0.233 | 0.233 | 0.274 |

DCG | 0.230 | 0.230 | ||

DCB | 0.230 | 0.230 | ||

Tl $(P1/22\u2192P3/22)$ | DC | 0.849 | 0.849 | 0.966 |

DCG | 0.839 | 0.839 | ||

DCB | 0.839 | 0.839 | ||

La $(D3/22\u2192D5/22)$ | DC | 0.086 | 0.086 | 0.131 |

DCG | 0.083 | 0.083 | ||

DCB | 0.083 | 0.083 | ||

Te $(P23\u2192P03)$ | DC | 0.589 | 0.589 | 0.584 |

DCG | 0.584 | 0.584 | ||

DCB | 0.584 | 0.584 |

. | . | 4C-CASSCF^{±}
. | 4C-DMRGSCF^{±}
. | Expt.^{55}
. |
---|---|---|---|---|

Ga $(P1/22\u2192P3/22)$ | DC | 0.087 | 0.087 | 0.102 |

DCG | 0.085 | 0.085 | ||

DCB | 0.085 | 0.085 | ||

In $(P1/22\u2192P3/22)$ | DC | 0.233 | 0.233 | 0.274 |

DCG | 0.230 | 0.230 | ||

DCB | 0.230 | 0.230 | ||

Tl $(P1/22\u2192P3/22)$ | DC | 0.849 | 0.849 | 0.966 |

DCG | 0.839 | 0.839 | ||

DCB | 0.839 | 0.839 | ||

La $(D3/22\u2192D5/22)$ | DC | 0.086 | 0.086 | 0.131 |

DCG | 0.083 | 0.083 | ||

DCB | 0.083 | 0.083 | ||

Te $(P23\u2192P03)$ | DC | 0.589 | 0.589 | 0.584 |

DCG | 0.584 | 0.584 | ||

DCB | 0.584 | 0.584 |

The computed results shown in Table III show a good agreement with experimental measurements, with an error ranging from 10 meV (or 1.7%) in Te to 127 meV (or 13%) in Tl. Both 4C-CASSCF^{±} and 4C-DMRGSCF^{±} converge to the same results as expected, given the large value of *m*. Even though the valence space CAS calculations lack dynamical correlation, the error could be reduced further by increasing the size of the active space. For Ga atom, increasing the size of the active space to a CAS(3,26) improves the Ga splitting to 0.093 eV with DC-DMRGSCF^{±}(*m* = 150) compared to the experimental value of 0.102 eV. We defer the discussion on correlation convergence using the Dirac–Coulomb–Breit MCSCF approach to a future study as the goal of this section is to analyze the importance of the Breit Hamiltonian beyond the Coulomb operator for spectroscopic properties.

As we go down the Periodic Table, the inclusion of the Gaunt or the Breit operator consistently decreases the fine-structure splitting, e.g., by 2 meV for Ga to 10 meV for Tl. The improvement brought about by the Gaunt or Breit operator seems to deviate away from experiments except for Te. This is mainly because the correlation effect is not fully converged in a small active space. In other words, the “better” Dirac–Coulomb results are false positive due to unconverged electron correlation treatment.

Comparing the Gaunt and the Breit electron correlations, the difference for the atomic systems tested here is less than a meV. This agrees with our alkaline earth metal findings regarding the importance of correlated Breit treatments for valence electrons.

### D. Lanthanide contraction

The decreased atomic radii across the lanthanide series, termed “lanthanide contraction,” depends on treatment of relativistic effects.^{60–64} We calculate lanthanide contraction by difference in equilibrium bond length between La- and Lu-containing molecules.^{61} Previous estimates of 10%–30%^{63} and 9%–23%^{64} of the lanthanide contraction were due to relativistic effects. However, contributions beyond the four-component Dirac–Coulomb have not been shown. We investigate the effect of relativistic Hamiltonian on lanthanide contraction by studying monofluoride diatomics, LaF and LuF. For the *f*-block diatomics, an uncontracted cc-pVDZ-X2C basis set was used for lanthanides (Ln, where Ln is La or Lu) and an uncontracted cc-pVDZ basis set was used for fluorine.^{65,66} We used state-specific CAS(8,14) for LnF. The LnF active space consists of the orbitals derived from the Ln 6*s*, Ln 6*p*, and F 2*p*.

We report the equilibrium bond lengths for lanthanum fluoride diatomic molecules in Table IV. We find a lanthanide contraction of 0.109 Å for DC-CASSCF^{±}, a 0.003 Å difference from experiment. In going from DC-CASSCF^{±} to DCG-CASSCF^{±}, the lanthanide contraction error is improved by 0.001 Å. We estimate relativistic effects beyond the Dirac–Coulomb Hamiltonian to be 1% of the lanthanide contraction and less than 1% beyond the Gaunt interaction.

. | . | 4C-CASSCF^{±}
. | Expt.^{58,59}
. |
---|---|---|---|

LaF | DC | 2.059 | 2.023 |

DCG | 2.059 | ||

DCB | 2.059 | ||

LuF | DC | 1.950 | 1.917 |

DCG | 1.951 | ||

DCB | 1.951 | ||

Δ_{Re} | DC | 0.109 | 0.106 |

DCG | 0.108 | ||

DCB | 0.108 |

. | . | 4C-CASSCF^{±}
. | Expt.^{58,59}
. |
---|---|---|---|

LaF | DC | 2.059 | 2.023 |

DCG | 2.059 | ||

DCB | 2.059 | ||

LuF | DC | 1.950 | 1.917 |

DCG | 1.951 | ||

DCB | 1.951 | ||

Δ_{Re} | DC | 0.109 | 0.106 |

DCG | 0.108 | ||

DCB | 0.108 |

## IV. CONCLUSION

In this work, we introduce a fully correlated Dirac–Coulomb–Breit multiconfigurational self-consistent-field method within the frameworks of complete active space and density matrix renormalization group. In this approach, the Dirac–Coulomb–Breit Hamiltonian is included variationally in both the mean-field and correlated electron treatment.

Benchmark studies suggest that the positive-negative-energy orbital rotation in the no-virtual-pair approximation is important for deep-core electrons. Rotating only the positive-energy orbital space will lead to an overestimation of electron correlation by ∼10%, whereas no orbital rotation severely underestimates the correlation energy. Calculations on the correlation energy of valence electrons suggest that positive-negative orbital rotation becomes less important for electrons further away from the nucleus.

Calculations on He-like two-electron systems show that the Breit correlation becomes increasingly important as the atomic number increases and accounts for more than 50% of core electron correlation for heavier elements. On the other hand, we find no strong dependence of electron correlation on the Breit operator for valence electrons. Atomic fine structure study suggests that spectroscopic properties computed using only the correlated Dirac–Coulomb operator can produce false positive results. High-accuracy calculations require the fully correlated Gaunt or Breit operator and a large active space. Studies of lanthanide contraction due to the relativistic effect show excellent agreement with experiment.

This work lays the theoretical foundation for using the correlated Breit operator for high-accuracy computational chemistry within a multireference framework, representing the non-quantum-electrodynamics limit of the many-body theory. DMRG’s true potential in fast convergence toward the full CI limit is not fully explored in the proof-of-concept numerical examples presented herein and will be a subject of study for highly correlated chemical systems. The accuracy of the correlated Breit CASSCF and DMRG methods can be improved with additional dynamic correlations in the form of perturbation theory or multireference configuration interaction, which will be investigated in a future work.

## SUPPLEMENTARY MATERIAL

Convergence of correlation energies for Rn^{84+} and correlation energies for two-electron systems are given in the supplementary material.

## ACKNOWLEDGMENTS

This work was supported by 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 variational relativistic multi-reference methods and in the Computational and Theoretical program (Grant No. DE-SC0006863) for the development of the Breit operator for Dirac and quantum field theory.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

C.E.H., L.L., and H.H. contributed equally to this work.

**Chad E. Hoyer**: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Lixin Lu**: Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal). **Hang Hu**: Data curation (equal); Formal analysis (equal); Methodology (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (supporting). **Kirill D. Shumilov**: Data curation (equal); Formal analysis (equal). **Shichao Sun**: Methodology (equal). **Stefan Knecht**: Formal analysis (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Writing – original draft (supporting); Writing – review & editing (supporting). **Xiaosong Li**: Conceptualization (lead); Formal analysis (equal); Funding acquisition (lead); Investigation (equal); Methodology (lead); Project administration (lead); Resources (lead); Software (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead).

## DATA AVAILABILITY

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

## REFERENCES

_{3}

^{1}Π − X

^{1}Σ

^{+}transition

^{1}Σ states of ScH,YH, LaH, AcH, TmH, LuH and LrH