# Unified construction of relativistic Hamiltonians

Wenjian Liu<EMAIL_ADDRESS>Qingdao Institute for Theoretical and
Computational Sciences, Institute of Frontier and Interdisciplinary Science,
Shandong University, Qingdao, Shandong 266237, P. R. China

###### Abstract

It is shown that four-component (4C), quasi-four-component (Q4C), and exact
two-component (X2C) relativistic Hartree-Fock (HF) equations can be
implemented in an unified manner, by making use of the atomic nature of the
small components of molecular 4-spinors. A model density matrix approximation
can first be invoked for the small-component charge/current density functions,
which gives rise to a static, pre-molecular mean field (pmf) to be combined
with the one-electron term. As a result, only the nonrelativistic-like two-
electron term of the 4C/Q4C/X2C Fock matrix needs to be updated during the
iterations. A ‘one-center small-component’ approximation can then be invoked
in the evaluation of relativistic integrals. That is, all atom-centered small-
component basis functions are regarded as extremely localized nearby the
position of the atom to which they belong, such that they have vanishing
overlaps with all small- or large-component functions centered at other
nuclei. Under these approximations, the 4C, Q4C, and X2C mean-field and many-
electron Hamiltonians share precisely the same structure and accuracy. Beyond
these is the effective quantum electrodynamics Hamiltonian that can be
constructed in the same way. Such approximations lead to errors that are
orders of magnitude smaller than other sources of errors (e.g., truncation
errors in the one- and many-particle bases as well as uncertainties of
experimental measurements) and are hence safe to use for whatever purposes.
The quaternion forms of the 4C, Q4C, and X2C equations are also presented in
the most general way, based on which the corresponding Kramers-restricted
open-shell (RKOHF) variants are formulated for ‘high-spin’ open-shell systems.

##  I Introduction

The last two decades have witnessed fast progresses in relativistic quantum
chemistry, as evidenced by both theoretical developments[1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13, 14, 15] and efficient implementations in modern
computational software[16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29,
30, 31, 32, 33, 34, 35, 36, 37]. Much of the success is due to the advent of
the so-called exact through a one-step block-diagonalization of the matrix
Dirac equation, similar to the normalized elimination of the small component
(NESC)[38]. It is even simpler than approximate two-component (A2C)
approaches[39, 40, 41, 42] resulting from order-by-order expansions of the
operator Dirac equation and has hence become the main workhorse of
relativistic quantum chemistry. Noticeably, the various implementations[43,
44, 45, 46, 47, 48, 49] of the same X2C equation[50, 51] have adopted
seemingly different assumptions and even different terminologies. It is
therefore timely to establish an unified formulation to standardize the
implementations (see Sec. IV), such that different codes can produce precisely
the same results, just like the nonrelativistic ones. Not only so, the four-
component (4C)[52] and quasi-four-component (Q4C)[44] approaches also remain
to be unified into the same framework. To this end, the Dirac-Hartree-Fock
(DHF) equation represented in a restricted kinetically balanced (RKB)
basis[53] is first recapitulated in Sec. II, to establish the notations. For
completeness, the quaternion form of DHF is formulated in the most general
way, based on which a four-component Kramers-restricted open-shell HF (KROHF)
approach is proposed for the first time to describe ‘high-spin’ open-shell
systems with both double group and time-reversal symmetries (see Appendix A).
The projected four-component (P4C) equation, arising from a special
discretization[54, 55, 56, 57, 58, 59] of the operator DHF equation, is then
presented in Sec. III, which leads naturally to Q4C[60, 61] by invoking a
model density matrix (MDM) approximation for the small-component
charge/current density functions. This gives rise to a static, pre-molecular
mean field (pmf) that can be added to the one-electron term. As a result, only
the nonrelativistic-like two-electron term of Q4C needs to be updated during
the iterative cycles for self-consistency. That is, Q4C is four-component in
structure but two-component in computation, thereby justifying the name[44].
It turns out that the same MDM approximation can be applied to DHF and X2C as
well, so as to recast them to the same form as Q4C. The ‘extended atomic mean-
field’ (eamf) variant of X2C proposed recently by Knecht and coworkers[49] can
also be reproduced this way. Moreover, the effective, HF-like one-body
potential[62, 8] describing the leading quantum electrodynamics (QED) effect
can also be included from the outset, particularly when it is fitted into a
model spectral form[63, 64]. After having discussed in depth these mean-field
relativistic approaches, the corresponding many-body relativistic Hamiltonians
will be presented in Sec. V. The presentation is closed with concluding
remarks in Sec. VI.

Plain and boldface letters are used to denote operators and matrices,
respectively. The Einstein summation convention over repeated indices is
always employed.

##  II DHF

The Dirac operator-based relativistic mean-field description of a system of
$N$ electrons and $\tilde{N}$ positrons gives rise to the following energy
expression[65, 12]

$\displaystyle E_{\mathrm{ep}}$
$\displaystyle=\sum_{k}n_{k}h_{k}^{k}+\frac{1}{2}\sum_{k,l}n_{k}n_{l}\bar{g}^{kl}_{kl},\quad
k,l\in\mathrm{PES,NES},$ (1) $\displaystyle h_{p}^{q}$
$\displaystyle=h_{pq}=\langle\psi_{p}|h|\psi_{q}\rangle,$ (2) $\displaystyle
g_{pr}^{qs}$
$\displaystyle=g_{pq,rs}=(\psi_{p}\psi_{q}|V(1,2)|\psi_{r}\psi_{s}),\quad\bar{g}_{pr}^{qs}=g_{pr}^{qs}-g_{pr}^{sq},$
(3)

where $\\{n_{k}\\}$ are the occupation numbers of the 4-spinors
$\\{\psi_{k}\\}$: $+1$ for the $N$ occupied positive-energy states (PES), $-1$
for the $\tilde{N}$ occupied negative-energy states (NES), and 0 for the
remaining unoccupied PESs and NESs. The one-body operator $h$ is composed of
the Dirac operator $D$ and nuclear attraction $V$,

$\displaystyle h$ $\displaystyle=D+V,$ (4) $\displaystyle D$
$\displaystyle=c\bm{\alpha}\cdot\bm{p}+(\beta-1)c^{2},$ (5)
$\displaystyle\bm{\alpha}$ $\displaystyle=\begin{pmatrix}0&\bm{\sigma}\\\
\bm{\sigma}&0\end{pmatrix},\quad\beta=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},\quad\bm{p}=-\mathbbm{i}\hbar\bm{\partial},$ (6)
$\displaystyle\sigma_{x}$ $\displaystyle=\begin{pmatrix}0&1\\\
1&0\end{pmatrix},\quad\sigma_{y}=\begin{pmatrix}0&-\mathbbm{i}\\\
\mathbbm{i}&0\end{pmatrix},\quad\sigma_{z}=\begin{pmatrix}1&0\\\
0&-1\end{pmatrix},$ (7) $\displaystyle V(\bm{r})$
$\displaystyle=-\sum_{A}\frac{Z_{A}}{|\bm{r}-\bm{R}_{A}|}.$ (8)

The leading term in the electron-electron interaction $V(1,2)$ is the Coulomb
(C) interaction,

$V^{C}(1,2)=\frac{1}{r_{12}}=g_{12},$ (9)

which describes the electrostatic, charge-charge interaction between two
electrons. On top of this, the Gaunt (G) interaction,

$\displaystyle V^{G}(1,2)$
$\displaystyle=-\bm{\alpha}_{1}\cdot\bm{\alpha}_{2}g_{12},$ (10)

can further be included to account for the magnetic, current-current
interaction between two electrons. The combined CG interaction,
$V^{C}(1,2)+V^{G}(1,2)$, covers all two-electron spin-orbit, spin-spin, and
orbit-orbit couplings[66] and is therefore accurate enough for most chemical
problems. Nevertheless, the gauge term $V^{g}(1,2)$ in the Breit (B)
interaction $V^{B}(1,2)$,

$\displaystyle V^{B}(1,2)$ $\displaystyle=V^{G}(1,2)+V^{g}(1,2),$ (11)
$\displaystyle V^{g}(1,2)$
$\displaystyle=\frac{1}{2}g_{12}\bm{\alpha}_{1}\cdot\bm{\alpha}_{2}-\frac{1}{2}g_{12}^{3}(\bm{\alpha}_{1}\cdot\bm{r}_{12})(\bm{\alpha}_{2}\cdot\bm{r}_{12}),\quad\bm{r}_{12}=\bm{r}_{1}-\bm{r}_{2}$
(12)
$\displaystyle=-\frac{1}{2}V^{G}(1,2)-\frac{1}{2}\bm{\alpha}_{1}\cdot\mathbf{b}_{12}\cdot\bm{\alpha}_{2},$
(13) $\displaystyle\mathbf{b}_{12}$
$\displaystyle=\frac{\bm{r}_{12}\otimes\bm{r}_{12}}{r_{12}^{3}},\quad
b_{12}^{ij}=\frac{(\bm{r}_{12})_{i}(\bm{r}_{12})_{j}}{r_{12}^{3}},\quad i,j\in
x,y,z,$ (14)

is known to be important for deep core electrons[67, 68] and should hence be
further included in accurate descriptions of properties sampling the wave
function in the vicinity of a nucleus. Overall, the electron-electron
interaction can be written as

$\displaystyle V(1,2)$
$\displaystyle=g_{12}+\mathrm{c_{g}}g_{12}\bm{\alpha}_{1}\cdot\bm{\alpha}_{2}+\mathrm{c_{b}}\bm{\alpha}_{1}\cdot\mathbf{b}_{12}\cdot\bm{\alpha}_{2},$
(15)

which reduces to the C, CG, and CB interactions by setting
$(\mathrm{c_{g}},\mathrm{c_{b}})$ to $(0,0)$, $(-1,0)$, and $(-1/2,-1/2)$,
respectively. To avoid over notation, we will refer to the second term in Eq.
(13) simply as the gauge term.

To minimize the energy $E_{\mathrm{ep}}$ (1) subject to the orthonormal
conditions $\langle\psi_{p}|\psi_{q}\rangle=\delta_{pq}$, we introduce the
following canonical Lagrangian

$\displaystyle
L=E_{\mathrm{ep}}-\sum_{k}n_{k}[\langle\psi_{k}|\psi_{k}\rangle-1]\epsilon_{k},\quad
k\in\mathrm{PES,NES}.$ (16)

The condition $\frac{\delta L}{\delta\psi_{i}^{\dagger}}=0$ then gives rise to

$\displaystyle Fn_{i}|\psi_{i}\rangle$
$\displaystyle=\epsilon_{i}n_{i}|\psi_{i}\rangle,\quad i\in\mathrm{PES,NES},$
(17)

where the Fock operator $F$ reads[60]

$\displaystyle F$ $\displaystyle=h+\sum_{k}n_{k}\bar{g}_{\cdot k}^{\cdot
k},\quad k\in\mathrm{PES,NES}$ (18)
$\displaystyle=\begin{pmatrix}h^{LL}&h^{LS}\\\
h^{SL}&h^{SS}\end{pmatrix}+\begin{pmatrix}G^{LL}&G^{LS}\\\
G^{SL}&G^{SS}\end{pmatrix}=\begin{pmatrix}F^{LL}&F^{LS}\\\
F^{SL}&F^{SS}\end{pmatrix},$ (19) $\displaystyle h^{LL}$
$\displaystyle=V,\quad h^{LS}=c\bm{\sigma}\cdot\bm{p},\quad
h^{SL}=c\bm{\sigma}\cdot\bm{p},\quad h^{SS}=V-2c^{2},$ (20) $\displaystyle
G^{LL}$
$\displaystyle=J^{C,LL}+J^{C,SS}-K^{C,LL}-\mathrm{c_{g}}K^{G,SS}-\mathrm{c_{b}}K^{g,SS}=G^{LL{\dagger}},$
(21) $\displaystyle G^{LS}$
$\displaystyle=-K^{C,LS}+\mathrm{c_{g}}(J^{G,LS}+J^{G,SL}-K^{G,SL})+\mathrm{c_{b}}(J^{g,LS}+J^{g,SL}-K^{g,SL})=G^{SL{\dagger}},$
(22) $\displaystyle G^{SL}$
$\displaystyle=-K^{C,SL}+\mathrm{c_{g}}(J^{G,LS}+J^{G,SL}-K^{G,LS})+\mathrm{c_{b}}(J^{g,LS}+J^{g,SL}-K^{g,LS})=G^{LS{\dagger}},$
(23) $\displaystyle G^{SS}$
$\displaystyle=J^{C,LL}+J^{C,SS}-K^{C,SS}-\mathrm{c_{g}}K^{G,LL}-\mathrm{c_{b}}K^{g,LL}=G^{SS{\dagger}}.$
(24)

As it stands, Eq. (17) determines only the occupied PESs and NESs but which
can be extended to the unoccupied ones, viz.,

$\displaystyle F|\psi_{p}\rangle$
$\displaystyle=\epsilon_{p}|\psi_{p}\rangle,\quad p\in\mathrm{PES,NES}.$ (25)

The validity of the canonical Lagrangian (16) stems from the fact that the
Fock operator $F$ (18) is Hermitian, such that
$\langle\psi_{p}|\psi_{q}\rangle=0$ holds automatically for $p\neq q$. The
energetically lowest $N$ PESs and highest $\tilde{N}$ NESs are to be occupied
in each iteration when solving Eq. (25) iteratively. That is, the non-Aufbau,
interior roots of Eq. (25) are to be sought in a self-consistent manner. Note
in passing that Eq. (1) and hence Eq. (25) are approximations to the mean-
field QED theory[8, 12] of electrons and positrons by neglecting the vacuum
polarization and electron self-energy. Further setting the occupation numbers
to zero for all NESs, Eq. (1)/(25) reduces to the usual Dirac-Hartee-Fock
(DHF) energy/stationarity condition for $N$ electrons alone. In this case, the
expressions (1) and (25) can also be derived[52] by starting with the empty
Dirac picture (where all NESs are assumed to be unoccupied), instead of the
filled Dirac picture underlying the QED formulation[8, 12]. Yet, it should be
kept in mind that the agreement between the empty and filled Dirac pictures
holds only for one-electron properties, but not for any two-body property
(including electron correlation)[69, 62].

In practice, the molecular 4-spinors (M4S) $\\{\psi_{p}\\}$ have to be
expanded in a suitably chosen basis. Among the various prescriptions[70], the
so-called restricted kinetic balance (RKB)[53] is most recommended due to its
simplicity, viz.,

$\displaystyle|\psi_{p}\rangle$
$\displaystyle=\sum_{\mu=1}^{n}|\xi_{\mu}\rangle C_{\mu p},\quad p\in[1,4n],$
(26) $\displaystyle|\xi_{\mu}\rangle$
$\displaystyle=\begin{pmatrix}|\chi^{L}_{\mu}\rangle&0\\\
0&|\chi_{\mu}^{S}\rangle\end{pmatrix},$ (27)
$\displaystyle|\chi^{L}_{\mu}\rangle$
$\displaystyle=\begin{pmatrix}|g_{\mu}\rangle&0\\\
0&|g_{\mu}\rangle\end{pmatrix}=\sigma_{0}|g_{\mu}\rangle,\quad\sigma_{0}=\mathbf{I}_{2},$
(28) $\displaystyle|\chi^{S}_{\mu}\rangle$
$\displaystyle=\frac{1}{2c}\Pi|\chi^{L}_{\mu}\rangle=\frac{1}{2c}\Pi|g_{\mu}\rangle=-|\chi^{S}_{\mu}\rangle^{\dagger},\quad\Pi=\bm{\sigma}\cdot\bm{p},$
(29) $\displaystyle C_{\mu p}$ $\displaystyle=\begin{pmatrix}C_{\mu p}^{L}\\\
C_{\mu p}^{S}\end{pmatrix},\quad C^{X}_{\mu p}=\begin{pmatrix}C^{X\alpha}_{\mu
p}\\\ C^{X\beta}_{\mu p}\end{pmatrix},\quad X\in L,S,$ (30)

where $\\{g_{\mu}\\}_{\mu=1}^{n}$ are prechosen $n$ scalar functions. Note
that each of the four columns of $|\xi_{\mu}\rangle$ (27) corresponds to an
independent four-component (4C) basis function, such that there are in total
$4n$ basis functions and hence $4n$ M4Ss $\\{|\psi_{p}\rangle\\}$, half of
which are PESs and the other half of which are NESs. In terms of such a basis
expansion, Eq. (25) is converted to an algebraic equation,

$\displaystyle\mathbf{F}\mathbf{C}$
$\displaystyle=\mathbf{M}\mathbf{C}\mathbf{E},$ (31) $\displaystyle
F_{\mu\nu}$ $\displaystyle=h_{\mu\nu}+G_{\mu\nu}[\mathbf{D}],\quad
M_{\mu\nu}=\langle\xi_{\mu}|\xi_{\nu}\rangle,$ (32) $\displaystyle h_{\mu\nu}$
$\displaystyle=\langle\xi_{\mu}|h|\xi_{\nu}\rangle,$ (33) $\displaystyle
G_{\mu\nu}[\mathbf{D}]$
$\displaystyle=(J^{C}_{\mu\nu}[\mathbf{D}]-K^{C}_{\mu\nu}[\mathbf{D}])+\mathrm{c_{g}}(J^{G}_{\mu\nu}[\mathbf{D}]-K^{G}_{\mu\nu}[\mathbf{D}])+\mathrm{c_{b}}(J^{g}_{\mu\nu}[\mathbf{D}]-K^{g}_{\mu\nu}[\mathbf{D}])$
(34)
$\displaystyle=[(\Omega_{\mu\nu}|g_{12}|\rho^{4c})-(\Omega_{\mu\lambda}|g_{12}|\mathbf{D}_{\lambda\kappa}\Omega_{\kappa\nu})]$
$\displaystyle+\mathrm{c_{g}}[(\bm{\Xi}_{\mu\nu}|g_{12}\cdot|\bm{j}^{4c})-(\bm{\Xi}_{\mu\lambda}|g_{12}\cdot|\mathbf{D}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu})]$
$\displaystyle+\mathrm{c_{b}}[(\bm{\Xi}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\bm{j}^{4c})-(\bm{\Xi}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\mathbf{D}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu})],$
(35) $\displaystyle\mathbf{D}$
$\displaystyle=\mathbf{C}_{+}\mathbf{n}\mathbf{C}_{+}^{\dagger}=\begin{pmatrix}\mathbf{C}_{+}^{L}\\\
\mathbf{C}_{+}^{S}\end{pmatrix}\mathbf{n}\begin{pmatrix}\mathbf{C}_{+}^{L{\dagger}}&\mathbf{C}_{+}^{S{\dagger}}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{D}^{LL}&\mathbf{D}^{LS}\\\
\mathbf{D}^{SL}&\mathbf{D}^{SS}\end{pmatrix},\quad\mathbf{D}^{XY}=\mathbf{C}^{X}_{+}\mathbf{n}C^{Y{\dagger}}_{+},\quad
X,Y=L,S,$ (36)

where $\mathbf{n}$ is the diagonal occupation number matrix, whereas the
overlap charge and current distribution functions are defined
as[ReSpectJCP2020]

$\displaystyle\Omega_{\mu\nu}$
$\displaystyle=\xi_{\mu}^{\dagger}\xi_{\nu}=\begin{pmatrix}\Omega_{\mu\nu}^{LL}&0\\\
0&\Omega_{\mu\nu}^{SS}\end{pmatrix}=\begin{pmatrix}\chi_{\mu}^{L{\dagger}}\chi_{\nu}^{L}&0\\\
0&\chi_{\mu}^{S{\dagger}}\chi_{\nu}^{S}\end{pmatrix}$ (37)
$\displaystyle=\frac{1}{4c^{2}}\begin{pmatrix}4c^{2}\sigma_{0}g_{\mu}g_{\nu}&0\\\
0&\sigma_{0}g_{\mu}^{i}g_{\nu}^{i}+\mathbbm{i}\epsilon_{ijk}\sigma_{i}g_{\mu}^{j}g_{\nu}^{k}\end{pmatrix},\quad
i,j,k\in x,y,z,$ (38) $\displaystyle\Xi_{\mu\nu}^{i}$
$\displaystyle=\xi_{\mu}^{\dagger}\alpha_{i}\xi_{\nu}=\begin{pmatrix}0&\Xi_{\mu\nu}^{i,LS}\\\
\Xi_{\mu\nu}^{i,SL}&0\end{pmatrix}=\begin{pmatrix}0&\chi_{\mu}^{L{\dagger}}\sigma_{i}\chi_{\nu}^{S}\\\
\chi_{\mu}^{S{\dagger}}\sigma_{i}\chi_{\nu}^{L}&0\end{pmatrix}$ (39)
$\displaystyle=\frac{\mathbbm{i}}{2c}\begin{pmatrix}0&-\sigma_{0}g_{\mu}g_{\nu}^{i}-\mathbbm{i}\epsilon_{ijk}g_{\mu}g_{\nu}^{j}\sigma_{k}\\\
\sigma_{0}g_{\mu}^{i}g_{\nu}-\mathbbm{i}\epsilon_{ijk}g_{\mu}^{j}g_{\nu}\sigma_{k}&0\end{pmatrix},\quad
i,j,k\in x,y,z,$ (40)

where $g_{\mu}^{i}=\partial_{i}g_{\mu}$ ($i\in x,y,z$). The number charge and
current densities can then be calculated as

$\displaystyle\rho^{4C}$
$\displaystyle=n_{i}\psi_{i}^{{\dagger}}\psi_{i}=n_{i}\psi_{i}^{L{\dagger}}\psi_{i}^{L}+n_{i}\psi_{i}^{S{\dagger}}\psi_{i}^{S}=\mathrm{Tr}[\Omega^{LL}\mathbf{D}^{LL}]+\mathrm{Tr}[\Omega^{SS}\mathbf{D}^{SS}]=\mathrm{Tr}[\Omega\mathbf{D}],$
(41) $\displaystyle\bm{j}^{4C}$
$\displaystyle=n_{i}\psi_{i}^{{\dagger}}\bm{\alpha}\psi_{i}=n_{i}\psi_{i}^{L{\dagger}}\bm{\sigma}\psi_{i}^{S}+n_{i}\psi_{i}^{S{\dagger}}\bm{\sigma}\psi_{i}^{L}=\mathrm{Tr}[\bm{\Xi}^{LS}\mathbf{D}^{SL}]+\mathrm{Tr}[\bm{\Xi}^{SL}\mathbf{D}^{LS}]=\mathrm{Tr}[\bm{\Xi}\mathbf{D}].$
(42)

As already indicated by Eqs. (37) and (39), Eq. (31) can be recast into block
form

$\displaystyle\begin{pmatrix}\mathbf{F}^{LL}&\mathbf{F}^{LS}\\\
\mathbf{F}^{SL}&\mathbf{F}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{C}^{L}_{+}&\mathbf{C}^{L}_{-}\\\
\mathbf{C}^{S}_{+}&\mathbf{C}^{S}_{-}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{M}^{LL}&\mathbf{0}\\\
\mathbf{0}&\mathbf{M}^{SS}\end{pmatrix}\begin{pmatrix}\mathbf{C}^{L}_{+}&\mathbf{C}^{L}_{-}\\\
\mathbf{C}^{S}_{+}&\mathbf{C}^{S}_{-}\end{pmatrix}\begin{pmatrix}\mathbf{E}_{+}&\mathbf{0}\\\
\mathbf{0}&\mathbf{E}_{-}\end{pmatrix},$ (43)
$\displaystyle\begin{pmatrix}\mathbf{F}^{LL}&\mathbf{F}^{LS}\\\
\mathbf{F}^{SL}&\mathbf{F}^{SS}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mathbf{h}^{LL}&\mathbf{h}^{LS}\\\
\mathbf{h}^{SL}&\mathbf{h}^{SS}\end{pmatrix}+\begin{pmatrix}\mathbf{G}^{LL}[\mathbf{D}]&\mathbf{G}^{LS}[\mathbf{D}]\\\
\mathbf{G}^{SL}[\mathbf{D}]&\mathbf{G}^{SS}[\mathbf{D}]\end{pmatrix},$ (44)

where

$\displaystyle M^{LL}_{\mu\nu}$
$\displaystyle=\langle\chi^{L}_{\mu}|\chi^{L}_{\nu}\rangle=\sigma_{0}S_{\mu\nu},\quad
S_{\mu\nu}=\langle g_{\mu}|g_{\nu}\rangle,$ (45) $\displaystyle
M^{SS}_{\mu\nu}$
$\displaystyle=\langle\chi^{S}_{\mu}|\chi^{S}_{\nu}\rangle=\frac{1}{2c^{2}}\sigma_{0}T_{\mu\nu},\quad
T_{\mu\nu}=\frac{1}{2}\langle g_{\mu}|\bm{p}^{2}|g_{\nu}\rangle,$ (46)
$\displaystyle h^{LL}_{\mu\nu}$
$\displaystyle=\langle\chi_{u}^{L}|V|\chi_{\nu}^{L}\rangle=\sigma_{0}V_{\mu\nu},\quad
V_{\mu\nu}=\langle g_{\mu}|V|g_{\nu}\rangle,$ (47) $\displaystyle
h^{LS}_{\mu\nu}$
$\displaystyle=\langle\chi_{u}^{L}|c\Pi|\chi_{\nu}^{S}\rangle=\sigma_{0}T_{\mu\nu}=h^{SL*}_{\nu\mu},$
(48) $\displaystyle h^{SS}_{\mu\nu}$
$\displaystyle=\langle\chi_{u}^{S}|V-2c^{2}|\chi_{\nu}^{S}\rangle=\frac{1}{4c^{2}}W_{\mu\nu}-\sigma_{0}T_{\mu\nu},\quad
W_{\mu\nu}=\langle\Pi g_{\mu}|V|\Pi g_{\nu}\rangle,$ (49) $\displaystyle
G^{LL}_{\mu\nu}[\mathbf{D}]$
$\displaystyle=(\mu^{L}\nu^{L}|g_{12}|i^{L}i^{L})-(\mu^{L}i^{L}|g_{12}|i^{L}\nu^{L})+(\mu^{L}\nu^{L}|g_{12}|i^{S}i^{S})$
$\displaystyle-\mathrm{c_{g}}(\mu^{L}\bm{\sigma}i^{S}|g_{12}\cdot|i^{S}\bm{\sigma}\nu^{L})-\mathrm{c_{b}}(\mu^{L}\bm{\sigma}i^{S}|\cdot\mathbf{b}_{12}\cdot|i^{S}\bm{\sigma}\nu^{L})$
(50)
$\displaystyle=(\Omega_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\mathbf{D}^{LL}_{\lambda\kappa}\Omega^{LL}_{\kappa\lambda}])-(\Omega_{\mu\lambda}^{LL}|g_{12}|\mathbf{D}^{LL}_{\lambda\kappa}\Omega^{LL}_{\kappa\nu})+(\Omega_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\mathbf{D}^{SS}_{\lambda\kappa}\Omega^{SS}_{\kappa\lambda}])$
$\displaystyle-\mathrm{c_{g}}(\bm{\Xi}_{\mu\lambda}^{LS}|g_{12}\cdot|\mathbf{D}^{SS}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu}^{SL})-\mathrm{c_{b}}(\bm{\Xi}_{\mu\lambda}^{LS}|\cdot\mathbf{b}_{12}\cdot|\mathbf{D}^{SS}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu}^{SL}),$
(51) $\displaystyle G^{LS}_{\mu\nu}[\mathbf{D}]$
$\displaystyle=-(\mu^{L}i^{L}|g_{12}|i^{S}\nu^{S})$
$\displaystyle+\mathrm{c_{g}}[(\mu^{L}\bm{\sigma}\nu^{S}|g_{12}\cdot|i^{L}\bm{\sigma}i^{S})-(\mu^{L}\bm{\sigma}i^{S}|g_{12}\cdot|i^{L}\bm{\sigma}\nu^{S})+(\mu^{L}\bm{\sigma}\nu^{S}|g_{12}\cdot|i^{S}\bm{\sigma}i^{L})]$
$\displaystyle+\mathrm{c_{b}}[(\mu^{L}\bm{\sigma}\nu^{S}|\cdot\mathbf{b}_{12}\cdot|i^{L}\bm{\sigma}i^{S})-(\mu^{L}\bm{\sigma}i^{S}|\cdot\mathbf{b}_{12}\cdot|i^{L}\bm{\sigma}\nu^{S})+(\mu^{L}\bm{\sigma}\nu^{S}|\cdot\mathbf{b}_{12}\cdot|i^{S}\bm{\sigma}i^{L})]$
(52)
$\displaystyle=-(\Omega^{LL}_{\mu\lambda}|g_{12}|\mathbf{D}^{LS}_{\lambda\kappa}\Omega^{SS}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[(\bm{\Xi}_{\mu\nu}^{LS}|g_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\lambda}])-(\bm{\Xi}^{LS}_{\mu\lambda}|g_{12}\cdot|\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\nu})+(\bm{\Xi}^{LS}_{\mu\nu}|g_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\lambda}])]$
$\displaystyle+\mathrm{c_{b}}[\bm{\Xi}_{\mu\nu}^{LS}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\lambda}])-(\bm{\Xi}^{LS}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\nu})+(\bm{\Xi}^{LS}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\lambda}])],$
(53) $\displaystyle G^{SL}_{\mu\nu}[\mathbf{D}]$
$\displaystyle=-(\mu^{S}i^{S}|g_{12}|i^{L}\nu^{L})$
$\displaystyle+\mathrm{c_{g}}[(\mu^{S}\bm{\sigma}\nu^{L}|g_{12}\cdot|i^{S}\bm{\sigma}i^{L})-(\mu^{S}\bm{\sigma}i^{L}|g_{12}\cdot|i^{S}\bm{\sigma}\nu^{L})+(\mu^{S}\bm{\sigma}\nu^{L}|g_{12}\cdot|i^{L}\bm{\sigma}i^{S})]$
$\displaystyle+\mathrm{c_{b}}[(\mu^{S}\bm{\sigma}\nu^{L}|\cdot\mathbf{b}_{12}\cdot|i^{S}\bm{\sigma}i^{L})-(\mu^{S}\bm{\sigma}i^{L}|\cdot\mathbf{b}_{12}\cdot|i^{S}\bm{\sigma}\nu^{L})+(\mu^{S}\bm{\sigma}\nu^{L}|\cdot\mathbf{b}_{12}\cdot|i^{L}\bm{\sigma}i^{S})]$
(54)
$\displaystyle=-(\Omega^{SS}_{\mu\lambda}|g_{12}|\mathbf{D}^{SL}_{\lambda\kappa}\Omega^{LL}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[\bm{\Xi}_{\mu\nu}^{SL}|g_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\lambda}])-(\bm{\Xi}^{SL}_{\mu\lambda}|g_{12}\cdot|\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\nu})+(\bm{\Xi}^{SL}_{\mu\nu}|g_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\lambda}])]$
$\displaystyle+\mathrm{c_{b}}[(\bm{\Xi}_{\mu\nu}^{SL}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\lambda}])-(\bm{\Xi}^{SL}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\mathbf{D}^{LS}_{\lambda\kappa}\bm{\Xi}^{SL}_{\kappa\nu})+(\bm{\Xi}^{SL}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\mathbf{D}^{SL}_{\lambda\kappa}\bm{\Xi}^{LS}_{\kappa\lambda}])],$
(55) $\displaystyle G^{SS}_{\mu\nu}[\mathbf{D}]$
$\displaystyle=(\mu^{S}\nu^{S}|g_{12}|i^{S}i^{S})-(\mu^{S}i^{S}|g_{12}|i^{S}\nu^{S})+(\mu^{S}\nu^{S}|g_{12}|i^{L}i^{L})$
$\displaystyle-\mathrm{c_{g}}(\mu^{S}\bm{\sigma}i^{L}|g_{12}\cdot|i^{L}\bm{\sigma}\nu^{S})-\mathrm{c_{b}}(\mu^{S}\bm{\sigma}i^{L}|\cdot\mathbf{b}_{12}\cdot|i^{L}\bm{\sigma}\nu^{S})$
(56)
$\displaystyle=(\Omega_{\mu\nu}^{SS}|g_{12}|\mathrm{Tr}_{2}[\mathbf{D}^{SS}_{\lambda\kappa}\Omega^{SS}_{\kappa\lambda}])-(\Omega_{\mu\lambda}^{SS}|g_{12}|\mathbf{D}^{SS}_{\lambda\kappa}\Omega^{SS}_{\kappa\nu})+(\Omega_{\mu\nu}^{SS}|g_{12}|\mathrm{Tr}_{2}[\mathbf{D}^{LL}_{\lambda\kappa}\Omega^{LL}_{\kappa\lambda}])$
$\displaystyle-\mathrm{c_{g}}(\bm{\Xi}_{\mu\lambda}^{SL}|g_{12}\cdot|\mathbf{D}^{LL}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu}^{LS})-\mathrm{c_{b}}(\bm{\Xi}_{\mu\lambda}^{SL}|\cdot\mathbf{b}_{12}\cdot|\mathbf{D}^{LL}_{\lambda\kappa}\bm{\Xi}_{\kappa\nu}^{LS}).$
(57)

In the above notation, $\mathbf{D}^{XY}$, $\Omega^{XY}$, and $\bm{\Xi}^{XY}$
($X,Y=L,S$) are all 2-by-2 block matrices, each block of which is a $n$-by-$n$
matrix. Alternatively, every element $\mathbf{D}^{XY}_{\mu\nu}$,
$\Omega^{XY}_{\mu\nu}$ or $\bm{\Xi}^{XY}_{\mu\nu}$ can be understood as a
2-by-2 matrix over spin labels. Therefore, the symbol $\mathrm{Tr}_{2}$ means
trace over the spin degrees of freedom, e.g.,

$\displaystyle\mathrm{Tr}_{2}[\mathbf{D}^{XY}_{\lambda\kappa}\Omega^{YX}_{\kappa\lambda}]$
$\displaystyle=\mathrm{Tr}_{2}[\Omega^{YX}_{\kappa\lambda}\mathbf{D}^{XY}_{\lambda\kappa}],\quad\lambda,\kappa\in[1,n]$
(58) $\displaystyle=D^{X\alpha Y\alpha}_{\lambda\kappa}\Omega^{Y\alpha
X\alpha}_{\kappa\lambda}+D^{X\alpha Y\beta}_{\lambda\kappa}\Omega^{Y\beta
X\alpha}_{\kappa\lambda}+D^{X\beta Y\alpha}_{\lambda\kappa}\Omega^{Y\alpha
X\beta}_{\kappa\lambda}+D^{X\beta Y\beta}_{\lambda\kappa}\Omega^{Y\beta
X\beta}_{\kappa\lambda}.$ (59)

Note in particular that $\mathbf{D}^{XY}_{\lambda\kappa}$ and
$\Omega^{YZ}_{\kappa\lambda}/\bm{\Xi}^{XY}_{\kappa\lambda}$ do not commute in
block form (i.e., $\lambda,\kappa\in[1,n]$). Nevertheless, they do commute
with each other when fully expanded (i.e., $\lambda,\kappa\in[1,2n]$). Given
the density matrix $\mathbf{D}$ (36), the energy $E_{\mathrm{ep}}$ (1) can be
calculated as

$\displaystyle E_{\mathrm{ep}}$
$\displaystyle=\frac{1}{2}\mathrm{Tr}(\mathbf{h}+\mathbf{F})\mathbf{D}.$ (60)

Some remarks are in order:

  1. (I)

The RKB prescription (29) for generating $2n$ small-component basis functions
$\\{\chi_{\mu}^{S}\\}$ directly from the $n$-scalar functions
$\\{g_{\mu}\\}_{\mu=1}^{n}$ guarantees the correct nonrelativistic limit (nrl)
of the PESs but the NESs are still in error of $\mathcal{O}(c^{0})$[70]. As
such, it does not provide full variational safety. Depending very much on the
construction of $\\{g_{\mu}\\}_{\mu=1}^{n}$, some bounds failures (or
prolapse[71]) of $O(c^{-4})$ may occur. Nevertheless, such bounds failures
will diminish when approaching to the basis set limit, at a rate that is not
much different from the nonrelativistic counterpart[72]. Another point that
deserves to be mentioned is that, when $\\{g_{\mu}\\}_{\mu=1}^{n}$ are
spherical Gaussians, the principal quantum number must be set to the angular
momentum $l$ plus one (i.e., $1s$, $2p$, $3d$, $4f$, $5g$, etc.). Otherwise,
terrible variational collapse would occur[70]. Although this is usually the
default option, it is not mandatory in the nonrelativistic case.

  2. (II)

It is the appearance of small-component basis functions that renders the DHF
calculation very expensive. As shown in Appendix A, there are in total 25
(325) real-valued scalar integrals to evaluate and process just for a single
term of the two-electron matrix element $G_{\mu\nu}[\mathbf{D}]$ under the
Coulomb (Gaunt/Breit) interaction. However, the situation is not really that
bad. Since the small component $\psi_{p}^{S}(\bm{r})$
($=\sum_{A}\sum_{\mu_{A}}\chi_{\mu_{A}}^{S}(\bm{r}-\bm{R}_{A})C_{\mu_{A}p}^{S}$)
of a M4S $\psi_{p}(\bm{r})$ is appreciable only nearby the positions of the
nuclei[73], the products
$C_{\mu_{A}p}^{S{\dagger}}\chi_{\mu_{A}}^{S{\dagger}}(\bm{r}-\bm{R}_{A})\chi_{\nu_{B}}^{S}(\bm{r}-\bm{R}_{B})C_{\nu_{B}q}^{S}$
in $\psi_{p}^{S{\dagger}}(\bm{r})\psi_{q}^{S}(\bm{r})$ are negligibly small
when $A$ and $B$ refer to different atoms in the molecule (NB:
$\chi_{\mu_{A}}^{S}$ denotes a small-component function $\chi_{\mu}^{S}$
centered at the position $\bm{R}_{A}$ of atom $A$). An immediate deduction is
that the small-component distribution functions $\Omega^{SS}_{\mu\nu}$ (of
$\mathcal{O}(c^{-2})$) and density matrix elements $\mathbf{D}^{SS}_{\nu\mu}$
can be confined to one-centered only, i.e.,

$\displaystyle\Omega^{SS}_{\mu_{A}\nu_{B}}(\bm{r})\approx\Omega^{SS}_{\mu_{A}\nu_{B}}(\bm{r})\delta_{AB},\quad\mathbf{D}^{SS}\approx\sum^{\oplus}_{A}\mathbf{D}^{SS}_{A},\quad\mathbf{D}^{SS}_{A}=\mathbf{C}_{A,+}^{S}\mathbf{n}_{A}\mathbf{C}_{A,+}^{S{\dagger}},$
(61)

thereby leading to

$\displaystyle\rho^{S}=\mathrm{Tr}[\Omega^{SS}\mathbf{D}^{SS}]\approx\tilde{\rho}^{S}=\sum_{A}\rho^{S}_{A}$
(62)

for the molecular small-component density. Such approximations imply the
following relations[74]

$\displaystyle(L_{A}L_{B}|S_{C}S_{D})$
$\displaystyle\approx(L_{A}L_{B}|S_{C}S_{D})\delta_{CD},$
$\displaystyle(S_{A}S_{B}|L_{C}L_{D})$
$\displaystyle\approx(S_{A}S_{B}|L_{C}L_{D})\delta_{AB},$
$\displaystyle(S_{A}S_{B}|S_{C}S_{D})$
$\displaystyle\approx(S_{A}S_{B}|S_{C}S_{D})\delta_{AB}\delta_{CD}.$ (63)

for the electron repulsion integrals (ERI). The $(S_{A}S_{A}|S_{B}S_{B})$ type
of ERIs, which are of $\mathcal{O}(c^{-4})$, can further be confined to one-
centered only, i.e.,

$\displaystyle(S_{A}S_{A}|S_{B}S_{B})\approx(S_{A}S_{A}|S_{B}S_{B})\delta_{AB}.$
(64)

The situation is even better when atomic 4-spinors (A4S; from DHF calculations
of averaged and nonpolarized atomic configurations) are taken as the basis[54,
55, 56, 57, 58, 59], for their small components do not overlap each other
discernibly when they come from different atoms. Similarly, the Gaunt/Breit
integrals $(LS|SL)$, $(SL|LS)$, $(LS|LS)$, and $(SL|SL)$ can be confined to at
most two-centered,

$\displaystyle(X_{A}Y_{B}|Y_{C}X_{D})$
$\displaystyle\approx(X_{A}Y_{B}|Y_{C}X_{D})\delta_{AB}\delta_{CD},\quad
X,Y=L,S,\quad X\neq Y,$ $\displaystyle(X_{A}Y_{B}|X_{C}Y_{D})$
$\displaystyle\approx(X_{A}Y_{B}|X_{C}Y_{D})\delta_{AB}\delta_{CD},\quad
X,Y=L,S,\quad X\neq Y.$ (65)

Eqs. (63) to (65) can be termed collectively ‘one-center small-component’
(1CSC) approximation (i.e., a small-component function is always regarded as
extremely localized nearby the position of a nucleus, such that it has
vanishing overlaps with all functions centered at other nuclei). Since such
ERIs involve at most three centers, the use of local symmetries (planar,
cylindrical, and spherical symmetries for three-, two-, and one-center
integrals) would be very beneficial. Moreover, since only the spin-free part
of the whole gauge term $V^{g}$ in Eq. (13) contributes to the energy at
$\mathcal{O}(c^{-2})$ (which is even simpler than that of $V^{G}$)[66], It can
be anticipated that only the spin-free part of $V^{g}$ needs to be included in
practical applications[75]. The surving ERIs can further be evaluated
efficiently with the RKB-based resolution of the identity (RKB-RI)[76].

  3. (III)

Only the $N_{o}$ ($=N+\tilde{N}\ll 4n$) occupied states need to be determined
when solving Eq. (31)/(43) iteratively. Such states appear as a narrow,
interior portion of the whole orbital spectrum. Therefore, it is of great
value to invoke an algorithm that can directly access such interior roots. The
recently proposed iterative vector interaction (iVI)[77, 78] is one of such
algorithms. It is not only very robust but also memory efficient, by working
with a search space of fixed dimension that is determined automatically by the
number of target states or a preset energy window. A speedup factor of
$(4n)^{3}/[N_{o}(4n)^{2}]=4n/N_{o}\gg 1$ can be gained as compared with the
full matrix diagonalization.

The above presentation provides the basics for an efficient implementation of
DHF. Additional gain in efficiency stems from the full use of double point
group and time reversal symmetries, which can be achieved in two distinct
ways, by constructing Kramers-paired double group symmetry functions according
to fermion irreducible representations (irrep)[79, 80, 81, 82, 83] or by
combining quaternion algebra (which incorporates time reversal symmetry[84,
85]) with corresponding boson irreps for the real and imaginary parts of each
component of a quaternion spinor[86, 87]. The former is suited for a
$j$-adapted spinor basis, whereas the latter is suited for a scalar basis. The
quaternion form of the DHF equation (31) underlying the latter was
formulated[86] only for closed-shell systems under the Coulomb interaction. As
a matter of fact, it can be formulated more generally (see Appendix A) via a
quaternion unitary transformation[88] that can block-diagonalize any matrix
(e.g., open-shell Fock matrix with the Gaunt/Breit interaction) with identical
diagonal blocks. Based on such a general quaternion DHF, a four-component
relativistic Kramers restricted open-shell DHF (KROHF) scheme for ‘high-spin’
open-shell systems can be formulated (see Appendix A), in the same way as the
two-component counterpart[89].

##  III Q4C

The DHF equation (31)/(43) generates both PESs and NESs. The latter are
usually discarded at a correlated level under the no-pair approximation (NPA).
The question is how to avoid the molecular NESs from the outset if the NPA is
doomed to make. Actually, this can be done in two different ways[45]: freeze
or remove the NESs. The former stays in the four-component framework but can
be made operationally the same as a two-component theory, whereas the latter
works with a two-component, electron-only Hamiltonian. In essence, the two
paradigms stem from the same physical origin and can hence be made and have
been made completely equivalent in terms of simplicity, accuracy, and
efficiency[45].

To realize the first paradigm, we first examine the $S/L$ ratio between the
small and large components of a PES $\psi_{i}$, which can be obtained from the
second row of Eq. (17),

$\displaystyle\psi_{i}^{S}$
$\displaystyle=(\epsilon_{i}-F^{SS})^{-1}F^{SL}\psi_{i}^{L}$ (66)
$\displaystyle\triangleq(\epsilon_{i}-V^{SS}_{\mathrm{eff}}+2c^{2})^{-1}(c\bm{\sigma}\cdot\bm{p}+V^{SL}_{\mathrm{eff}})\psi_{i}^{L}$
(67)
$\displaystyle\approx\frac{1}{2c}R_{i}\bm{\sigma}\cdot\bm{p}\psi_{i}^{L},$
(68) $\displaystyle R_{i}(\bm{r})$
$\displaystyle=[1+\frac{1}{2c^{2}}(\epsilon_{i}-V^{SS}_{\mathrm{eff}}(\bm{r}))]^{-1}\stackrel{{\scriptstyle
c\rightarrow\infty}}{{\rightarrow}}1,$ (69)

where the potential terms in $G^{SL}$ (23) and $G^{SS}$ (24) have been grouped
into $V^{SL}_{\mathrm{eff}}$ and $V^{SS}_{\mathrm{eff}}$, respectively, when
going from Eq. (66) to Eq. (67). The former is further neglected when going
from Eq. (67) to Eq. (68), for it is of $\mathcal{O}(c^{-1})$. The major
effect of $\bm{\sigma}\cdot\bm{p}$ is to change the parity of the large
component to that of the small component. So the $S/L$ ratio is determined
mainly by the $R_{i}(\bm{r})$ operator (69), which is extremely short
ranged[90]: Each $R_{i}(\bm{r})$ becomes a constant factor just slightly away
from the position of a heavy atom (e.g., $0.05$ a.u. (roughly the radii of
$2s$ and $2p$ shells) in the case of Rn; cf. Fig. 1 in Ref. 45). Imagine we
have first solved the (radial) Dirac equation for each isolated (spherical and
unpolarized) atom $A$ and thus obtained the corresponding A4Ss
$\\{\varphi_{Ap}\\}$. Then, the atoms are brought together to synthesize the
molecule. While both the large and small components of $\varphi_{Ap}$ will
change, the $S/L$ ratio will _not_![45, 59]. The mathematical realization[54,
55, 56, 57, 58, 59] of such a physical picture is to expand the positive-
energy M4Ss $\psi_{i}$ in terms only of the _positive-energy_ A4Ss
$\\{\varphi_{Ap}\\}$, viz.,

$\displaystyle|\psi_{i}\rangle$ $\displaystyle=\sum_{A}^{N_{A}}\sum_{p\in
A}|\varphi_{Ap}\rangle\bar{C}_{+,Ap,i}=\sum_{A}^{N_{A}}\sum_{p\in
A}\begin{pmatrix}|\varphi_{Ap}^{L}\rangle\\\
|\varphi_{Ap}^{S}\rangle\end{pmatrix}\bar{C}_{+,Ap,i}$ (70)
$\displaystyle=\sum_{A}^{N_{A}}\sum_{p\in A}\sum_{\mu\in
A}\begin{pmatrix}|\chi^{L}_{A\mu}\rangle a_{A\mu,Ap}\\\
|\chi^{S}_{A\mu}\rangle
b_{A\mu,Ap}\end{pmatrix}\bar{C}_{+,Ap,i}=\sum_{A}^{N_{A}}\sum_{\mu\in
A}\begin{pmatrix}|\chi^{L}_{A\mu}\rangle\bar{A}_{A\mu,i}\\\
|\chi^{S}_{A\mu}\rangle\bar{B}_{A\mu,i}\end{pmatrix},$ (71)

which gives rise to the following projected four-component (P4C) approach

$\displaystyle\mathbf{F}_{+}^{\mathrm{P4C}}\bar{\mathbf{C}}_{+}$
$\displaystyle=\bar{\mathbf{M}}_{+}\bar{\mathbf{C}}_{+}\bar{\mathbf{E}}_{+},$
(72) $\displaystyle\mathbf{F}_{+}^{\mathrm{P4C}}$
$\displaystyle=\bar{\mathbf{h}}^{\mathrm{P4C}}+\bar{\mathbf{G}}^{\mathrm{P4C}}[\bar{\mathbf{P}}],$
(73) $\displaystyle\bar{\mathbf{h}}^{\mathrm{P4C}}$
$\displaystyle=\sum_{X,Y=L,S}\langle[\varphi^{X}]|h^{XY}|[\varphi^{Y}]\rangle=\sum_{X,Y=L,S}\bar{\mathbf{h}}^{XY}$
(74)
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{h}^{LL}\mathbf{a}+\mathbf{a}^{\dagger}\mathbf{h}^{LS}\mathbf{b}+\mathbf{b}^{\dagger}\mathbf{h}^{SL}\mathbf{a}+\mathbf{b}^{\dagger}\mathbf{h}^{SS}\mathbf{b},$
(75) $\displaystyle\bar{\mathbf{G}}^{\mathrm{P4C}}[\bar{\mathbf{P}}]$
$\displaystyle=\sum_{X,Y=L,S}\langle[\varphi^{X}]|G^{XY}|[\varphi^{Y}]\rangle=\sum_{X,Y=L,S}\bar{\mathbf{G}}^{XY}[\bar{\mathbf{D}}]$
(76)
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{G}^{LL}[\bar{\mathbf{D}}]\mathbf{a}+\mathbf{a}^{\dagger}\mathbf{G}^{LS}[\bar{\mathbf{D}}]\mathbf{b}+\mathbf{b}^{\dagger}\mathbf{G}^{SL}[\bar{\mathbf{D}}]\mathbf{a}+\mathbf{b}^{\dagger}\mathbf{G}^{SS}[\bar{\mathbf{D}}]\mathbf{b},$
(77) $\displaystyle\mathbf{a}$
$\displaystyle=\sum_{A}^{\oplus}\mathbf{a}_{A},\quad\mathbf{b}=\sum_{A}^{\oplus}\mathbf{b}_{A},$
(78) $\displaystyle\bar{\mathbf{D}}$
$\displaystyle=\begin{pmatrix}\bar{\mathbf{D}}^{LL}&\bar{\mathbf{D}}^{LS}\\\
\bar{\mathbf{D}}^{SL}&\bar{\mathbf{D}}^{SS}\end{pmatrix}=\begin{pmatrix}\mathbf{a}\bar{\mathbf{P}}\mathbf{a}^{\dagger}&\mathbf{a}\bar{\mathbf{P}}\mathbf{b}^{\dagger}\\\
\mathbf{b}\bar{\mathbf{P}}\mathbf{a}^{\dagger}&\mathbf{b}\bar{\mathbf{P}}\mathbf{b}^{\dagger}\end{pmatrix},\quad\bar{\mathbf{P}}=\bar{\mathbf{C}}_{+}\mathbf{n}\bar{\mathbf{C}}_{+}^{\dagger},$
(79) $\displaystyle\bar{\mathbf{M}}_{+}$
$\displaystyle=\bar{\mathbf{M}}^{LL}+\bar{\mathbf{M}}^{SS}=\mathbf{a}^{\dagger}\mathbf{M}^{LL}\mathbf{a}+\mathbf{b}^{\dagger}\mathbf{M}^{SS}\mathbf{b}.$
(80)

The dimension of $\mathbf{F}_{+}^{\mathrm{P4C}}$ (73) in the $j$-adapted
spinor basis $\\{\varphi_{p}\\}_{p=1}^{2n}$ is $2n$ instead of $4n$. That is,
molecular NESs are excluded completely. What is neglected here is rotations
between the PESs of an atom and the NESs of the other atoms in the molecule, a
kind of polarization of the atomic vacua. Being of $\mathcal{O}(c^{-4})$, such
an approximation introduces no discernible errors to molecular spectra [58,
44, 45]. Yet, it can be envisaged that P4C will break down when the
interatomic distance between two heavy atoms is very short (which may occur in
highly charged molecular systems). In this case, the A4Ss may be replaced with
diatomic 4-spinors (DA4S) obtained by diagonalizing

$\displaystyle\mathbf{F}^{\mathrm{Frag}}_{\mu\nu}=h^{\mathrm{Frag}}_{\mu\nu}+G^{\mathrm{Frag}}_{\mu\nu}[\sum_{A\in
F}^{\oplus}\mathbf{D}_{A}],\quad\mu,\nu\in F,$ (81)

for every pair $F$ of atoms. Here, $\mathbf{D}_{A}$ is the four-component
density matrix of atom $A$. In essence, it is the interatomic interaction
strength that is taken here as a perturbation parameter to expand the
projector for the molecular no-pair relativistic Hamiltonian[45, 91].

The structures of $\bar{\mathbf{h}}^{\mathrm{P4C}}$ (75),
$\bar{\mathbf{G}}^{\mathrm{P4C}}$ (77), and $\bar{\mathbf{M}}_{+}$ (80) are
reminiscent[12] of the NESC (normalized elimination of the small component)
approach[38], which becomes more transparent by introducing the following
formal relations

$\displaystyle\mathbf{b}_{A}$
$\displaystyle=\mathbf{X}_{A}\mathbf{a}_{A},\quad\mathbf{b}=\mathbf{X}\mathbf{a},\quad\mathbf{X}=\sum_{A}^{\oplus}\mathbf{X}_{A},$
(82)

such that

$\displaystyle\bar{\mathbf{h}}^{\mathrm{P4C}}$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,1e}}\mathbf{a},$
(83) $\displaystyle\bar{\mathbf{G}}^{\mathrm{P4C}}[\bar{\mathbf{P}}]$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\bar{\mathbf{D}}]\mathbf{a},$
(84) $\displaystyle\mathbf{L}_{+}^{\mathrm{NESC,1e}}$
$\displaystyle=\mathbf{V}+\mathbf{T}\mathbf{X}+\mathbf{X}^{\dagger}\mathbf{T}+\frac{1}{4c^{2}}\mathbf{X}^{\dagger}\mathbf{W}\mathbf{X}-\mathbf{X}^{\dagger}\mathbf{T}\mathbf{X},$
(85) $\displaystyle\mathbf{L}^{\mathrm{NESC,2e}}_{+}[\mathbf{Z}]$
$\displaystyle=\mathbf{G}^{LL}[\mathbf{Z}]+\mathbf{G}^{LS}[\mathbf{Z}]\mathbf{X}+\mathbf{X}^{\dagger}\mathbf{G}^{SL}[\mathbf{Z}]+\mathbf{X}^{\dagger}\mathbf{G}^{SS}[\mathbf{Z}]\mathbf{X}$
(86)
$\displaystyle\triangleq\mathbf{G}^{\mathrm{NR}}[\mathbf{Z}^{LL}]+\mathbf{G}^{\mathrm{NESC}}[\mathbf{Z}],$
(87) $\displaystyle G^{\mathrm{NR}}_{\mu\nu}[\mathbf{Z}]$
$\displaystyle=(\Omega^{LL}_{\mu\nu}|g_{12}|\mathrm{Tr}_{2}[\mathbf{Z}_{\lambda\kappa}\Omega^{LL}_{\kappa\lambda}])-(\Omega^{LL}_{\mu\lambda}|g_{12}|\mathbf{Z}_{\lambda\kappa}\Omega^{LL}_{\kappa\nu}),$
(88) $\displaystyle\bar{\mathbf{M}}_{+}$
$\displaystyle=\mathbf{a}^{\dagger}\tilde{\mathbf{M}}_{+}\mathbf{a},\quad\tilde{\mathbf{M}}_{+}=\mathbf{M}^{LL}+\mathbf{X}^{\dagger}\mathbf{M}^{SS}\mathbf{X}.$
(89)

It looks like that P4C is just NESC working with the A4S basis and meanwhile
the atomic approximation (AtomX) (82) to the molecular decoupling (transfer)
matrix $\mathbf{X}$,

$\displaystyle\mathbf{C}_{+}^{S}=\mathbf{X}\mathbf{C}_{+}^{L}.$ (90)

However, P4C and NESC were introduced in completely different ways. The former
is by construction a direct four-component approximation to the DHF equation
(31), whereas the latter is a two-component theory for the large components of
the PESs, viz.,

$\displaystyle\mathbf{L}^{\mathrm{NESC}}_{+}\mathbf{C}^{L}_{+}$
$\displaystyle=\tilde{\mathbf{M}}_{+}\mathbf{C}^{L}_{+}\mathbf{E}_{+},\quad\mathbf{L}^{\mathrm{NESC}}_{+}=\mathbf{L}^{\mathrm{NESC,1e}}_{+}+\mathbf{L}^{\mathrm{NESC,2e}}_{+}.$
(91)

The AtomX (82) is not mandatory but can be introduced[59] to
$\mathbf{L}^{\mathrm{NESC}}_{+}$ later on. On the practical side, it is not
mandatory for NESC to use a generally contracted basis, even with the AtomX
(e.g., $\mathbf{a}_{A}$ and $\mathbf{b}_{A}$ can be set to unit matrices for
an uncontracted RKB basis). On the other hand, P4C need not be limited to a
RKB basis. Rather, it can also adopt, e.g., numerical A4Ss by solving the
atomic radial DHF equations with grids (cf. the second equality of Eq. (70)
and Eqs. (74), and (76)). A more important distinction between P4C and NESC
lies in that the eigenvectors of the P4C eigenequation (72) give the total
density represented in the A4S basis $\\{\varphi_{Ap}\\}$, whereas those of
the NESC eigenequation (91) give only the large-component density represented
in the basis $\\{\chi^{L}_{\mu}\\}$.

At this stage, the computational efficiency is gained only in the matrix
diagonalization step (by a factor of 8 in the case of full diagonalization or
by a factor of 4 in the case of partial diagonalization by iVI[77, 78]), which
is very little for a moderate basis. The real gain in efficiency can only be
achieved by further invoking a ‘model small component approximation’ (MSCA),
which was proposed[44] originally in the context of density functional theory
(DFT). Therein, the molecular small-component density $\rho^{S}$ is first
approximated as the superposition of the atomic ones (cf. Eq. (62)),

$\displaystyle\rho^{S}=\mathrm{Tr}[\Omega^{SS}\bar{\mathbf{P}}]\approx\tilde{\rho}^{S}=\sum_{A}\tilde{\rho}^{S}_{A}=\mathrm{Tr}[\Omega^{SS}\bar{\mathbf{P}}^{0}],\quad\Omega^{SS}_{pq}=\varphi_{p}^{S{\dagger}}\varphi_{q}^{S},\quad\bar{\mathbf{P}}^{0}=\sum_{A}^{\oplus}\mathbf{n}_{A},$
(92)

such that the molecular density $\rho^{4c}$ can be calculated as

$\displaystyle\rho^{4c}$
$\displaystyle\approx\tilde{\rho}^{4c}=\rho^{L}+\tilde{\rho}^{S},\quad\rho^{L}=\mathrm{Tr}[\Omega^{LL}\bar{\mathbf{P}}],\quad\Omega^{LL}_{pq}=\varphi_{p}^{L{\dagger}}\varphi_{q}^{L}.$
(93)

The matrix elements of the local Kohn-Sham (KS) potential
$V_{\mathrm{eff}}[\rho^{4c}]$ are then approximated as

$\displaystyle\langle\varphi_{p}^{L}|V_{\mathrm{eff}}[\rho^{4c}]|\varphi_{q}^{L}\rangle+\langle\varphi_{p}^{S}|V_{\mathrm{eff}}[\rho^{4c}]|\varphi_{q}^{S}\rangle$
$\displaystyle\approx\langle\varphi_{p}^{L}|V_{\mathrm{eff}}[\tilde{\rho}^{4c}]|\varphi_{q}^{L}\rangle+\langle\varphi_{p}^{S}|V_{\mathrm{eff}}[\rho^{4c}_{\mathrm{mod}}]|\varphi_{q}^{S}\rangle,$
(94) $\displaystyle\rho^{4c}_{\mathrm{mod}}$
$\displaystyle=\sum_{A}\rho^{4c}_{A}=\mathrm{Tr}[(\Omega^{LL}+\Omega^{SS})\bar{\mathbf{P}}^{0}],$
(95)

such that, under the LCA4S and MSCA, only the first iteration of the
calculation is four-component but subsequent iterations are just two-component
(i.e., only $\rho^{L}$ need to be updated), thereby justifying the name quasi-
four-component (Q4C)[44]. Even the first iteration is much cheaper than a
regular four-component iteration, for the sparsity associated with
$\bar{\mathbf{P}}^{0}$ renders the evaluation of the two-electron matrix
elements only at a fractional cost of a regular evaluation of the two-electron
matrix elements. It has been shown[44, 61] that Q4C-KS is indeed very
accurate, e.g., with errors being only a few milli-Hartrees for the energy
levels of the innermost shells of E1172 and completely negligible for
molecular spectroscopic constants. The MSCA[45], which can also be termed
‘model density approximation’, can trivially be generalized to a model density
matrix (MDM) approximation to incorporate HF exchange. That is, the molecular
density matrix elements $\bar{\mathbf{D}}_{\lambda\kappa}$ (cf. Eq. (79)) in
all small-component-containing terms of
$\mathbf{L}^{\mathrm{NESC,2e}}_{+}[\bar{\mathbf{D}}]$ (cf. Eq. (87) and
$\mathbf{G}^{XY}[\bar{\mathbf{D}}]$ in Eqs. (51), (53), (55), and (57)) are
replaced with those of the model density matrix
$\bar{\mathbf{D}}^{0}=\sum_{A}^{\oplus}\bar{\mathbf{D}}_{A}$. This amounts to
replacing $\mathbf{G}^{\mathrm{NESC}}[\bar{\mathbf{D}}]$ in Eq. (87) with
$\mathbf{G}^{\mathrm{NESC}}[\bar{\mathbf{D}}^{0}]$, thereby giving rise to the
Q4C-HF equation[60, 61] represented in the RKB basis employed to expand the
A4Ss $\\{\varphi_{Ap}\\}$,

$\displaystyle{\mathbf{F}}_{+}^{\mathrm{Q4C}}\bar{\mathbf{C}}_{+}$
$\displaystyle=\bar{\mathbf{M}}_{+}\bar{\mathbf{C}}_{+}\bar{\mathbf{E}}_{+},$
(96) $\displaystyle{\mathbf{F}}_{+}^{\mathrm{Q4C}}$
$\displaystyle=\bar{\mathbf{h}}^{\mathrm{Q4C}}_{\mathrm{eff}}+\bar{\mathbf{G}}^{\mathrm{NR}}[\bar{\mathbf{P}}],$
(97) $\displaystyle\bar{\mathbf{h}}^{\mathrm{Q4C}}_{\mathrm{eff}}$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,1e}}\mathbf{a}+\bar{\mathbf{G}}^{\mathrm{Q4C}}_{\mathrm{eff}},\quad\bar{\mathbf{G}}^{\mathrm{Q4C}}_{\mathrm{eff}}=\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NESC}}[\bar{\mathbf{D}}^{0}]\mathbf{a},$
(98) $\displaystyle\bar{\mathbf{G}}^{\mathrm{NR}}[\bar{\mathbf{P}}]$
$\displaystyle=\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\bar{\mathbf{D}}^{LL}]\mathbf{a},\quad\bar{\mathbf{D}}^{LL}=\mathbf{a}^{\dagger}\bar{\mathbf{P}}\mathbf{a},$
(99)

where $\bar{\mathbf{G}}^{\mathrm{Q4C}}_{\mathrm{eff}}$ can be interpreted as a
correction of two-electron picture-change errors (2ePCE)[92], although
everything is done here within the four-component framework. Its construction
is very cheap due to the sparsity associated with $\bar{\mathbf{D}}^{0}$.

Usually only the occupied and low-lying virtual A4Ss
$\\{\varphi_{p}\\}_{p=1}^{2n_{B}}$ are needed to form the backbone of the
basis, which is to be augmented with some uncontracted flat functions
$\\{f_{\nu}\\}_{\nu=2n_{B}+1}^{2n}$ for describing the deformation and
polarization of the atoms when forming the molecule. As
$R_{\nu}(\vec{r})\approx 1$ in the valence region (cf. Fig. 1 in Ref. 45),
such flat functions can simply be taken as

$\displaystyle f_{\nu}=\begin{pmatrix}\chi^{L}_{\nu}\\\
\chi_{\nu}^{S}\end{pmatrix},\quad\chi_{\nu}^{S}=\frac{1}{2c}\Pi\chi^{L}_{\nu}\approx
0,\quad\nu\in[2n_{B}+1,2n].$ (100)

Unlike Eq. (27), the large and small component basis functions are here
combined together. Such a ‘HF+P’ type of spinor basis, single-zeta (or double-
zeta) for core shells and multiple-zeta for valence shells, is very
effective[45, 44]. The elements of $\mathbf{F}^{\mathrm{Q4C}}_{+}$ among the
added flat functions $\\{f_{\nu}\\}_{\nu=2n_{B}+1}^{2n}$ can be treated
nonrelativistically,

$\displaystyle F^{\mathrm{Q4C}}_{+,\mu\nu}$
$\displaystyle=\langle\chi^{L}_{\mu}|T|\chi^{L}_{\nu}\rangle+\langle\chi^{L}_{\mu}|V|\chi^{L}_{\nu}\rangle+G^{\mathrm{NR}}_{\mu\nu}[\bar{\mathbf{P}}],\quad\mu,\nu\in[2n_{B}+1,2n],$
(101)

while those between $\\{\varphi_{p}\\}_{p=1}^{2n_{B}}$ and
$\\{f_{\nu}\\}_{\nu=2n_{B}+1}^{2n}$ can be treated in the same way as Eq. (97)
or simply approximated as

$\displaystyle F^{\mathrm{Q4C}}_{+,p\nu}$
$\displaystyle=\bar{h}^{\mathrm{Q4C}}_{\mathrm{eff},p\nu}+\bar{G}^{\mathrm{NR}}_{p\nu}[\bar{\mathbf{P}}]=[F^{\mathrm{Q4C}}_{+,\nu
p}]^{*},\quad p\in[1,2n_{B}],\quad\nu\in[2n_{B}+1,2n],$ (102)
$\displaystyle\bar{G}^{\mathrm{NR}}_{p\nu}[\bar{\mathbf{P}}]$
$\displaystyle=(\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\bar{\mathbf{D}}^{LL}])_{p\nu},$
(103) $\displaystyle\bar{h}^{\mathrm{Q4C}}_{\mathrm{eff},p\nu}$
$\displaystyle=\langle\varphi_{p}^{L}|T|\chi^{L}_{\nu}\rangle+\langle\varphi_{p}^{L}|V|\chi^{L}_{\nu}\rangle+\bar{G}^{\mathrm{Q4C}}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$
(104)
$\displaystyle=[\mathbf{a}^{\dagger}(\mathbf{T}+\mathbf{V}+\mathbf{G}^{\mathrm{NESC}}[\bar{\mathbf{D}}^{0}])]_{p\nu},$
(105) $\displaystyle G^{\mathrm{NESC}}_{\mu\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=G^{\mathrm{NESC,LL}}_{\mu\nu}[\bar{\mathbf{D}}^{0}]+G^{\mathrm{NESC},SL}_{\mu\nu}[\bar{\mathbf{D}}^{0}],$
(106)
$\displaystyle\bar{G}^{\mathrm{Q4C}}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=\bar{G}^{\mathrm{Q4C},LL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]+\bar{G}^{\mathrm{Q4C},SL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]=(\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NESC}}[\bar{\mathbf{D}}^{0}])_{p\nu},$
(107)
$\displaystyle\bar{G}^{\mathrm{Q4C},LL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=(\Omega_{p\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\bar{P}^{0}_{rs}\Omega^{SS}_{sr}])-\mathrm{c_{g}}(\Xi_{pr}^{LS}|g_{12}\cdot|\bar{P}^{0}_{rs}\Xi_{s\nu}^{SL})$
$\displaystyle-\mathrm{c_{b}}(\Xi_{pr}^{LS}|\cdot\mathbf{b}_{12}\cdot|\bar{P}^{0}_{rs}\Xi_{s\nu}^{SL})$
(108)
$\displaystyle=(\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NESC},LL}[\bar{\mathbf{D}}^{0}])_{p\nu},$
(109)
$\displaystyle\mathbf{G}^{\mathrm{NESC},LL}_{\mu\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=(\Omega_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\bar{D}^{0,SS}_{\lambda\kappa}\Omega^{SS}_{\kappa\lambda}])-\mathrm{c_{g}}(\Xi_{\mu\lambda}^{LS}|g_{12}\cdot|\bar{D}^{0,SS}_{\lambda\kappa}\Xi_{\kappa\nu}^{SL})$
$\displaystyle-\mathrm{c_{b}}(\Xi_{\mu\lambda}^{LS}|\cdot\mathbf{b}_{12}\cdot|\bar{D}^{0,SS}_{\lambda\kappa}\Xi_{\kappa\nu}^{SL}),$
(110)
$\displaystyle\bar{G}^{\mathrm{Q4C},SL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=-(\Omega_{pr}^{SS}|g_{12}|\bar{P}^{0}_{rs}\Omega^{LL}_{s\nu})$
$\displaystyle+\mathrm{c_{g}}[(\Xi_{p\nu}^{SL}|g_{12}\cdot|\mathrm{Tr}_{2}[\bar{P}^{0}_{rs}(\Xi_{sr}^{SL}+\Xi_{sr}^{LS})])-(\Xi_{pr}^{SL}|g_{12}\cdot|\bar{P}^{0}_{rs}\Xi_{s\nu}^{SL})]$
$\displaystyle+\mathrm{c_{b}}[(\Xi_{p\nu}^{SL}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\bar{P}^{0}_{rs}(\Xi_{sr}^{SL}+\Xi_{sr}^{LS})])-(\Xi_{pr}^{SL}|\cdot\mathbf{b}_{12}\cdot|\bar{P}^{0}_{rs}\Xi_{s\nu}^{SL})]$
(111)
$\displaystyle=(\mathbf{a}^{\dagger}\mathbf{G}^{\mathrm{NESC},SL}[\bar{\mathbf{D}}^{0}])_{p\nu},$
(112)
$\displaystyle\mathbf{G}^{\mathrm{NESC},SL}_{\sigma\nu}[\bar{\mathbf{D}}^{0}]$
$\displaystyle=X_{\mu\sigma}^{*}\\{-(\Omega_{\mu\lambda}^{SS}|g_{12}|\bar{D}^{0,SL}_{\lambda\kappa}\Omega^{LL}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[(\Xi_{\mu\nu}^{SL}|g_{12}\cdot|\mathrm{Tr}_{2}[\bar{D}^{0,LS}_{\lambda\kappa}\Xi_{\kappa\lambda}^{SL}+\bar{D}^{0,SL}_{\lambda\kappa}\Xi_{\kappa\lambda}^{LS}])-(\Xi_{\mu\lambda}^{SL}|g_{12}\cdot|\bar{D}^{0,LS}_{\lambda\kappa}\Xi_{\kappa\nu}^{SL})]$
$\displaystyle+\mathrm{c_{b}}[(\Xi_{\mu\nu}^{SL}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\bar{D}^{0,LS}_{\lambda\kappa}\Xi_{\kappa\lambda}^{SL}+\bar{D}^{0,SL}_{\lambda\kappa}\Xi_{\kappa\lambda}^{LS}])-(\Xi_{\mu\lambda}^{SL}|\cdot\mathbf{b}_{12}\cdot|\bar{D}^{0,LS}_{\lambda\kappa}\Xi_{\kappa\nu}^{SL})]\\}.$
(113)

The corresponding overlap matrix elements read

$\displaystyle\bar{M}_{+,\mu\nu}$
$\displaystyle=M^{LL}_{\mu\nu}+M^{SS}_{\mu\nu},\quad\mu,\nu\in[2n_{B}+1,2n],$
(114) $\displaystyle\bar{M}_{+,p\nu}$
$\displaystyle=[\mathbf{a}^{\dagger}(\mathbf{M}^{LL}+\mathbf{M}^{SS})]_{p\nu},\quad
p\in[1,2n_{B}].$ (115)

The ‘rule’ for obtaining the elements
$\bar{G}^{\mathrm{Q4C},XY}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$ is to
set the $\chi^{S}_{\nu}$-containing terms of $G^{XY}_{p\nu}$ (cf. Eqs. (50),
(52), (54), and (56)) to zero, thereby leaving only
$\bar{G}^{\mathrm{Q4C},LL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$ (108)
and $\bar{G}^{\mathrm{Q4C},SL}_{\mathrm{eff},p\nu}[\bar{\mathbf{D}}^{0}]$
(111) stemming from $G^{LL}_{p\nu}$ (50) and $G^{SL}_{p\nu}$ (54),
respectively. The use of the same expansion coefficients
$\\{\bar{\mathbf{C}}_{+,i}\\}$ for the small and large components of the M4Ss
$\\{\psi_{i}\\}$ is also made here (cf. Eq. (70)). Eq. (100) along with the
expressions (101) and (102) has been termed RKB0[45], which is very accurate
for spectroscopic constants even of the heaviest molecular systems[93, 94, 95,
96, 97].

Additional remarks on Q4C should be made here:

  1. (a)

The above ‘HF+P’ type of spinor basis, subject to the RKB and RKB0 conditions,
is most compact.

  2. (b)

Once the A4Ss are symmetrized according to both double point group and time
reversal symmetries[82], the quaternion form of Q4C-HF can readily be obtained
by the $\mbox{}^{\mathrm{Q}}\mathbf{U}$ transformation (295) of
$\mathbf{F}^{\mathrm{Q4C}}_{+}$ (cf. Eqs. (97), (101), and (102)). The Q4C
variant of KROHF can then be formulated (see Appendix A).

  3. (c)

Since the large and small components of $\varphi_{p}$ share the same
coefficients, the A4S integrals
$(\varphi_{p}\varphi_{q}|V(1,2)|\varphi_{r}\varphi_{s})$ can be transformed to
the MO representation as a whole (instead of component-wise). Therefore, both
the integral transformation and correlation steps in Q4C are computationally
the same as those in two-component approaches[60, 61].

  4. (d)

Unlike two-component approaches, Q4C does not suffer from 2ePCEs in both the
mean-field and correlation steps.

  5. (e)

The MDM approximation can also be applied to the DHF Fock matrix (43), leading
to

$\displaystyle\mathbf{F}$
$\displaystyle=\mathbf{h}_{\mathrm{eff}}^{\mathrm{4C}}+\mathbf{G}^{\mathrm{4C}}[\mathbf{D}],$
(116) $\displaystyle\mathbf{h}_{\mathrm{eff}}^{\mathrm{4C}}$
$\displaystyle=\begin{pmatrix}\mathbf{h}^{LL}+\mathbf{G}^{\prime
LL}[\mathbf{D}^{0}]&\mathbf{h}^{LS}+\mathbf{G}^{LS}[\mathbf{D}^{0}]\\\
\mathbf{h}^{SL}+\mathbf{G}^{SL}[\mathbf{D}^{0}]&\mathbf{h}^{SS}+\mathbf{G}^{SS}[\mathbf{D}^{0}]\end{pmatrix},\quad\mathbf{D}^{0}=\sum_{A}^{\oplus}\mathbf{D}_{A},$
(117) $\displaystyle\mathbf{G}^{\prime LL}[\mathbf{D}^{0}]$
$\displaystyle=\mathbf{G}^{LL}[\mathbf{D}^{0}]-\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{0,LL}],$
(118) $\displaystyle\mathbf{G}^{\mathrm{4C}}[\mathbf{D}]$
$\displaystyle=\begin{pmatrix}\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}\end{pmatrix},$ (119)

where only $\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]$ needs to be updated in
each iteration.

##  IV X2C

By definition, a two-component relativistic theory is a Schrödinger-like
equation that describes only electrons relativistically. In this sense,
NESC[38] is not yet a genuine two-component relativistic theory, since the
eigenequation (91), along with the relativistic metric (89), determines only
the large components of the M4Ss, whereas the small components have to be
constructed explicitly via relation (90). That is, NESC is still within the
Dirac picture, _unlike_ P4C and Q4C. The real significance of NESC lies in the
initiative of matrix formulations of two-component relativistic theories. That
is, the starting point is the matrix DHF equation (43) instead of the operator
DHF equation (25). The former can be block-diagonalized exactly, so as to
achieve exact decoupling of the PESs from the NESs, ending up with an “exact
two-component” (X2C)[98] relativistic theory for electrons. This is different
from approximate two-component (A2C) theories[39, 40, 41, 42] that are correct
only to a finite order in relativity. The block-diagonalization of Eq. (43)
can be done in one step[50, 51, 99, 43], two steps[100, 101, 102], and
multiple steps[103, 104, 105]. It has been shown[60] that the three types of
matrix formulations share the same decoupling condition and differ only in the
renormalization. There exist even closed mappings in between[60]. Since the
initio free-particle transformation invoked in the two-step and multiple-step
formulations is only necessary to ensure variational stability and
regularization of finite-order expansions[41, 42] but not needed for an exact
decoupling, and becomes very clumsy in the presence of magnetic interactions,
it is clear that it is the one-step formulation[50, 51, 99, 43] that should be
advocated. It deserves to be mentioned that the same one-step matrix X2C
equation[50, 51] can also be obtained by the (formal) Foldy-Wouthuysen (FW)
transformation[106] of the operator DHF equation (25), followed by making use
of the RKB-RI to convert the operators therein to matrix forms[60]. That is,
the matrix and operator (more precisely, operator-like) formulations of the
one-step X2C are actually identical. This applies also to the two-step[100,
101, 102] and multiple-step[103, 104, 105] formulations: instead of the usual
operator formulations followed by the use of the RKB-RI, the same matrix
equations therein can also be obtained[107] by starting with the partitioned
matrix DHF equation (43). In essence, the use of the RKB-RI in the
_unexpanded_ one-step, two-step, and multiple-step operator formulations is
merely a formal step and does not introduce any error even with a finite RKB
basis. This is different from A2C approaches, where the use of the RKB-RI for
their matrix representations is indeed an approximation [12].

The one-step block-diagonalization of the matrix DHF equation (43) starts with
the formal relation (90) between the small- and large-component coefficients
of the PESs, and a similar relation

$\displaystyle\mathbf{C}^{L}_{-}=\tilde{\mathbf{X}}\mathbf{C}^{S}_{-}$ (120)

between the small (upper)- and large (lower)-component coefficients of the
NESs. In terms of such relations, the following transformation matrix
$\mathbf{U}$ can be constructed[60]

$\displaystyle\mathbf{U}$
$\displaystyle=\begin{pmatrix}\mathbf{U}^{LL}&\mathbf{U}^{LS}\\\
\mathbf{U}^{SL}&\mathbf{U}^{SS}\end{pmatrix}=\bm{\Omega}_{D}\bm{\Omega}_{N}=\begin{pmatrix}\mathbf{R}_{+}&\tilde{\mathbf{X}}\mathbf{R}_{-}\\\
\mathbf{X}\mathbf{R}_{+}&\mathbf{R}_{-}\end{pmatrix},$ (121)
$\displaystyle\bm{\Omega}_{D}$
$\displaystyle=\begin{pmatrix}\mathbf{I}&\tilde{\mathbf{X}}\\\
\mathbf{X}&\mathbf{I}\end{pmatrix},\quad\bm{\Omega}_{N}=\begin{pmatrix}\mathbf{R}_{+}&\mathbf{0}\\\
\mathbf{0}&\mathbf{R}_{-}\end{pmatrix},$ (122)

where $\bm{\Omega}_{D}$ does the decoupling, whereas $\bm{\Omega}_{N}$
establishes the renormalization. Note that the matrix $\mathbf{U}$ defined
here is just the matrix representation[51] of the FW transformation[106] in a
finite RKB basis. The decoupling matrix $\mathbf{X}$ (cf. Eq. (90)) is to be
determined by the condition

$\displaystyle(\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U})^{SL}=\sum_{X,Y=L,S}(\mathbf{U}^{\dagger})^{SX}\mathbf{F}^{XY}\mathbf{U}^{YL}=\mathbf{0}.$
(123)

More specifically[50],

$\displaystyle\mathbf{F}^{SL}+\mathbf{F}^{SS}\mathbf{X}$
$\displaystyle=-\tilde{\mathbf{X}}^{\dagger}\mathbf{L}_{+}^{\mathrm{UESC}},\quad\mathbf{L}_{+}^{\mathrm{UESC}}=\mathbf{F}^{LL}+\mathbf{F}^{LS}\mathbf{X},$
(124)

where $\mathbf{L}_{+}^{\mathrm{UESC}}$ is the UESC (unnormalized elimination
of the small component) Hamiltonian[38] associated with the following
eigenvalue problem,

$\displaystyle\mathbf{L}_{+}^{\mathrm{UESC}}\mathbf{C}^{L}_{+}$
$\displaystyle=\mathbf{M}^{LL}\mathbf{C}^{L}_{+}\mathbf{E}_{+}.$ (125)

Further combined with Eq. (91), we have[108]

$\displaystyle\mathbf{L}_{+}^{\mathrm{UESC}}=\mathbf{M}^{LL}\tilde{\mathbf{M}}^{-1}_{+}\mathbf{L}_{+}^{\mathrm{NESC}}.$
(126)

Eq. (124) can hence be rewritten as

$\displaystyle\mathbf{F}^{SL}+\mathbf{F}^{SS}\mathbf{X}$
$\displaystyle=-\tilde{\mathbf{X}}\mathbf{M}^{LL}\tilde{\mathbf{M}}^{-1}_{+}\mathbf{L}_{+}^{\mathrm{NESC}},$
(127)

which turns out to be more robust than Eq. (124) when solved iteratively[99].
It deserves to be emphasized that both Eqs. (124) and (127), proposed in Refs.
50 and 99, respectively, in different forms though, are state-independent and
hence fundamentally different from the state-dependent decoupling conditions
proposed by Dyall[38]. The orthogonality condition
$\mathbf{C}_{+}^{\dagger}\mathbf{M}\mathbf{C}_{-}=\mathbf{0}$ leads to

$\tilde{\mathbf{X}}=-(\mathbf{M}^{LL})^{-1}\mathbf{X}^{\dagger}\mathbf{M}^{SS}.$
(128)

It follows that $\tilde{\mathbf{X}}$ is determined directly by $\mathbf{X}$.
In other words, the two sets solutions, $\mathbf{X}$ and $\tilde{\mathbf{X}}$,
of the quadratic decoupling condition (124)/(127) are mutually related.
Without going into further details, the renormalization matrices read[51]

$\displaystyle\mathbf{R}_{+}$
$\displaystyle=[(\mathbf{M}^{LL})^{-1}\tilde{\mathbf{M}}_{+}]^{-\frac{1}{2}}=(\mathbf{M}^{LL})^{-\frac{1}{2}}\mathbf{K}_{+}(\mathbf{M}^{LL})^{\frac{1}{2}},$
$\displaystyle\mathbf{K}_{+}$
$\displaystyle=[(\mathbf{M}^{LL})^{-\frac{1}{2}}\tilde{\mathbf{M}}_{+}(\mathbf{M}^{LL})^{-\frac{1}{2}}]^{-\frac{1}{2}},$
(129) $\displaystyle\mathbf{R}_{-}$
$\displaystyle=[(\mathbf{M}^{SS})^{-1}\tilde{\mathbf{M}}_{-}]^{-\frac{1}{2}}=(\mathbf{M}^{SS})^{-\frac{1}{2}}\mathbf{K}_{-}(\mathbf{M}^{SS})^{\frac{1}{2}},$
$\displaystyle\mathbf{K}_{-}$
$\displaystyle=[(\mathbf{M}^{SS})^{-\frac{1}{2}}\tilde{\mathbf{M}}_{-}(\mathbf{M}^{SS})^{-\frac{1}{2}}]^{-\frac{1}{2}},$
(130)

where the relativistic metric for the NESs reads (cf. Eq. (89))

$\displaystyle\tilde{\mathbf{M}}_{-}$
$\displaystyle=\mathbf{M}^{SS}+\tilde{\mathbf{X}}^{\dagger}\mathbf{M}^{LL}\tilde{\mathbf{X}}.$
(131)

It has been proven[109] that the particular renormalization $\mathbf{R}_{+}$
(129) renders the resulting two-component spinors closest to the large
components of the positive-energy M4Ss in the least-squares sense, whereas the
two-step[100, 101, 102] and multiple-step[103, 104, 105] transformations
amount to adopting different renormalizations. The $\mathbf{U}$-transformation
(121) of Eq. (43) leads to

$\displaystyle(\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U})(\mathbf{U}^{-1}\mathbf{C})$
$\displaystyle=(\mathbf{U}^{\dagger}\mathbf{M}\mathbf{U})(\mathbf{U}^{-1}\mathbf{C})\mathbf{E}\Leftrightarrow\tilde{\mathbf{F}}\tilde{\mathbf{C}}=\mathbf{M}\tilde{\mathbf{C}}\mathbf{E},$
(132)

where

$\displaystyle\tilde{\mathbf{F}}$
$\displaystyle=\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U}=\begin{pmatrix}\mathbf{F}_{+}^{\mathrm{X2C}}&\mathbf{0}\\\
\mathbf{0}&\mathbf{F}_{-}^{\mathrm{X2C}}\end{pmatrix},\quad\mathbf{U}^{\dagger}\mathbf{M}\mathbf{U}=\begin{pmatrix}\mathbf{M}^{LL}&\mathbf{0}\\\
\mathbf{0}&\mathbf{M}^{SS}\end{pmatrix},$ (133)
$\displaystyle\tilde{\mathbf{C}}=\mathbf{U}^{-1}\mathbf{C}=\mathbf{M}^{-1}\mathbf{U}^{\dagger}\mathbf{M}\mathbf{C}=\begin{pmatrix}\mathbf{R}_{+}^{-1}\mathbf{C}_{+}^{L}&\mathbf{0}\\\
\mathbf{0}&\mathbf{R}_{-}^{-1}\mathbf{C}_{-}^{S}\end{pmatrix}=\begin{pmatrix}\tilde{\mathbf{C}}_{+}&\mathbf{0}\\\
\mathbf{0}&\tilde{\mathbf{C}}_{-}\end{pmatrix}.$ (134)

The upper-left block of Eq. (132) defines the unique X2C equation[50, 51] for
the PESs in the Schrödinger picture,

$\displaystyle\mathbf{F}_{+}^{\mathrm{X2C}}\tilde{\mathbf{C}}_{+}$
$\displaystyle=\mathbf{M}^{LL}\tilde{\mathbf{C}}_{+}\mathbf{E}_{+},$ (135)
$\displaystyle\mathbf{F}_{+}^{\mathrm{X2C}}$
$\displaystyle=(\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U})^{LL}=\sum_{X,Y=L,S}(\mathbf{U}^{\dagger})^{LX}\mathbf{F}^{XY}\mathbf{U}^{YL}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{X}}_{+}\mathbf{R}_{+},\quad\mathrm{X}=\mathrm{NESC,SESC},$
(136) $\displaystyle\mathbf{L}^{\mathrm{SESC}}_{+}$
$\displaystyle=\frac{1}{2}[\tilde{\mathbf{M}}_{+}(\mathbf{M}^{LL})^{-1}\mathbf{L}^{\mathrm{UESC}}_{+}+\mathrm{c.c.}].$
(137)

Here, the SESC (symmetrized elimination of the small component)[45, 99]
Hamiltonian $\mathbf{L}^{\mathrm{SESC}}_{+}$ arises from the identity
$\mathbf{L}^{\mathrm{SESC}}_{+}=\frac{1}{2}(\mathbf{L}^{\mathrm{NESC}}_{+}+\mathbf{L}^{\mathrm{NESC}}_{+})$
and the relation (126). The one- ($\mathbf{F}^{\mathrm{X2C,1e}}_{+}$) and two-
electron ($\mathbf{F}^{\mathrm{X2C,2e}}_{+}$) terms of
$\mathbf{F}_{+}^{\mathrm{X2C}}$ (136) read

$\displaystyle\mathbf{F}^{\mathrm{X2C,1e}}_{+}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{X,1e}}_{+}\mathbf{R}_{+},\quad\mathrm{X=NESC,SESC},$
(138) $\displaystyle\mathbf{L}^{\mathrm{SESC,1e}}_{+}$
$\displaystyle=\frac{1}{2}[\tilde{\mathbf{M}}_{+}(\mathbf{M}^{LL})^{-1}\mathbf{L}^{\mathrm{UESC,1e}}_{+}+\mathrm{c.c.}],\quad\mathbf{L}^{\mathrm{UESC,1e}}_{+}=\mathbf{V}+\mathbf{TX},$
(139) $\displaystyle\mathbf{F}^{\mathrm{X2C,2e}}_{+}[\mathbf{D}]$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{X,2e}}_{+}[\mathbf{D}]\mathbf{R}_{+},\quad\mathrm{X=NESC,SESC},$
(140) $\displaystyle\mathbf{L}^{\mathrm{SESC,2e}}_{+}[\mathbf{D}]$
$\displaystyle=\frac{1}{2}[\tilde{\mathbf{M}}_{+}(\mathbf{M}^{LL})^{-1}\mathbf{L}^{\mathrm{UESC,2e}}_{+}+\mathrm{c.c.}],$
(141) $\displaystyle\mathbf{L}^{\mathrm{UESC,2e}}_{+}$
$\displaystyle=\mathbf{G}^{LL}[\mathbf{D}]+\mathbf{G}^{LS}[\mathbf{D}]\mathbf{X}.$
(142)

At this stage, the two-electron term $\mathbf{F}^{\mathrm{X2C,2e}}_{+}$ (140)
still involves explicitly the four-component molecular density matrix
$\mathbf{D}$ (36), which can be eliminated as follows. In view of the relation
$\mathbf{C}=\mathbf{U}\tilde{\mathbf{C}}$ given by the first equality of Eq.
(134), i.e.,

$\displaystyle\mathbf{C}^{X}_{+}=\mathbf{U}^{XL}\tilde{\mathbf{C}}_{+},\quad
X=L,S,$ (143)

the components $\psi_{i}^{X}$ ($X=L,S$) of a positive-energy M4S $\psi_{p}$
can be expressed as

$\displaystyle\psi_{p}^{X}$
$\displaystyle=\sum_{\mu=1}^{n}\chi_{\mu}^{X}C_{+,\mu
p}^{X}=\sum_{\mu=1}^{n}\chi_{\mu}^{X}(\mathbf{U}^{XL}\tilde{\mathbf{C}}_{+})_{\mu
p}=\sum_{\nu=1}^{n}\tilde{\chi}_{\nu}^{X}\tilde{\mathbf{C}}_{+,\nu p},$ (144)
$\displaystyle\tilde{\chi}_{\nu}^{X}$
$\displaystyle=\sum_{\mu=1}^{n}\chi_{\mu}^{X}U^{XL}_{\mu\nu},\quad\tilde{\chi}_{\nu}^{X{\dagger}}=\sum_{\mu=1}^{n}(U^{XL{\dagger}})_{\nu\mu}\chi_{\mu}^{X{\dagger}}=\sum_{\mu=1}^{n}[(U^{\dagger})^{LX}]_{\nu\mu}\chi_{\mu}^{X{\dagger}},$
(145)

where $\\{\tilde{\chi}^{X}_{\nu}\\}_{\nu=1}^{n}$ ($X=L,S$) can be referred to
as renormalized two-component spinor functions. Likewise, the components
$\mathbf{D}^{XY}$ ($X,Y=L,S$) of $\mathbf{D}$ (36) can be written as

$\displaystyle\mathbf{D}^{XY}$
$\displaystyle=\mathbf{U}^{XL}\tilde{\mathbf{D}}\mathbf{U}^{YL{\dagger}}=\mathbf{U}^{XL}\tilde{\mathbf{D}}(\mathbf{U}^{\dagger})^{LY},\quad\tilde{\mathbf{D}}=\tilde{\mathbf{C}}_{+}\mathbf{n}\tilde{\mathbf{C}}^{\dagger}_{+}=\begin{pmatrix}\tilde{\mathbf{D}}^{\alpha\alpha}&\tilde{\mathbf{D}}^{\alpha\beta}\\\
\tilde{\mathbf{D}}^{\beta\alpha}&\tilde{\mathbf{D}}^{\beta\beta}\end{pmatrix}.$
(146)

It then follows that $\mathbf{F}^{\mathrm{X2C,2e}}_{+}[\mathbf{D}]$ (140) can
be expressed in terms entirely of the renormalized quantities, viz.,

$\displaystyle\mathbf{F}^{\mathrm{X2C,2e}}_{+}[\mathbf{D}]$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{X,2e}}_{+}[\tilde{\mathbf{D}}]\mathbf{R}_{+},\quad\mathrm{X=NESC,SESC},$
(147)
$\displaystyle\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{NESC,2e}}_{+}[\tilde{\mathbf{D}}]\mathbf{R}_{+}$
$\displaystyle=\tilde{G}^{LL}[\tilde{\mathbf{D}}]+\tilde{G}^{LS}[\tilde{\mathbf{D}}]+\tilde{G}^{SL}[\tilde{\mathbf{D}}]+\tilde{G}^{SS}[\tilde{\mathbf{D}}],$
(148)
$\displaystyle\mathbf{R}_{+}^{\dagger}\mathbf{L}^{\mathrm{SESC,2e}}_{+}[\tilde{\mathbf{D}}]\mathbf{R}_{+}$
$\displaystyle=\frac{1}{2}[(\bar{G}^{LL}[\tilde{\mathbf{D}}]+\bar{G}^{LS}[\tilde{\mathbf{D}}])+\mathrm{c.c.}],$
(149)

where

$\displaystyle\tilde{G}^{LL}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=(\tilde{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\lambda}])-(\tilde{\Omega}_{\mu\lambda}^{LL}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\nu})+(\tilde{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\lambda}])$
$\displaystyle-\mathrm{c_{g}}(\tilde{\bm{\Xi}}_{\mu\lambda}^{LS}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{SL})-\mathrm{c_{b}}(\tilde{\bm{\Xi}}_{\mu\lambda}^{LS}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{SL}),$
(150) $\displaystyle\tilde{G}^{LS}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=-(\tilde{\Omega}^{LL}_{\mu\lambda}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[(\tilde{\bm{\Xi}}_{\mu\nu}^{LS}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])-(\tilde{\bm{\Xi}}^{LS}_{\mu\lambda}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\nu})+(\tilde{\bm{\Xi}}^{LS}_{\mu\nu}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])]$
$\displaystyle+\mathrm{c_{b}}[\tilde{\bm{\Xi}}_{\mu\nu}^{LS}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])-(\tilde{\bm{\Xi}}^{LS}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\nu})+(\tilde{\bm{\Xi}}^{LS}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])],$
(151) $\displaystyle\tilde{G}^{SL}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=-(\tilde{\Omega}^{SS}_{\mu\lambda}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[\tilde{\bm{\Xi}}_{\mu\nu}^{SL}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])-(\tilde{\bm{\Xi}}^{SL}_{\mu\lambda}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\nu})+(\tilde{\bm{\Xi}}^{SL}_{\mu\nu}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])]$
$\displaystyle+\mathrm{c_{b}}[(\tilde{\bm{\Xi}}_{\mu\nu}^{SL}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])-(\tilde{\bm{\Xi}}^{SL}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\nu})+(\tilde{\bm{\Xi}}^{SL}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])]$
$\displaystyle=\tilde{G}^{LS*}_{\nu\mu},$ (152)
$\displaystyle\tilde{G}^{SS}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=(\tilde{\Omega}_{\mu\nu}^{SS}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\lambda}])-(\tilde{\Omega}_{\mu\lambda}^{SS}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\nu})+(\tilde{\Omega}_{\mu\nu}^{SS}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\lambda}])$
$\displaystyle-\mathrm{c_{g}}(\tilde{\bm{\Xi}}_{\mu\lambda}^{SL}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{LS})-\mathrm{c_{b}}(\tilde{\bm{\Xi}}_{\mu\lambda}^{SL}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{LS}),$
(153) $\displaystyle\bar{G}^{LL}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=(\bar{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\lambda}])-(\bar{\Omega}_{\mu\lambda}^{LL}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\nu})+(\bar{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\lambda}])$
$\displaystyle-\mathrm{c_{g}}(\bar{\bm{\Xi}}_{\mu\lambda}^{LS}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{SL})-\mathrm{c_{b}}(\bar{\bm{\Xi}}_{\mu\lambda}^{LS}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}_{\kappa\nu}^{SL}),$
(154) $\displaystyle\bar{G}^{LS}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=-(\bar{\Omega}^{LL}_{\mu\lambda}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{SS}_{\kappa\nu})$
$\displaystyle+\mathrm{c_{g}}[(\bar{\bm{\Xi}}_{\mu\nu}^{LS}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])-(\bar{\bm{\Xi}}^{LS}_{\mu\lambda}|g_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\nu})+(\bar{\bm{\Xi}}^{LS}_{\mu\nu}|g_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])]$
$\displaystyle+\mathrm{c_{b}}[\bar{\bm{\Xi}}_{\mu\nu}^{LS}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\lambda}])-(\bar{\bm{\Xi}}^{LS}_{\mu\lambda}|\cdot\mathbf{b}_{12}\cdot|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{LS}_{\kappa\nu})+(\bar{\bm{\Xi}}^{LS}_{\mu\nu}|\cdot\mathbf{b}_{12}\cdot|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\bm{\Xi}}^{SL}_{\kappa\lambda}])],$
(155) $\displaystyle\tilde{\Omega}^{LL}_{\mu\lambda}$
$\displaystyle=\tilde{\chi}^{L{\dagger}}_{\mu}\tilde{\chi}^{L}_{\nu},\quad\tilde{\Omega}^{SS}_{\mu\lambda}=\tilde{\chi}^{S{\dagger}}_{\mu}\tilde{\chi}^{S}_{\nu},$
(156) $\displaystyle\tilde{\bm{\Xi}}^{LS}_{\mu\lambda}$
$\displaystyle=\tilde{\chi}^{L{\dagger}}_{\mu}\bm{\sigma}\tilde{\chi}^{S}_{\nu},\quad\tilde{\bm{\Xi}}^{SL}_{\mu\lambda}=\tilde{\chi}^{S{\dagger}}_{\mu}\bm{\sigma}\tilde{\chi}^{L}_{\nu},$
(157) $\displaystyle\bar{\Omega}^{LL}_{\mu\lambda}$
$\displaystyle=\bar{\chi}^{L{\dagger}}_{\mu}\tilde{\chi}^{L}_{\nu},\quad\bar{\chi}^{L}_{\mu}=\sum_{\nu=1}^{n}\chi^{L}_{\nu}\bar{U}^{LL}_{\nu\mu},\quad\bar{\mathbf{U}}^{LL}=(\mathbf{M}^{LL})^{-1}\tilde{\mathbf{M}}_{+}\mathbf{R}_{+},$
(158) $\displaystyle\bar{\bm{\Xi}}^{LS}_{\mu\lambda}$
$\displaystyle=\bar{\chi}^{L{\dagger}}_{\mu}\bm{\sigma}\tilde{\chi}^{S}_{\nu}.$
(159)

Note that $\tilde{\mathbf{G}}^{XY}[\tilde{\mathbf{D}}]$ in Eqs. (150)-(153)
are completely parallel to $\mathbf{G}^{XY}[\mathbf{D}]$ in Eqs. (51)-(57), in
the sense that they can be obtained simply by replacing $\Omega^{VV}$,
$\bm{\Xi}^{VW}$, and $\mathbf{D}^{VW}$ ($V,W=L,S$) with $\tilde{\Omega}^{VV}$,
$\tilde{\bm{\Xi}}^{VW}$, and $\tilde{\mathbf{D}}$, respectively, in the
latter. However, as far as implementation is concerned, it is more
advantageous to first construct
$\mathbf{L}_{+}^{\mathrm{NESC/SESC,2e}}[\mathbf{D}]$ by transforming
$\tilde{\mathbf{D}}$ back to $\mathbf{D}^{XY}$ (cf. Eq. (146)), and then do
the matrix transformation (140).

So far no approximation has been made, for what has been done is merely to
convert a single $4n$-by-$4n$ matrix DHF equation (43) for $2n$ PESs and $2n$
NESs to two mutually coupled $2n$-by-$2n$ equations, Eqs. (135) and
(124)/(127), for the $2n$ PESs alone. Once these equations are solved via a
dual-level iteration scheme[50, 99], the PESs of the parent matrix DHF
equation (43) can be reproduced up to machine accuracy, thereby justifying the
name “exact two-component”[98]. However, this is neither computationally
favorable (even for one-electron systems) nor necessary in practice. What
really matters is to find an accurate and easily accessible approximation to
the decoupling matrix $\mathbf{X}$, so as to solve the decoupling condition
(124)/(127) indirectly. Taking $\mathbf{X}$ as a matrix functional
$\mathbf{X}[V_{\mathrm{eff}}]$ of the effective potential $V_{\mathrm{eff}}$
(which itself is a functional of the four-component molecular density matrix
$\mathbf{D}$ (36)), and the interatomic interaction strength as a formal
expansion parameter, various approximations to $\mathbf{X}$ can readily be
envisaged[45]. The first approximation is obviously the AtomX (82)[59], where
each atomic $\mathbf{X}_{A}$ is derived from the eigenvectors of the matrix
DHF equation (43) for a spherically averaged and unpolarized atomic
configuration. Note that the renomalization $\mathbf{R}_{+}$ (129) is still of
full dimension. This approximation for $\mathbf{L}^{\mathrm{NESC}}_{+}$ in
$\mathbf{F}_{+}^{\mathrm{X2C}}$ stays obviously in the same spirit as
P4C/Q4C[12], and is hence unsurprisingly very accurate, not only for ground
state energies of molecular systems[44, 45], but also for electric[110, 111,
112, 113, 114, 115] and magnetic[116, 117] response properties, analytic
energy gradient and Hessian[118], as well as periodic systems[119]. There have
been attempts[120, 121, 122] to approximate the renormalization matrix
$\mathbf{R}_{+}$ (129) also as the superposition of the atomic ones (in
conjunction with a ‘diagonal local X’ obtained by solving the one-electron
Dirac equation that is block-diagonal in atoms but including all nuclear
attractions), so as to make $\mathbf{U}$ (121) block diagonal in atoms. Since
$\mathbf{R}_{+}$ is much less local than $\mathbf{X}$, such approximation does
introduce discernible errors[123], which are only tolerable for large systems
in view of the dramatic gain in computational efficiency (especially in
gradient and Hessian calculations[118]). An obvious improvement of the AtomX
(82) is a diatomic approximation (DAX), i.e.,

$\displaystyle\mathbf{X}\approx\sum_{F}^{\oplus}\mathbf{X}_{F},$ (160)

where $\mathbf{X}_{F}$ is derived from the diagonalization of
$\mathbf{F}^{\mathrm{Frag}}$ in Eq. (81) for every pair $F$ of atoms. The
common-atom blocks of fragmental matrices can simply be averaged in this case.
This option is not only necessary for $\mathbf{L}^{\mathrm{NESC}}_{+}$ in
$\mathbf{F}_{+}^{\mathrm{X2C}}$ for situations where two heavy atoms are
located too closely (e.g., at a distance shorter than half of a regular bond),
such that the off-diagonal blocks of the molecular $\mathbf{X}$ become
significant[12], but also a must[45, 61] for $\mathbf{F}_{+}^{\mathrm{X2C}}$
in conjunction with $\mathbf{L}^{\mathrm{SESC}}_{+}$ (which is less
accurate[38] than $\mathbf{L}^{\mathrm{NESC}}_{+}$ for an approximate
$\mathbf{X}$). More generally, the DAX can be extended to a fragmental
approximation (FragX)[45], where a fragment can be chosen to include one, two,
or multiple atoms. It is just that the chosen fragments should be kept fixed
for all molecular geometries. The limiting case, where the whole molecule is
treated as a single fragment, has been called model-potential approximation
(ModX)[45]. Here, the following four-component Fock matrix

$\displaystyle\mathbf{F}^{\mathrm{pmf}}=\mathbf{h}+\mathbf{G}[\mathbf{D}^{0}],\quad\mathbf{D}^{0}=\sum_{A}^{\oplus}\mathbf{D}_{A}$
(161)

is to be constructed and diagonalized. The term ‘model potential’ (arising
from the superposition of atomic densities) was first introduced by van
Wüllen[124] to fix the gauge problem of the zeroth-order regular
approximation[39, 40] in the context of DFT. Very recently, this ansatz was
also employed in the extended atomic mean-field (eamf) approach[49]. Since the
ModX was already introduced[45] long before, and there is no difference
between X2C-KS[45] and X2C-HF[49] in the context of constructing the unitary
transformation (121), a new name is hardly justified. Nonetheless, the term
‘pre-molecular mean-field’ (pmf) seems to be a better characterization of Eq.
(161) than both ‘model potential’ and ‘extended atomic mean-field’, for a pre-
molecule is by definition the superposition of atoms. Moreover, the
abbreviation pmf is in closer analogy with mmf (molecular mean-field)[46],
where it is the full, converged molecular four-component Fock matrix that is
converted to $\mathbf{F}_{+}^{\mathrm{X2C}}$ (136). Therefore, the ModX will
be renamed to pmfX from now on. As an approximation of pmfX, the 1eX[43]
obtained by diagonalizing only the first, one-electron term of Eq. (161) is
widely used. However, it is not accurate enough for magnetic properties[116,
117], and cannot be applied to periodic systems. As an alternative
approximation of pmfX, Eq. (161) was replaced by

$\displaystyle\mathbf{F}^{\mathrm{amf}}=\mathbf{h}+\sum_{A}^{\oplus}\mathbf{G}_{A}[\mathbf{D}_{A}]$
(162)

in the so-called atomic mean-field (amf) approach[49] for the correction of
2ePCEs. However, $\mathbf{F}^{\mathrm{amf}}$ (162) is not a legitimate
Hamiltonian, for there is no physical justification for ignoring the two-
electron but retaining the one-electron interatomic couplings. In particular,
the (scalar) electrostatic interaction between electrons is known to be long-
ranged. Numerical experimentations[125] do reveal that
$\mathbf{F}^{\mathrm{amf}}$ may have spurious solutions, especially when
uncontracted RKB basis sets are used. As such, this amfX should be dumped
definitely. In contrast, the Hamiltonians $\sum_{A}^{\oplus}\mathbf{F}_{A}$,
$\mathbf{F}^{\mathrm{Frag}}$ (81), $\mathbf{F}^{\mathrm{pmf}}$ (161), and
$\mathbf{h}$ (33) employed for deriving the AtomX, FragX, pmfX, and 1eX,
respectively, are all well behaved.

Having discussed the various ways for approximating the decoupling matrix
$\mathbf{X}$, possible approximations to the elements of
$\mathbf{F}_{+}^{\mathrm{X2C}}$ (136) should further be pursued. The very
first of these is the MDM approximation (an equivalent of MSCA[45]) underlying
Q4C (see Sec. III). That is, $\mathbf{G}^{\mathrm{NESC}}[\mathbf{D}]$ in
$\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\mathbf{D}]$ (87) is replaced with
$\mathbf{G}^{\mathrm{NESC}}[\mathbf{D}^{0}]$, so as to rewrite the NESC-based
$\mathbf{F}_{+}^{\mathrm{X2C}}$ (136) as

$\displaystyle\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$
$\displaystyle=\tilde{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{NESC}}+\tilde{\mathbf{G}}^{\mathrm{NR}}[\tilde{\mathbf{D}}],$
(163) $\displaystyle\tilde{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{NESC}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,1e}}\mathbf{R}_{+}+\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}},\quad\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}=\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NESC}}[\mathbf{D}^{0}]\mathbf{R}_{+},$
(164) $\displaystyle\tilde{G}^{\mathrm{NR}}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=(\tilde{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\lambda}])-(\tilde{\Omega}_{\mu\lambda}^{LL}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\nu})$
(165)
$\displaystyle=(\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]\mathbf{R}_{+})_{\mu\nu}.$
(166)

As a matter of fact, $\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (163) can be
obtained directly by the $\mathbf{U}$-transformation (121) of pmfDHF (116). It
has the same structure as $\mathbf{F}_{+}^{\mathrm{Q4C}}$ (97). However, there
exists a subtle but important difference in between: unlike the
renormalization $\mathbf{R}_{+}$ (129) in $\tilde{\mathbf{G}}^{\mathrm{NR}}$
(166), the $\mathbf{a}$ matrix in $\bar{\mathbf{G}}^{\mathrm{NR}}$ (99)/(103)
is merely collection of the contraction coefficients of the A4Ss (cf. Eq.
(71)). That is, $\bar{\mathbf{G}}^{\mathrm{NR}}$ (99)/(103) is a true analog
of the nonrelativistic Coulomb interaction, whereas
$\tilde{\mathbf{G}}^{\mathrm{NR}}$ (166) represents a renormalization of such
interaction, which has significant impact on the energy levels of heavy
elements. By virtue of the identity (cf. Eq (87)),

$\displaystyle\mathbf{G}^{\mathrm{NESC}}[\mathbf{D}^{0}]=\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\mathbf{D}^{0}]-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\mathbf{D}^{LL}_{A}],$
(167)

$\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$ in Eq. (164) can be
written as

$\displaystyle\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\mathbf{D}^{0}]\mathbf{R}_{+}-\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\mathbf{D}^{LL}_{A}]\mathbf{R}_{+}$
(168)
$\displaystyle=\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}-\Delta\mathbf{G}^{\mathrm{NESC}},$
(169)

where

$\displaystyle\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\mathbf{D}^{0}]\mathbf{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\tilde{\mathbf{D}}_{A}]),$
(170) $\displaystyle\Delta\mathbf{G}^{\mathrm{NESC}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\mathbf{D}^{LL}_{A}]\mathbf{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\tilde{\mathbf{D}}_{A}]$
(171)
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]\mathbf{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}]$
$\displaystyle-\mathbf{R}_{+}^{\dagger}(\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\mathbf{D}^{LL}_{A}])\mathbf{R}_{+}+(\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}]-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\tilde{\mathbf{D}}_{A}])$
(172)
$\displaystyle=\tilde{\mathbf{G}}^{\mathrm{NR}}[\tilde{\mathbf{D}}]-\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}]-\Delta\Delta\mathbf{G}^{\mathrm{NR}},$
(173) $\displaystyle\Delta\Delta\mathbf{G}^{\mathrm{NR}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{NR}}[\Delta\mathbf{D}^{LL}]\mathbf{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\Delta\tilde{\mathbf{D}}],$
(174) $\displaystyle\Delta{\mathbf{D}}^{LL}$
$\displaystyle=\mathbf{D}^{LL}-\sum_{A}^{\oplus}\mathbf{D}_{A}^{LL},\quad\Delta\tilde{\mathbf{D}}=\tilde{\mathbf{D}}-\sum_{A}^{\oplus}\tilde{\mathbf{D}}_{A}.$
(175)

Since both $\Delta{\mathbf{D}}^{LL}$ and $\Delta\tilde{\mathbf{D}}$
characterize the deformation and polarization when going from free atoms to
the molecule (due to bonding interactions between valence shells), the two
terms of Eq. (174) essentially cancel each other, thereby leading to

$\displaystyle\Delta\mathbf{G}^{\mathrm{NESC}}\approx\tilde{\mathbf{G}}^{\mathrm{NR}}[\tilde{\mathbf{D}}]-\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}].$
(176)

Eq. (163) can hence be reduced to

$\displaystyle\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$
$\displaystyle=\bar{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{NESC}}+\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}],$
(177) $\displaystyle\bar{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{NESC}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,1e}}\mathbf{R}_{+}+\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}},$
(178)

which amounts to shifting $\Delta\mathbf{G}^{\mathrm{NESC}}$ (176) from
$\tilde{\mathbf{G}}^{\mathrm{NR}}$ (165) to
$\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$ (164). That is,
$\Delta\mathbf{G}^{\mathrm{NESC}}$ (176) behaves as an additional correction
of the 2e-PCEs. Eq. (177), in conjunction with the pmfX derived from Eq.
(161), has been dubbed eamfX2C[49] (or preferably pmfX2C), where the
expressions (177) and (178) arise naturally by regarding Eq. (170) as the pmf
approximation to the true term
$\mathrm{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{NESC,2e}}[\mathbf{D}]\mathrm{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}]$.
Note that any approximation to the decoupling matrix $\mathbf{X}$ can be used
in $\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (177), although the pmfX appears to
be a natural choice, for $\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$
(170) can be obtained for free in this case.

Similarly, if the SESC-based $\mathbf{F}_{+}^{\mathrm{X2C}}$ (136) is to be
adopted, we will have under the MDM approximation

$\displaystyle\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$
$\displaystyle=\tilde{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{SESC}}+\tilde{\mathbf{g}}^{\mathrm{NR}}[\tilde{\mathbf{D}}],$
(179) $\displaystyle\tilde{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{SESC}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{SESC,1e}}\mathbf{R}_{+}+\tilde{\mathbf{G}}^{\mathrm{SESC}}_{\mathrm{eff}},\quad\tilde{\mathbf{G}}^{\mathrm{SESC}}_{\mathrm{eff}}=\mathbf{R}_{+}^{\dagger}\mathbf{G}^{\mathrm{SESC}}[\mathbf{D}^{0}]\mathbf{R}_{+},$
(180) $\displaystyle\tilde{g}^{\mathrm{NR}}_{\mu\nu}[\tilde{\mathbf{D}}]$
$\displaystyle=\frac{1}{2}\\{[(\bar{\Omega}_{\mu\nu}^{LL}|g_{12}|\mathrm{Tr}_{2}[\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\lambda}])-(\bar{\Omega}_{\mu\lambda}^{LL}|g_{12}|\tilde{\mathbf{D}}_{\lambda\kappa}\tilde{\Omega}^{LL}_{\kappa\nu})]+\mathrm{c.c.}\\}$
(181)
$\displaystyle=\frac{1}{2}\\{(\mathbf{R}_{+}^{\dagger}\tilde{\mathrm{M}}_{+}(\mathbf{M}^{LL})^{-1}\mathbf{G}^{\mathrm{NR}}[\mathbf{D}^{LL}]\mathbf{R}_{+})_{\mu\nu}+\mathrm{c.c.}\\},$
(182)

whereas $\mathbf{G}^{\mathrm{SESC}}_{\mathrm{eff}}[\mathbf{D}^{0}]$ collects
all small-component-containing terms of
$\mathbf{L}_{+}^{\mathrm{SESC,2e}}[\mathbf{D}^{0}]$ (141). Following the same
procedure going from Eq. (167) to Eq. (176), Eq. (179) can be reduced to

$\displaystyle\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$
$\displaystyle=\bar{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{SESC}}+\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}],$
(183) $\displaystyle\bar{\mathbf{h}}_{\mathrm{eff}}^{\mathrm{SESC}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{SESC,1e}}\mathbf{R}_{+}+\bar{\mathbf{G}}^{\mathrm{SESC}}_{\mathrm{eff}},$
(184) $\displaystyle\bar{\mathbf{G}}^{\mathrm{SESC}}_{\mathrm{eff}}$
$\displaystyle=\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{SESC,2e}}[\mathbf{D}^{0}]\mathbf{R}_{+}-\mathbf{G}^{\mathrm{NR}}[\sum_{A}^{\oplus}\tilde{\mathbf{D}}_{A}]),$
(185)

which is the SESC variant of pmfX2C.

It has been attempted to replace the eamf/pmf-2ePCE correction
$\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$ (170) in Eq. (178) with the
superposition of the atomic ones, viz.,

$\displaystyle\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}=\sum_{A}^{\oplus}\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff},A},\quad\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff},A}=\mathrm{R}_{+,A}^{\dagger}\mathbf{L}_{+,A}^{\mathrm{NESC,2e}}[\mathbf{D}_{A}]\mathrm{R}_{+,A}-\mathbf{G}^{\mathrm{NR}}_{A}[\tilde{\mathbf{D}}_{A}],$
(186)

so as to reduce eamfX2C/pmfX2C to amfX2C[49]. However, even when combined with
a proper approximation (e.g., 1eX) to the decoupling matrix $\mathbf{X}$ (in
lieu of the original amfX[49] derived from the ill-behaved Hamiltonian (162)),
this amfX2C has sizeable errors (up to 5 $m\mathrm{E}_{h}$) for the energy
levels of frontier spinors of polar systems (e.g., HI) with extended basis
sets[125], a point that was not observed in Ref. 49 due to the use of low-
quality basis sets. The reason should be ascribed to the long-range nature of
electrostatic electron-electron interactions. In contrast, when
$\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff},A}$ is composed only of one-
center two-electron spin-orbit (2eSO) integrals (contracted with the four-
component atomic density matrix $\mathbf{D}_{A}$), the resulting SOX2CAMF
approach[47] turns out to be very successful. It appears that two-electron
scalar relativistic (2eSC) PCEs are either ignored completely or accounted for
by the pmf scheme. For the same reason, the scalar Gaunt/Breit integrals
should also be removed from SOX2CAMF that incorporates the full atomic mean-
field Gaunt/Breit integrals[48].

It has been a common practice to take the spin-free part of
$\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (177), in conjunction with
$\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}=0$, as the Hamiltonian in
scalar relativistic mean-field calculations. Since the spin-free part of
$\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$ (170) can readily be
obtained, it is strongly recommended to correct the 2eSC-PCEs even in scalar
relativistic calculations, so as to achieve a balanced description of core and
valence electrons.

At this stage, some general remarks should be made.

  1. (1)

Although precisely the same in structure, Q4C[45, 60, 61] is much simpler than
eamfX2C/pmfX2C[49] (defined by Eqs. (161), (170), and (177)), especially for
response properties (due to the presence of pmfX and renormalization
$\mathbf{R}_{+}$ in eamfX2C; see Ref. 118). Moreover, unlike pmfX2C, Q4C does
not suffer from 2cPCEs even in the fluctuation potential governing electron
correlation. It is hence unclear why eamX2C, instead of Q4C, was characterized
as a ‘fundamental milestone toward a universal and reliable relativistic two-
component quantum-chemical approach’[49].

  2. (2)

Since $\tilde{\mathbf{G}}^{\mathrm{NR}}[\tilde{\mathbf{D}}]$ has little
computational overhead compared to
$\mathbf{G}^{\mathrm{NR}}[\tilde{\mathbf{D}}]$,
$\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (163)/(179), in conjunction with any
acceptable approximation to the decoupling matrix $\mathbf{X}$, is also a
valuable variant of X2C. Conceptually, it is even preferred over
$\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (177)/(183), for it is the only
formulation that is consistent with the renormalization procedure when going
from pmfDHF (116) to pmfX2C. In particular, the difference between
$\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$ and
$\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$, i.e.,
$\Delta\Delta\mathbf{G}^{\mathrm{NR}}$ (174), has discernible effects on the
innermost shells of heavy elements and is hence relevant for properties (e.g.,
contact densities and shifts[49]) that are very sensitive to such electronic
shells.

  3. (3)

SESC is much simpler than NESC. As can be seen from Eq. (149), SESC does not
require the very expensive two-electron term $\tilde{\mathbf{G}}^{SS}$ (153)
that enters Eq. (148) for NESC. Even the one-electron term
$\mathbf{L}_{+}^{\mathrm{SESC,1e}}$ (139) of SESC is much simpler than that
[Eq. (85)] of NESC. In particular, for a pure density functional,
$\mathbf{R}_{+}^{\dagger}\mathbf{L}_{+}^{\mathrm{SESC,1e}}\mathbf{R}_{+}$
would be the full SESC-KS Hamiltonian, which does not have an explicit spin-
orbit operator (since
$\mathbf{L}_{+}^{\mathrm{UESC,1e}}=\mathbf{T}+\mathbf{V}_{\mathrm{KS}}[\rho]\mathbf{X}$
in $\mathbf{L}_{+}^{\mathrm{SESC,1e}}$), but does not miss any spin-orbit
couplings described by the given functional[45]. However, when the pmfX is
employed, the simplicity of SESC in the form of Eq. (179)/(183) over NESC
(163)/(177) is lost, for the expensive
$\mathbf{G}^{\mathrm{NESC}}[\sum_{A}^{\oplus}\mathbf{D}_{A}]$ required for
assembling
$\tilde{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$/$\bar{\mathbf{G}}^{\mathrm{NESC}}_{\mathrm{eff}}$
is already available in this case. Nonetheless, SESC is still much simpler
than NESC if the MDM approximation is not employed.

  4. (4)

The 1CSC approximation of the ERIs [cf. Eqs. (63) to (65)] can also be applied
to $\mathbf{F}^{\mathrm{Frag}}$ (81), $\mathbf{F}_{+}^{\mathrm{Q4C}}$
(97)/(102), $\mathbf{F}_{+}^{\mathrm{X2C}}$ (136), $\mathbf{F}^{\mathrm{pmf}}$
(161), $\tilde{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (163)/(179), and
$\bar{\mathbf{F}}_{+}^{\mathrm{X2C}}$ (177)/(183), so as to render them share
exactly the same ERIs as the parent four-component Fock matrix $\mathbf{F}$
(43).

  5. (5)

Although only defined algebraically, it is still possible to separate the
various forms of the X2C Hamiltonian into spin-free and spin-dependent terms.
As a matter of fact, taking the spin-free part of X2C (sf-X2C) as the zeroth
order, a family of spin-dependent operators of finite orders in spin-orbit
(SO) interaction can be obtained by means of matrix perturbation theory[107,
109] or derivative technique[126, 127] (see Ref. 128 for a comprehensive
review). Among these, the Douglas-Kroll-Hess type of SO operators (so-
DKHn)[107, 109] are particularly simple and variationally stable (so-DKH1 is
even computationally the same as the Breit-Pauli spin-orbit operator).

  6. (6)

As a conceptual point, it should be pointed that, while the decoupling matrix
$\mathbf{X}$ can in principle be obtained within the two-component
framework[50, 99], it is the borrowing of easily accessible four-component
information that renders X2C really effective. The same applies also to P4C
and Q4C.

  7. (7)

If wanted, the unoccupied NESs can be obtained by a single diagonalization of
Eq. (187) upon convergence of the PESs,

$\displaystyle\mathbf{F}_{-}^{\mathrm{X2C}}\tilde{\mathbf{C}}_{-}$
$\displaystyle=\mathbf{M}^{SS}\tilde{\mathbf{C}}_{-}\mathbf{E}_{-},$ (187)
$\displaystyle\mathbf{F}_{-}^{\mathrm{X2C}}$
$\displaystyle=(\mathbf{U}^{\dagger}\mathbf{F}\mathbf{U})^{SS}=\sum_{X,Y=L,S}(\mathbf{U}^{\dagger})^{SX}\mathbf{F}^{XY}\mathbf{U}^{YS}$
$\displaystyle=\mathbf{R}_{-}^{\dagger}\mathbf{L}^{\mathrm{X}}_{-}\mathbf{R}_{-},\quad\mathrm{X}=\mathrm{NESC,SESC},$
(188) $\displaystyle\mathbf{L}^{\mathrm{NESC}}_{-}$
$\displaystyle=\mathbf{F}^{SS}+\mathbf{F}^{SL}\tilde{\mathbf{X}}+\tilde{\mathbf{X}}^{\dagger}\mathbf{F}^{LS}+\tilde{\mathbf{X}}^{\dagger}\mathbf{F}^{LL}\tilde{\mathbf{X}},$
(189) $\displaystyle\mathbf{L}^{\mathrm{SESC}}_{-}$
$\displaystyle=\frac{1}{2}[\tilde{\mathbf{M}}_{-}(\mathbf{M}^{SS})^{-1}\mathbf{L}^{\mathrm{UESC}}_{-}+\mathrm{c.c.}],$
(190) $\displaystyle\mathbf{L}^{\mathrm{UESC}}_{-}$
$\displaystyle=\mathbf{F}^{SS}+\mathbf{F}^{SL}\tilde{\mathbf{X}},$ (191)
$\displaystyle\mathbf{F}^{XY}$
$\displaystyle=\mathbf{h}^{XY}+\mathbf{G}^{XY}[\mathbf{D}].$ (192)

This is a cheap step, for the relevant two-electron ERIs
$G^{XY}_{\mu\nu}[\mathbf{D}]$ are already available. Here, the four-component
density matrix $\mathbf{D}$ is back-transformed from the converged two-
component density matrix $\tilde{\mathbf{D}}$ (cf. Eq. (146)). Note that the
same equation (187) can also be employed to obtain the NESs of Q4C. It is just
that the four-component density matrix is $\bar{\mathbf{D}}$, back-transformed
from the Q4C density matrix $\bar{\mathbf{P}}$ (cf. Eq. (79)). Alternatively,
the four-component DHF equation (43) can be assembled upon convergence of the
X2C/Q4C calculation and is then diagonalized to obtain the PESs and NESs
simultaneously.

  8. (8)

The quaternion form of X2C can be obtained by the
$\mbox{}^{\mathrm{Q}}\mathbf{U}$ transformation (295) of
$\mathbf{F}^{\mathrm{X2C}}_{+}$/$\bar{\mathbf{F}}^{\mathrm{X2C}}_{+}$. The X4C
variant of KROHF can then readily be formulated (see Appendix A).

  9. (9)

Since the leading quantum electrodynamics effect (i.e., vacuum polarization
and electron self-energy) can be described by an effective, HF-like one-body
potential [62, 8] that can be added to the one-body Dirac operator (19), it is
readily accounted for in DHF, Q4C or X2C calculations. The situation is
particularly simple when the eQED potential is fitted into a model spectral
form[63, 64]. For a deep understanding of the fundamentals of QED, see
comprehensive reviews[11, 129, 12, 130].

##  V Many-electron Relativistic Hamiltonians

Having determined the PESs (and NESs) by DHF, Q4C or X2C, a generic normal-
ordered many-electron Hamiltonian can be written down

$\displaystyle H_{n}=H-\langle
0|H|0\rangle=F^{\mathrm{X}}_{pq}\\{a^{p}_{q}\\}_{n}+\frac{1}{2}g_{pq}^{rs}\\{a^{pq}_{rs}\\}_{n},\quad\mathrm{X=DHF,Q4C,X2C,eQED},$
(193)

for post-HF calculations. Here, $|0\rangle$ is a reference state, with respect
to which the normal ordering of the one-body
($a^{p}_{q}=a_{p}^{\dagger}a_{q}$) and two-body
($a^{pq}_{rs}=a_{p}^{\dagger}a_{q}^{\dagger}a_{s}a_{r}$) excitation operators
has been made. Like the 4C case, the two-body operator of Q4C incorporates
automatically the full Coulomb-Gaunt/Breit interaction (15). Since the large
and small components of the spinor basis functions in Q4C share the same
expansion coefficients, the AO ERIs can be transformed to the MO
representation as a whole, instead of component-wise as in the 4C case. When
the AO ERIs refer to those of the renormalized basis functions (cf. Eq.
(145)), the two-body operator of X2C will also incorporate the full Coulomb-
Gaunt/Breit interaction, which is equivalent to working with the M4Ss
recovered from the X2C spinors (cf. Eq. (143)). As such, the index
transformation of the two-body operator of X2C is as expensive as that of 4C
and more expensive than that of Q4C. The computation of X2C is simplified
greatly by using only the untransformed two-body operator (i.e.,
$\\{g_{pq}^{rs}\\}$ are just the ERIs of the Coulomb interaction over the X2C
spinors). However, X2C will then suffer from 2ePCEs in electron-electron
interaction, which are particularly prominent for deep-core properties of
heavy elements[49, 131]. Not only so, the lack of genuine 2e-SO interaction
will fail to describe SO splittings between states differing by two electrons.
The simplest remedy[132, 133] of this failure is to include the one-center
2e-SO AO ERIs (resulting from spin separation of the transformed X2C two-body
operator[107, 109]) and transform them only to the active orbitals spanning
the target states. Beyond these no-pair relativistic Hamiltonians is the
effective QED Hamiltonian[62, 8], where the normal ordering is taken with
respective to the filled Dirac sea but meanwhile incorporating charge
conjugation symmetry[12]. This is the only correct and complete QED
Hamiltonian[130] in the same form of Eq. (193). It is just that the one-body
operator further includes the effective potential describing vacuum
polarization and electron self-energy and that the orbital indices refer to
both PESs and NESs. Unlike no-pair correlation, the full QED correlation
energy is independent of the ways of generating the orbitals (like
nonrelativistic full configuration interaction)[8], for the _filled_ NESs are
also correlated therein. Note that the frequency-dependent Breit interaction
must be employed for this purpose. For possible means of treating relativity,
correlation, and QED effects simultaneously, see Ref. 15.

The above spinor-based Hamiltonians are imperative for core properties of
heavy atoms or valence properties involving $np$ ($n\geq 5$) orbitals.
However, for chemical systems with moderate spin-orbit couplings, it is more
appealing to invoke a two-step Hamiltonian that treats scalar relativity and
spin-orbit coupling separately. The sf-X2C+so-DKHn variant[107, 109] is
arguably the best one for this purpose. The great advantage here lies in that
real-valued orbitals can be used, so as to facilitate the treatment of
electron correlation, on top of which spin-orbit couplings can further be
added in one way or another. Methodologies and applications along this line
have recently been summarized[134, 135, 136] and are hence not repeated here.

##  VI Conclusions and outlook

It has been shown that the DHF, Q4C, and X2C equations can be recast into the
same mean-field form by making use of the MDM approximation for the small-
component charge/current density functions. They also share exactly the same
relativistic integrals that can be simplified by using the 1CSC approximation.
As such, it is a matter of taste which variant is to be adopted for real-life
applications. Nevertheless, Q4C is manifestly most efficient due to the use of
most compact bases and all-together integral transformation for subsequent
treatment of electron correlation. In particular, it has no 2ePCEs even in
electron-election interaction and is simpler than X2C for response properties.
QED effects, including not only the one-body vacuum polarization and electron
self-energy but also the two-body correlation of NESs, can readily be
accounted for in precision spectroscopic calculations.

###### Acknowledgements.

This work was supported by the National Natural Science Foundation of China
(Grant Nos. 22373057, 21833001, and 21973054) and Mount Tai Scholar Climbing
Project of Shandong Province.

## Data Availability Statement

The data that support the findings of this study are available within the
article.

## Conflicts of interest

There are no conflicts to declare.

##  Appendix A Quaternion DHF, KUHF, and KROHF

Although the first (third) and second (fourth) columns of the RKB basis
$|\xi_{\mu}\rangle$ (27) are manifestly time-reversal related (cf. Eq. (261)),
it turns out that the overall time-reversal structure of the 4-by-4 matrix
(27) can be made more transparent only by the following unitary
transformation[85],

$\displaystyle\mathbf{Q}$ $\displaystyle=\begin{pmatrix}1&0&0&0\\\ 0&0&1&0\\\
0&1&0&0\\\ 0&0&0&1\end{pmatrix}=\mathbf{Q}^{-1}=\mathbf{Q}^{\dagger},$ (194)
$\displaystyle\mbox{}^{\mathrm{Q}}\xi_{\mu}$
$\displaystyle=\mathbf{Q}\xi_{\mu}\mathbf{Q}^{\dagger}=\begin{pmatrix}\xi^{\alpha}_{\mu}&\bar{\xi}^{\beta}_{\mu}\\\
\bar{\xi}^{\alpha}_{\mu}&\xi^{\beta}_{\mu}\end{pmatrix}=\begin{pmatrix}\xi^{\alpha}_{\mu}&\bar{\xi}^{\beta}_{\mu}\\\
-\bar{\xi}^{\beta*}_{\mu}&\xi^{\alpha*}_{\mu}\end{pmatrix},\quad\mu\in[1,n],$
(195) $\displaystyle\xi^{\alpha}_{\mu}$
$\displaystyle=\begin{pmatrix}g_{\mu}&0\\\
0&-\frac{\mathbbm{i}}{2c}g^{z}_{\mu}\end{pmatrix},\quad\bar{\xi}^{\beta}_{\mu}=\begin{pmatrix}0&0\\\
0&-\frac{\mathbbm{i}}{2c}(g^{x}_{\mu}-\mathbbm{i}g^{y}_{\mu})\end{pmatrix}.$
(196)

$\mbox{}^{\mathrm{Q}}\xi_{\mu}$ (195) is clearly time-reversal symmetric (cf.
Eq. (262)) and can hence be block-diagonalized by the quaternion unitary
transformation $\mbox{}^{\mathrm{Q}}\mathbf{U}$ (295), viz.,

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}\xi_{\mu}\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{q}}\xi_{\mu}&0\\\
0&\mbox{}^{\mathrm{q}}\xi_{\mu}\end{pmatrix},$ (197)
$\displaystyle\mbox{}^{\mathrm{q}}\xi_{\mu}$
$\displaystyle=\xi^{\alpha}_{\mu}+\bar{\xi}^{\beta}_{\mu}\check{j}=\mbox{}^{0}X_{\mu}+\mathbf{X}_{\mu},\quad\mathbf{X}_{\mu}=\sum_{i=1}^{3}\mbox{}^{i}X_{\mu}e_{i},$
(198) $\displaystyle\mbox{}^{0}X_{\mu}$
$\displaystyle=\begin{pmatrix}g_{\mu}&0\\\
0&0\end{pmatrix},\quad\mbox{}^{i}X_{\mu}=-\frac{1}{2c}\begin{pmatrix}0&0\\\
0&g^{i}_{\mu}\end{pmatrix},$ (199)

where the quaternion units $(e_{1},e_{2},e_{3})$ refer to the ordering
$(z,y,x)$ instead of the usual ordering $(x,y,z)$ ($=-(z,y,x)$), see Eq.
(271). It is said that the 4-by-4 matrix $\mbox{}^{\mathrm{Q}}\xi_{\mu}$ (195)
is algebraically isomorphic to the 2-by-2 (real) quaternion matrix
$\mbox{}^{\mathrm{q}}\xi_{\mu}$ (198) [cf. Eq. (270)].

The $\mathbf{Q}$-transformation (194) of an M4S $\psi_{p}$ amounts to swapping
its second and third components,

$\displaystyle\mbox{}^{\mathrm{Q}}\psi_{p}$
$\displaystyle=\mathbf{Q}\psi_{p}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\psi_{p}^{\alpha}\\\
\mbox{}^{\mathrm{Q}}\psi_{p}^{\beta}\end{pmatrix}=\mbox{}^{\mathrm{Q}}\xi_{\mu}\mbox{}^{\mathrm{Q}}C_{\mu
p},\quad\mbox{}^{\mathrm{Q}}\psi_{p}^{\sigma}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\psi_{p}^{\sigma
L}\\\ \mbox{}^{\mathrm{Q}}\psi_{p}^{\sigma S}\end{pmatrix},$ (200)
$\displaystyle\mbox{}^{\mathrm{Q}}C_{\mu p}$ $\displaystyle=\mathbf{Q}C_{\mu
p}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}C^{\alpha}_{\mu p}\\\
\mbox{}^{\mathrm{Q}}C^{\beta}_{\mu
p}\end{pmatrix},\quad\mbox{}^{\mathrm{Q}}C_{\mu
p}^{\sigma}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}C^{\sigma L}_{\mu p}\\\
\mbox{}^{\mathrm{Q}}C^{\sigma S}_{\mu p}\end{pmatrix},$ (201)

such that the resulting 4-spinor $\mbox{}^{\mathrm{Q}}\psi_{p}$ is grouped
according to spin labels $(\alpha,\beta)$ instead of large and small
components $(L,S)$. The Dirac operator (4) is transformed accordingly,

$\displaystyle\mbox{}^{\mathrm{Q}}h$
$\displaystyle=\mathbf{Q}h\mathbf{Q}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}h^{\beta\alpha}&\mbox{}^{\mathrm{Q}}h^{\beta\beta}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta}\\\
-(\mbox{}^{\mathrm{Q}}h^{\alpha\beta})^{*}&(\mbox{}^{\mathrm{Q}}h^{\alpha\alpha})^{*}\end{pmatrix},$
(202) $\displaystyle\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\alpha,LL}&\mbox{}^{\mathrm{Q}}h^{\alpha\alpha,LS}\\\
\mbox{}^{\mathrm{Q}}h^{\alpha\alpha,SL}&\mbox{}^{\mathrm{Q}}h^{\alpha\alpha,SS}\end{pmatrix}=\begin{pmatrix}V&cp_{z}\\\
cp_{z}&V-2c^{2}\end{pmatrix}=(\mbox{}^{\mathrm{Q}}h^{\beta\beta})^{*},$ (203)
$\displaystyle\mbox{}^{\mathrm{Q}}h^{\alpha\beta}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\beta,LL}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta,LS}\\\
\mbox{}^{\mathrm{Q}}h^{\alpha\beta,SL}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta,SS}\end{pmatrix}=\begin{pmatrix}0&c(p_{x}-\mathbbm{i}p_{y})\\\
c(p_{x}-\mathbbm{i}p_{y})&0\end{pmatrix}=-(\mbox{}^{\mathrm{Q}}h^{\beta\alpha})^{*},$
(204)

which can also be block-diagonalized by the quaternion unitary transformation
(295),

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}h\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{q}}h&0\\\
0&\mbox{}^{\mathrm{q}}h\end{pmatrix},$ (205)

where

$\displaystyle\mbox{}^{\mathrm{q}}h$
$\displaystyle=\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}+\mbox{}^{\mathrm{Q}}h^{\alpha\beta}\check{j}=\sum_{i=0}^{3}\mbox{}^{i}he_{i},$
(206) $\displaystyle\mbox{}^{0}h_{D}$ $\displaystyle=\begin{pmatrix}V&0\\\
0&V-2c^{2}\end{pmatrix},\quad\mbox{}^{i}h=\begin{pmatrix}0&-c\partial_{i}\\\
-c\partial_{i}&0\end{pmatrix}.$ (207)

Precisely in the same way, the one-body matrix elements $h_{\mu\nu}$ can be
transformed to quaternion form,

$\displaystyle\mbox{}^{\mathrm{Q}}h_{\mu\nu}$
$\displaystyle=\mathbf{Q}h_{\mu\nu}\mathbf{Q}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta}_{\mu\nu}\\\
\mbox{}^{\mathrm{Q}}h^{\beta\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}h^{\beta\beta}_{\mu\nu}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}h^{\alpha\beta}_{\mu\nu}\\\
-\mbox{}^{\mathrm{Q}}h^{\alpha\beta*}_{\mu\nu}&\mbox{}^{\mathrm{Q}}h^{\alpha\alpha*}_{\mu\nu}\end{pmatrix},$
(208)
$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}h_{\mu\nu}\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{q}}h_{\mu\nu}&0\\\
0&\mbox{}^{\mathrm{q}}h_{\mu\nu}\end{pmatrix},$ (209)
$\displaystyle\mbox{}^{\mathrm{q}}h_{\mu\nu}$
$\displaystyle=\langle\mbox{}^{\mathrm{q}}\xi_{\mu}|\mbox{}^{\mathrm{q}}h|\mbox{}^{\mathrm{q}}\xi_{\nu}\rangle=\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}_{\mu\nu}+\mbox{}^{\mathrm{Q}}h^{\alpha\beta}_{\mu\nu}\check{j}=\sum_{i=0}^{3}\mbox{}^{i}h_{\mu\nu}e_{i},$
(210) $\displaystyle\mbox{}^{\mathrm{Q}}h^{\alpha\alpha}_{\mu\nu}$
$\displaystyle=\begin{pmatrix}V_{\mu\nu}&T_{\mu\nu}\\\
T_{\mu\nu}&\frac{1}{4c^{2}}\langle g_{\mu}^{j}|V_{N}|g_{\nu}^{j}\rangle-
T_{\mu\nu}\end{pmatrix}+\frac{1}{4c^{2}}\begin{pmatrix}0&0\\\
0&\mathbbm{i}\epsilon_{zjk}\langle
g_{\mu}^{j}|V|g_{\nu}^{k}\rangle\end{pmatrix},\quad j,k\in x,y,z,$ (211)
$\displaystyle\mbox{}^{\mathrm{Q}}h^{\alpha\beta}_{\mu\nu}$
$\displaystyle=\frac{1}{4c^{2}}\begin{pmatrix}0&0\\\
0&(\epsilon_{yjk}+\mathbbm{i}\epsilon_{xjk})\langle
g_{\mu}^{j}|V|g_{\nu}^{k}\rangle\end{pmatrix},\quad j,k\in x,y,z.$ (212)

To obtain the quaternion form
$\mbox{}^{\mathrm{q}}G_{\mu\nu}[\mbox{}^{\mathrm{q}}\mathbf{D}]$ of the two-
electron terms $G_{\mu\nu}[\mathbf{D}]$ (32), we first calculate the
quaternion form $\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}$ of the overlap charge
distribution function $\Omega_{\mu\nu}$ (38):

$\displaystyle\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}=\mathbf{Q}\Omega_{\mu\nu}\mathbf{Q}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\beta\alpha}&\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\beta\beta}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\beta}\\\
-(\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\beta})^{*}&(\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\alpha})^{*}\end{pmatrix},$
(213)
$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}&0\\\
0&\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}\end{pmatrix},$ (214)

where

$\displaystyle\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\alpha}+\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\beta}\check{j}=(\mbox{}^{\mathrm{q}}\xi_{\mu})^{\dagger}\mbox{}^{\mathrm{q}}\xi_{\nu}=\sum_{i=0}^{3}\mbox{}^{i}\Omega_{\mu\nu}e_{i},$
(215) $\displaystyle\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\alpha}$
$\displaystyle=\begin{pmatrix}g_{\mu}g_{\nu}&0\\\
0&\frac{1}{4c^{2}}g^{j}_{\mu}g^{j}_{\nu}\end{pmatrix}+\frac{1}{4c^{2}}\begin{pmatrix}0&0\\\
0&\mathbbm{i}\epsilon_{zjk}g^{j}_{\mu}g^{k}_{\nu}\end{pmatrix},\quad j,k\in
x,y,z,$ (216) $\displaystyle\mbox{}^{\mathrm{Q}}\Omega_{\mu\nu}^{\alpha\beta}$
$\displaystyle=\frac{1}{4c^{2}}\begin{pmatrix}0&0\\\
0&(\epsilon_{yjk}+\mathbbm{i}\epsilon_{xjk})g^{j}_{\mu}g^{k}_{\nu}\end{pmatrix},\quad
j,k\in x,y,z.$ (217)

Likewise, the complex quaternion form $\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}$
(cf. Eq. (280)) of the time-reversal antisymmetric overlap current
distribution vectors $\Xi^{i}_{\mu\nu}$ (40) can be obtained by block-
diagonalizing $\mbox{}^{\mathrm{Q}}\Xi^{i}_{\mu\nu}$,

$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{i}_{\mu\nu}=\mathbf{Q}\Xi^{i}_{\mu\nu}\mathbf{Q}^{\dagger}=h\begin{pmatrix}\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\beta}_{\mu\nu}\\\
\mbox{}^{\mathrm{Q}}\Xi^{i,\beta\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}\Xi^{i,\beta\beta}_{\mu\nu}\end{pmatrix}=h\begin{pmatrix}\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\alpha}_{\mu\nu}&\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\beta}_{\mu\nu}\\\
-(\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\beta}_{\mu\nu})^{*}&(\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\alpha}_{\mu\nu})^{*}\end{pmatrix},$
(218)
$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}\Xi^{i}_{\mu\nu}\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}&0\\\
0&\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}\end{pmatrix},$ (219)

where

$\displaystyle\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}$
$\displaystyle=h(\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\alpha}_{\mu\nu}+\mbox{}^{\mathrm{Q}}\Xi^{i,\alpha\beta}_{\mu\nu}\check{j}),\quad
i\in x,y,z,$ (220)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{x,\alpha\alpha}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&-g_{\mu}g_{\nu}^{x}-\mathbbm{i}g_{\mu}g_{\nu}^{y}\\\
g_{\mu}^{x}g_{\nu}-\mathbbm{i}g_{\mu}^{y}g_{\nu}&0\end{pmatrix},$ (221)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{x,\alpha\beta}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&g_{\mu}g_{\nu}^{z}\\\
g_{\mu}^{z}g_{\nu}&0\end{pmatrix},$ (222)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{y,\alpha\alpha}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&-g_{\mu}g_{\nu}^{y}+\mathbbm{i}g_{\mu}g_{\nu}^{x}\\\
g_{\mu}^{y}g_{\nu}+\mathbbm{i}g_{\mu}^{x}g_{\nu}&0\end{pmatrix},$ (223)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{y,\alpha\beta}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&-\mathbbm{i}g_{\mu}g_{\nu}^{z}\\\
-\mathbbm{i}g_{\mu}^{z}g_{\nu}&0\end{pmatrix},$ (224)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{z,\alpha\alpha}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&-g_{\mu}g_{\nu}^{z}\\\
g_{\mu}^{z}g_{\nu}&0\end{pmatrix},$ (225)
$\displaystyle\mbox{}^{\mathrm{Q}}\Xi^{z,\alpha\beta}_{\mu\nu}$
$\displaystyle=\frac{1}{2c}\begin{pmatrix}0&-g_{\mu}g_{\nu}^{x}+\mathbbm{i}g_{\mu}g_{\nu}^{y}\\\
-g_{\mu}^{x}g_{\nu}+\mathbbm{i}g_{\mu}^{y}g_{\nu}&0\end{pmatrix}.$ (226)

Note that the symbol $h$ in Eqs. (218) and (220) is an alternative notation
for the square root of minus one, to emphasize that it commutes with the
quaternion units $e_{i}$ in the context of biquaternion algebra. For a double
check, we recalculate directly $\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}$ (220),
by noting first that

$\displaystyle Q\widetilde{\alpha}_{i}Q^{\dagger}$
$\displaystyle=\widetilde{\sigma}_{i}\otimes\sigma_{x}\simeq\mbox{}^{\mathrm{q}}\alpha_{i}=-h\sigma_{x}e_{i},\quad
he_{i}=e_{i}h,\quad i\in 1,2,3,$ (227)

where
$(\widetilde{\alpha}_{1},\widetilde{\alpha}_{2},\widetilde{\alpha}_{3})=(\alpha_{z},\alpha_{y},\alpha_{x})$.
We then have (cf. (331))

$\displaystyle\mbox{}^{\mathrm{q}}\xi^{\alpha{\dagger}}_{\mu}\mbox{}^{\mathrm{q}}\alpha_{i}$
$\displaystyle=-h(\mbox{}^{0}X_{\mu}-\mathbf{X}_{\mu})\sigma_{x}e_{i}$ (228)
$\displaystyle=-h(\mbox{}^{i}X_{\mu}\sigma_{x}+\mbox{}^{0}X_{\mu}\sigma_{x}e_{i}-e_{j}\epsilon_{jki}\mbox{}^{k}X_{\mu}\sigma_{x}),\quad
i,j,k\in 1,2,3.$ (229)

Further considering the fact that among the 16 products
$\mbox{}^{i}X_{\mu}\sigma_{x}\mbox{}^{j}X_{\nu}$ ($i,j\in[0,3]$), only
$\mbox{}^{i}X_{\mu}\sigma_{x}\mbox{}^{0}X_{\nu}$ and
$\mbox{}^{0}X_{\mu}\sigma_{x}\mbox{}^{i}X_{\nu}$ ($i\in[1,3]$) are
nonvanishing, we readily obtain

$\displaystyle\mbox{}^{\mathrm{q}}\Xi^{i,\alpha\alpha}_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{q}}\xi^{\alpha{\dagger}}_{\mu}\mbox{}^{\mathrm{q}}\alpha_{i}\mbox{}^{\mathrm{q}}\xi^{\alpha}_{\nu}=\sum_{j=0}^{3}\mbox{}^{j}\Xi^{i}_{\mu\nu}e_{j},$
(230) $\displaystyle\mbox{}^{0}\Xi^{i,\alpha\alpha}_{\mu\nu}$
$\displaystyle=-h(\mbox{}^{i}X_{\mu}\sigma_{x}\mbox{}^{0}X_{\nu}-\mbox{}^{0}X_{\mu}\sigma_{x}\mbox{}^{i}X_{\nu}),\quad
i\in 1,2,3$ (231)
$\displaystyle=\frac{h}{2c}\begin{pmatrix}0&-g_{\mu}g^{i}_{\nu}\\\
g^{i}_{\mu}g_{\nu}&0\end{pmatrix},$ (232)
$\displaystyle\mbox{}^{j}\Xi^{i,\alpha\alpha}_{\mu\nu}$
$\displaystyle=h\epsilon_{jki}(\mbox{}^{k}X_{\mu}\sigma_{x}\mbox{}^{0}X_{\nu}+\mbox{}^{0}X_{\mu}\sigma_{x}\mbox{}^{k}X_{\nu}),\quad
i,j,k\in 1,2,3$ (233)
$\displaystyle=\frac{h}{2c}\epsilon_{jki}\begin{pmatrix}0&-g_{\mu}g^{k}_{\nu}\\\
-g^{k}_{\mu}g_{\nu}&0\end{pmatrix},\quad j\neq i.$ (234)

It can readily be verified that Eq. (230) agrees with Eq. (220).

In contrast, the density matrix $\mathbf{D}$ (36) is generally neither time-
reversal symmetric nor antisymmetric, but can be decomposed into a sum of
time-reversal symmetric ($\tilde{\mathbf{D}}$) and antisymmetric
($\mathbbm{i}\bar{\mathbf{D}}$) components, viz.,

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{D}$
$\displaystyle=\mathbf{Q}\mathbf{D}\mathbf{Q}^{\dagger}=\mbox{}^{\mathrm{Q}}\mathbf{C}\mathbf{n}\mbox{}^{\mathrm{Q}}\mathbf{C}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{D}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\mathbf{D}^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}\mathbf{D}^{\beta\alpha}&\mbox{}^{\mathrm{Q}}\mathbf{D}^{\beta\beta}\end{pmatrix}=\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}+\mathbbm{i}\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}},$
(235) $\displaystyle\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\beta\alpha}&\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\beta\beta}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\beta}\\\
-\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\beta*}&\mbox{}^{\mathrm{Q}}\tilde{\mathbf{D}}^{\alpha\alpha*}\end{pmatrix},$
(236) $\displaystyle\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\beta\alpha}&\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\beta\beta}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\beta}\\\
-\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\beta*}&\mbox{}^{\mathrm{Q}}\bar{\mathbf{D}}^{\alpha\alpha*}\end{pmatrix},$
(237)

such that

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{U}\mbox{}^{\mathrm{Q}}\mathbf{D}\mbox{}^{\mathrm{Q}}\mathbf{U}^{\dagger}=\begin{pmatrix}\mbox{}^{\mathrm{q}}\mathbf{D}&\mathbf{0}\\\
\mathbf{0}&\mbox{}^{\mathrm{q}}\mathbf{D}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}+h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}&\mathbf{0}\\\
\mathbf{0}&\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}+h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}\end{pmatrix}.$
(238)

Specific expressions for $\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}$ and
$\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}$ can be read out from Eq. (302) with
$\mbox{}^{\mathrm{Q}}\mathbf{D}$ for $\mathbf{M}$.

In terms of the above quantities, the $\mbox{}^{\mathrm{Q}}\mathbf{U}$
transformation (295) of the $\mathbf{Q}$-transformed DHF equation (31),

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{F}\mbox{}^{\mathrm{Q}}\mathbf{C}_{p}$
$\displaystyle=\mbox{}^{\mathrm{Q}}\mathbf{M}\mbox{}^{\mathrm{Q}}\mathbf{C}_{p}\epsilon_{p}\Leftrightarrow$
(239)
$\displaystyle\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{F}^{\alpha\alpha}&\mbox{}^{\mathrm{Q}}\mathbf{F}^{\alpha\beta}\\\
\mbox{}^{\mathrm{Q}}\mathbf{F}^{\beta\alpha}&\mbox{}^{\mathrm{Q}}\mathbf{F}^{\beta\beta}\end{pmatrix}\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{C}^{\alpha}_{p}\\\
\mbox{}^{\mathrm{Q}}\mathbf{C}^{\beta}_{p}\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{M}^{\alpha\alpha}&\mathbf{0}\\\
\mathbf{0}&\mbox{}^{\mathrm{Q}}\mathbf{M}^{\beta\beta}\end{pmatrix}\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{C}^{\alpha}_{p}\\\
\mbox{}^{\mathrm{Q}}\mathbf{C}^{\beta}_{p}\end{pmatrix}\epsilon_{p},\quad\epsilon_{p}\in\mathbb{R},$
(240)

gives rise to

$\displaystyle\mbox{}^{\mathrm{q}}\mathbf{F}\mbox{}^{\mathrm{q}}\mathbf{C}_{p}=\mbox{}^{\mathrm{q}}\mathbf{M}\mbox{}^{\mathrm{q}}\mathbf{C}_{p}\epsilon_{p},\quad\mbox{}^{\mathrm{q}}\mathbf{C}_{p}\in\mathbb{Q}^{2n\times
1},\quad\epsilon_{p}\in\mathbb{R},$ (241)

where $\mbox{}^{\mathrm{q}}F_{\mu\nu}$ has both time-reversal symmetric
($\mbox{}^{\mathrm{q}}\tilde{F}_{\mu\nu}$) and antisymmetric
($\mbox{}^{\mathrm{q}}\bar{F}_{\mu\nu}$) components, viz.,

$\displaystyle\mbox{}^{\mathrm{q}}F_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{q}}\tilde{F}_{\mu\nu}+\mbox{}^{\mathrm{q}}\bar{F}_{\mu\nu},$
(242) $\displaystyle\mbox{}^{\mathrm{q}}\tilde{F}_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{q}}h_{\mu\nu}+[(\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}|g_{12}|\mathrm{Tr}[\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Omega_{\kappa\lambda}])-(\mbox{}^{\mathrm{q}}\Omega_{\mu\lambda}|g_{12}|\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Omega_{\kappa\nu})]$
$\displaystyle-\mathrm{c_{g}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\lambda}|g_{12}|\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{i}_{\kappa\nu})-\mathrm{c_{b}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\lambda}|b_{12}^{ij}|\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{j}_{\kappa\nu})],$
(243) $\displaystyle\mbox{}^{\mathrm{q}}\bar{F}_{\mu\nu}$
$\displaystyle=\mathrm{c_{g}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}|g_{12}|\mathrm{Tr}[h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{i}_{\kappa\lambda}])+\mathrm{c_{b}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}|b_{12}^{ij}|\mathrm{Tr}[h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{j}_{\kappa\lambda}])$
$\displaystyle-(\mbox{}^{\mathrm{q}}\Omega_{\mu\lambda}|g_{12}|h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Omega_{\kappa\nu})-\mathrm{c_{g}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\lambda}|g_{12}|h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{i}_{\kappa\nu})$
$\displaystyle-\mathrm{c_{b}}(\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\lambda}|b_{12}^{ij}|h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}\mbox{}^{\mathrm{q}}\Xi^{j}_{\kappa\nu}),$
(244) $\displaystyle\mbox{}^{\mathrm{q}}M_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{Q}}M^{\alpha\alpha}_{\mu\nu}=\mbox{}^{\mathrm{Q}}M^{\beta\beta}_{\mu\nu}=\begin{pmatrix}S_{\mu\nu}&0\\\
0&\frac{1}{2c^{2}}T_{\mu\nu}\end{pmatrix}.$ (245)

Use of the fact that the time-reversal antisymmetric
($h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}$) and symmetric
($\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}$) parts of the density matrix
$\mbox{}^{\mathrm{q}}\mathbf{D}$ (238) do not contribute to the number density
(41) and current density (42), respectively, has been made to derive the Fock
matrix elements. Since the quaternion density matrix element
$\mbox{}^{\mathrm{q}}\tilde{\mathbf{D}}_{\lambda\kappa}$
($h\mbox{}^{\mathrm{q}}\bar{\mathbf{D}}_{\lambda\kappa}$) is generally
nonzero, each of the five nonzero scalar elements in
$\mbox{}^{\mathrm{q}}\Omega_{\mu\nu}$ (215) (3 in
$\mbox{}^{\mathrm{Q}}\Omega^{\alpha\alpha}_{\mu\nu}$ and 2 in
$\mbox{}^{\mathrm{Q}}\Omega^{\alpha\beta}_{\mu\nu}$) can ‘interact’ with each
of the five nonzero scalar elements in
$\mbox{}^{\mathrm{q}}\Omega_{\kappa\lambda}$, implying that there will be in
total 25 real-valued scalar integrals to evaluate and process just for a
single term of
$\mbox{}^{\mathrm{q}}G_{\mu\nu}[\mbox{}^{\mathrm{q}}\mathbf{D}]$ under the
Coulomb interaction[36]. On the other hand, each
$\mbox{}^{\mathrm{q}}\Xi^{i}_{\mu\nu}$ (220) has 18 nonzero scalar elements (6
in each of the three Cartesian components), there will be 324 real-valued
scalar integrals to evaluate and process just for a single term of
$\mbox{}^{\mathrm{q}}G_{\mu\nu}[\mbox{}^{\mathrm{q}}\mathbf{D}]$ under the
Gaunt or Breit interaction. Use of such observations has been made to achieve
very efficient implementations of the Gaunt[137] and Breit[138] integrals.

If wanted, the eigenvector $\mbox{}^{\mathrm{Q}}\mathbf{C}_{p}$ of Eq. (240)
can be mapped out from
$\mbox{}^{\mathrm{q}}\mathbf{C}_{p}=\sum_{i=0}^{3}\mbox{}^{i}\mathbf{C}_{p}$
(cf. Eq. (312)),

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{C}_{p}$
$\displaystyle=\psi(\mbox{}^{\mathrm{q}}\mathbf{C}_{p})\begin{pmatrix}\mathbf{I}_{2n}\\\
\mathbf{0}\end{pmatrix}=\begin{pmatrix}\mbox{}^{0}\mathbf{C}_{p}+\mathbbm{i}\mbox{}^{1}\mathbf{C}_{p}\\\
-\mbox{}^{2}\mathbf{C}_{p}+\mathbbm{i}\mbox{}^{3}\mathbf{C}_{p}\end{pmatrix}=\begin{pmatrix}\mbox{}^{\mathrm{Q}}\mathbf{C}^{\alpha}_{p}\\\
\mbox{}^{\mathrm{Q}}\mathbf{C}^{\beta}_{p}\end{pmatrix}.$ (246)

The quaternion density matrix $\mbox{}^{\mathrm{q}}\mathbf{D}$ can be
calculated as

$\displaystyle\mbox{}^{\mathrm{q}}\mathbf{D}$
$\displaystyle=\mbox{}^{\mathrm{q}}\mathbf{C}\mathbf{n}\mbox{}^{\mathrm{q}}\mathbf{C}^{\dagger}=\sum_{i=0}^{3}\mbox{}^{i}\mathbf{D}e_{i}\in\mathbb{Q}^{2n\times
2n},$ (247)

which is isomorphic to (cf. Eq. (301))

$\displaystyle\mbox{}^{\mathrm{Q}}\mathbf{D}=\psi(\mbox{}^{\mathrm{q}}\mathbf{D})=\mbox{}^{\mathrm{Q}}\mathbf{C}\mathbf{n}\mbox{}^{\mathrm{Q}}\mathbf{C}^{\dagger}\in\mathbb{C}^{4n\times
4n}.$ (248)

For a closed-shell system in the absence of external magnetic fields, the
time-reversal antisymmetric component $\mbox{}^{\mathrm{q}}\bar{\mathbf{F}}$
vanishes, such that
$\mbox{}^{\mathrm{q}}\mathbf{F}=\mbox{}^{\mathrm{q}}\tilde{\mathbf{F}}\in\mathbb{Q}^{2n\times
2n}$ can be diagonalized very efficiently[84, 139, 140, 86, 141], faster than
the diagonalization of
$\mbox{}^{\mathrm{Q}}\mathbf{F}=\mathbf{Q}\mathbf{F}\mathbf{Q}^{\dagger}=\mbox{}^{\mathrm{Q}}\tilde{\mathbf{F}}+\mbox{}^{\mathrm{Q}}\bar{\mathbf{F}}\in\mathbb{C}^{4n\times
4n}$ by up to a factor of two[141]. Yet, in the presence of
$\mbox{}^{\mathrm{q}}\bar{F}$, no gain in efficiency can be achieved in the
matrix diagonalization. Therefore, it is more appealing[36] to first construct
$\mbox{}^{\mathrm{q}}\tilde{\mathbf{F}}$ and
$\mbox{}^{\mathrm{q}}\bar{\mathbf{F}}$ and then map them to
$\mbox{}^{\mathrm{Q}}\tilde{\mathbf{F}}$ and
$\mbox{}^{\mathrm{Q}}\bar{\mathbf{F}}$ (cf. Eqs. (303) and (304)),
respectively.

Apart from the above generalized quaternion DHF equation (241), a Kramers-
unrestriced quaternion DHF (KUHF) scheme can also be formulated for open-shell
systems (see Ref. 142 for the two-component counterpart), by drawing analogy
with the spin-free UHF approach. The working equations read

$\displaystyle\mbox{}^{\mathrm{q}}\mathbf{F}^{\mathrm{x}}\mbox{}^{\mathrm{q}}\mathbf{C}^{\mathrm{x}}$
$\displaystyle=\mbox{}^{\mathrm{q}}\mathbf{M}\mbox{}^{\mathrm{q}}\mathbf{C}^{\mathrm{x}}\epsilon^{\mathrm{x}}_{p},\quad\mathrm{x}=\mathrm{u},\mathrm{d},$
(249) $\displaystyle\mbox{}^{\mathrm{q}}F^{\mathrm{x}}_{\mu\nu}$
$\displaystyle=\mbox{}^{\mathrm{q}}\tilde{F}_{\mu\nu}+\mbox{}^{\mathrm{q}}\bar{F}_{\mu\nu},$
(250) $\displaystyle\mbox{}^{\mathrm{q}}\tilde{F}_{\mu\nu}$