Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation Energy Functionals in the MGC-SDFT-UHFD Decoupling

Abstract

The Kohn-Sham formalism for the density functional theory (DFT) proposed a half-century ago has been the extensive motive force for the material science community, despite it is incomplete because of its problematic notion of eternally-unknown correlation energy functional including a separated part of kinetic energy. Here, we widely explain an alternative method recently discovered by us, i.e. the multiple grand canonical spin DFT (MGC-SDFT) in the unrestricted Hartree-Fock-Dirac (MGC-SDFT-UHFD) approximation. It is proved that the correlation energy functional consists of well-defined principal and secondary parts: the former yields the principal internal energy functional responsible for a set of the one-body quasi-particle spectra defined by the respective ground and excited states with each natural LCAO-MO as well as a set of the expected values of Heisenberg spin Hamiltonian, and the latter does a well-defined spin-dependent perturbation energy responsible for some many-body effects. An application will be made to explain why the water-splitting S1-state Mn4CaO5-clusters in photosystem II can exhibit two different EPR signals, called “g4.8” and “g12-multiline”. Moreover, the secondary correlation energy part will be shown to promote Cooper-pairings of Bloch-electrons near Fermi level in the superconductor, provided that their eigenstates might be exactly determined by the MGC-SDFT-UHFD method.

1. Introduction

In this chapter, we aim to explain why the predominant Kohn-Sham formalism of density functional theory (KS-DFT) based on the variational principle with respect to the electron density in a closed N-electron system [1, 2, 3, 4, 5, 6], must be stated as incomplete, during a number of active works motivated on it (e.g [7, 8, 9, 10, 11, 12]) still continuing, by pushing out the alternative electron density functional theory based on the multiple grand canonical quantum statistical variational principle capable of generating a large enough number of quantized energy levels of the ground and excited states in the unrestricted Hartree-Fock-Dirac approximation taking account of the explicit principal exchange-correlation energy functional -K_XC. This ultimate theory has been recently developed and called “multiple grand canonical spin DFT in the UHFD approximation (MGC-SDFT-UHFD method)” [13]. Moreover, we aim to present here a compact text of this ultimate MGC-SDFT-UHFD method in Sections 2.1 and 2.2 in order to help not only the reader’s understanding but also some program developers to challenge this painful-but-promising project to revise some codes associated with this paradigm shift from KS-DFT to MGC-SDFT-UHFD world in an extensive range of so far dedicated codes for predicting molecular and crystalline properties.

It is also important and exciting for us to be able to present as much as possible experimental evidence powerfully supporting the quantitative and systematic aspects of the MGC-SDFT-UHFD method to determine the one-body energy spectra, the quasi-particle’s wave functions, the magnetic property such as the mean isotropic spin-exchange coupling constants {J_i,j} and the total electronic internal energies, in Section 2.3. In [13], we provided the first experimental evidence for it; in Table 1, the derived formulas for J_1,2 demonstrated excellent quantitative agreements (less than 1% errors) with 10 experimental results from biomimetic binuclear transition metal complexes (TM: Cu, Mn, Fe), using the Mulliken’s atomic spin densities [14] and a set of the internal energies calculated by the UB3LYP/PBS/lacvp** method [15, 16]. Among many controversial problems that remained to be elucidated in photosynthesis research (see a recent review [17]), in Section 2.4 we discuss the second experimental supporting evidence provided by two broad EPR signals, named “g4.8” [18, 19] and “g12-multiline” [20], observed from the dark-stable S₁ state Mn₄CaO₅ clusters in the PSII having slightly different structures between thermophilic cyanobacteria in [18, 19] and higher-plant spinach in [20], respectively. At present, however, we have at hand only the structure of former’s PSII crystal at 1.95 Å high-resolution viewed by femtosecond XFEL pulse irradiation [21] but do not have any structure of the latter PSII crystal at least at similar high-resolution. It should be also noted that the super-brilliant femtosecond XFEL-pulse irradiation may generate high-density secondary photoelectrons to deoxidize nearby Mn₄ clusters with high probability during diffraction measurements. Then, the quantitative determination of the Heisenberg spin Hamiltonian involved in the principal exchange-correlation energy function can play a key role in the geometry optimization by the UB3LYP/PBS(ε)/lacvp** method to make the model Mn₄CaO_x cluster being thermally distributed in some isomeric substates of any Kok-S_i state.

Complexes	a; ${\mathrm{Cu}}_2^{\mathrm{II},\mathrm{II}}$	c; ${\mathrm{Mn}}_2^{\mathrm{IV},\mathrm{IV}}$	d; ${\mathrm{Mn}}_2^{\mathrm{III},\mathrm{IV}}$	e; ${\mathrm{Mn}}_2^{\mathrm{III},\mathrm{IV}}$	f; ${\mathrm{Mn}}_2^{\mathrm{IV},\mathrm{III}}$	g; ${\mathrm{Mn}}_2^{\mathrm{IV},\mathrm{IV}}$	h; ${\mathrm{Mn}}_2^{\mathrm{II},\mathrm{II}}$	i; ${\mathrm{Fe}}_2^{\mathrm{II},\mathrm{II}}$	j; ${\mathrm{Mn}}_2^{\mathrm{III},\mathrm{III}}$	k; ${\mathrm{Fe}}_2^{\mathrm{III},\mathrm{III}}$
NA	62	118	91	89	86	88	77	77	81	81
R1,2(Å)	2.659	2.745	2.591	2.591	3.230	2.296	3.370	3.315	3.140	3.202
bridge ligands	(μ−OAc)4	(μ−O)2	(μ−O)2(μ−OAc)	(μ−O)2(μ−OAc)	(μ−O)(μ−OAc)2	(μ−O)3	(μ−OH)(μ−OAc)2	(μ−OH)(μ−OAc)2	(μ−OH)(μ−OAc)2	(μ−OH)(μ−OAc)2
ε	10	20	20	20	40	20	40	20	20	10
${\mathrm{\Delta U}}_{\mathit{\text{UHFD}}}^{\left(2\right)}$(cm^−1)	−271.1	−1366.	−1579.	−1549.	−506.6	−3928.	−248.9	−223.8	141.4	−2522.
$n_{1\mathit{EF}}^{\left(1\right)}$; $n_{1\mathit{EF}}^{\left(2\right)}$	0.975; 0.980	3.064; 3.033	4.040; 4.031	4.087; 4.111	3.140; −3.22	3.106; 3.038	4.939; 4.912	3.913; 3.907	4.006; 3.882	4.686; 4.583
$n_{2\mathit{EF}}^{\left(1\right)}$; $n_{2\mathit{EF}}^{\left(2\right)}$	0.976; −0.981	3.107; −3.039	3.076; −3.189	3.034; −3.187	4.010; 4.010	3.103; −3.026	4.941; −4.911	3.912; −3.907	4.006; −3.886	4.683; −4.636
$S_{1\mathit{EF},2\mathit{EF}}^{\left(1\right)}$; $S_{1\mathit{EF},2\mathit{EF}}^{\left(2\right)}$; (×10^−2)	5.13; −3.76	−5.46; 2.38	−3.15; 5.50	−5.47; 8.21	−3.74; 5.69	−6.61; 2.13	2.44; −3.65	4.52; −4.80	0.280; −6.06	13.94; 16.46
$J_{1,2}^{\left(4,-,0\right)}$; (cm^−1)	−284.9	−143.5	−126.5	−123.4	−40.1	−407.6	−10.2	−14.6	8.8	−114.9
$J_{1,2}^{\mathit{exp}}$; (cm^−1)	−285	−144	−125	−125	−40	−407	−9	−14	9	−115

Table 1.

Benchmark-test results of the 2GC-UHFD-SD averaged ES-exchange coupling constants, designated $J_{1,2}^{\left(4,-,0\right)}$, for 10 biomimetic binuclear Cu, Mn and Fe complexes, using the UB3LYP/PBS(ε)/lacvp** (4^th XC) method combined with the inherent formulas of Eqs. (103), (113), derived in the present MGC-UHFD-SDFT (0^th XC) method. These model TM₂ complexes consist of N_A atoms, are imbedded in each ε dielectric constant medium, and exhibit a variety of TM1-TM2 distances, designated R_{1, 2}, which depend strongly on different bridge structures and the different TM valences, and also weakly on the different non-bridging ligations of paramagnetically-polarized O atoms and diamagnetically-/paramagnetically-polarized N atoms (not shown here). Here the data for b (di-μ-oxo bridged Mn^IV-Mn^IV dimer ligated by four picolinic anions) in Table I in [13] was omitted because of its optimized structure containing no solvent molecules. (see Supplemental Online Material of Ref. [13] about how to calculate the effective spin densities using the Mulliken’s atomic spin densities [14]).

Furthermore, in Section 2.5, we discuss the most interesting many-body effect induced by the secondary correlation energy term, which represents a spin-dependent attractive correlation interaction between a couple of conductive Bloch-electrons with antiparallel spins that could be generated only near the Fermi surface in the metallic crystal. This strong correlation interaction may accelerate the phase transition from the normal state to the superconductive state by promoting Cooper-pairings of conductive Bloch-electrons near the Fermi level in the superconductor against the common knowledge [22, 23, 24, 25, 26, 27].

A problematic idea underlining the KS-DFT formalism may be described in other words such that the ground state energy E of the one-particle self-consistent field Hamiltonian for N electron systems, which corresponds to the internal energy functional of the electron density determined in thermal equilibrium state, should be further minimized by “the exact variation principle” with respect to the electron density regarded as a variational variable to search for “the exact energy functional” of “the exact electron density”, subjected to the N-representability condition [4]. This wrong variational idea appears to have been widely accepted so far, although it may have been enforced by a special situation enforced by too strong expressions involving many exact’s: “the exact variational principle”, and “the existence theorem of an exact energy functional of the exact electron density” as well as “N-representability condition”. Especially, “the N-representability condition” seems to be too strong to consider any open quantum system, in which the total numbers (N_α, N_β) of (up, down) electrons in the system should be replaced by the mean values (${\overline{N}}_{\alpha }$, ${\overline{N}}_{\beta })$ of a pair of the expected values of their operators, (${\widehat{N}}_{\alpha }$, ${\widehat{N}}_{\alpha }),\text{respectively}$, in the context of applying the variational principle to the minimum grand potential including them, as will be shown in this chapter. So far, neither theoretical proof nor evidence for the K-S formalism could not be provided unless the exact correlation energy function is discovered.

In a GC ensemble, one may consider a much larger M-electron system (M » N) of atoms, molecules, and solids, which will be maximally realized with a finite probability in contact with a grand canonical heat/particle reservoir containing a much larger number of electrons at temperature θ/k_B (k_B is the Boltzmann constant). All the stational states of the M-electron system, which involve the ground and all kinds of excited states, may be assumed to be describable in terms of the time-independent non-relativistic Schrödinger wave equation in 3D-space:

$$H\mathrm{\Psi }=E\mathrm{\Psi }\left({\mathit{x}}_1,{\mathit{x}}_2,\cdots \cdots 、{\mathit{x}}_M\right);M»N,$$E1

where H is the Hamiltonian operator given by

$$H=T+V_{\mathit{ne}}+V_{\mathit{ee}},$$E2

$$=\sum _{i=1}^M\left(-\frac{1}{2}{\nabla }_i^2\right)+\sum _{i=1}^M\upsilon \left({\mathit{r}}_i,\mathrm{\epsilon }\right)+\sum _{i<j}^M\frac{1}{r_{i,j}}.$$E3

Here, x_i ≡ (r_i, s_i) represents the (orbital, spin) coordinates of the ith electron, T is the electron kinetic energy operator; V_ne is the electrostatic interaction operator of electrons with all nuclei and the surrounding medium of the dielectric constant ε if a convenient “the linear Poisson-Boltzmann equation Solver (PBS)” model [8] is augmented; V_ee is the electron Coulomb interaction operator; and r_i,j ≡ |r_i-r_j|. Since it is impossible to exactly solve Eq. (1) except for the case of a hydrogen atom, we have developed the ultimate MGC-SDFT formalism [7], which has been constructed by developing five new methodological concepts in Subsections 2.1 through 2.5 along the basic principles of quantum thermodynamics with the theory of open quantum systems, but not of closed quantum system as adopted in the Kohn-Sham formalism.

2. Multiple grand canonical spin density functional theory

2.1 Definition of a grand canonical ensemble: One-particle and two-particle reduced density matrices

The principally most general choice would be made for an extended antisymmetric Slater determinant wave function as the trial many-electron wave function $\mathrm{\Psi }$(x₁, x₂,…, x_M) in the Schrödinger Eq. (1), consisting of a complete large enough number, (M_α, M_β), of mutually-independent natural molecular spin-orbitals (NMO) wave functions, Eq. (6), which may be most-appropriately expandable in terms of a Linear Combination of gaussian-type or Slater-type Atomic Orbitals (LCAO), G_A≡ [$g_l^a$(r)] (a, the atomic order number from 1 to N_A; l, a set of plural AO quantum numbers). The maximum size of (M_α, M_β)-dimensional Hilbert space of the NMO set M′ of Eq. (4) must be as large enough as possible to satisfy the near-completeness condition of Eq. (5) and the orthonormality relations of Eqs. (6) and (7) as far as the highest NMO energy levels may not exceed a dissociation limit given by a work function, W.

$${\mathbf{M}}^{\prime }={}^{\alpha }{\mathbf{M}}^{\prime }+{}^{\beta }{\mathbf{M}}^{\prime }\equiv \left\{{\mathrm{\Psi }}_{\mathrm{\tau }}^{\prime }\left(\mathbf{x}\right),ϵ{\prime }_{\tau };\tau =1_{\mathrm{\alpha }},\dots ,M_{\alpha },1_{\mathrm{\beta }},\dots ,M_{\beta }\right\},$$E4

$$\sum _{m\mathrm{\sigma }}^{{}^{\sigma }{\mathbf{M}}^{\prime }}\left|\mathrm{\Psi }{\prime }_{\mathit{m\sigma }}⟩\right.⟨\left.\mathrm{\Psi }{\prime }_{\mathit{m\sigma }}\right|\cong \mathbf{1}\text{for}{M}_{\sigma }>>N_{\mathrm{\sigma }};\mathrm{\sigma }=\mathrm{\alpha },\mathrm{\beta }.$$E5

$$\mathrm{\Psi }{\prime }_{\mathit{m\sigma }}\left(\mathit{x}\right)\equiv \mathrm{\Phi }{\prime }_{\mathit{m\sigma }}\left(\mathit{r}\right)\xi {\prime }_{\mathit{m\sigma }}\left(\mathit{s}\right),$$E6

$$⟨\mathrm{\Psi }{\prime }_{\mathit{m\sigma }}|\mathrm{\Psi }{\prime }_{\mathit{m\sigma }\prime \prime }⟩=⟨\mathrm{\Phi }{\prime }_{\mathit{m\sigma }}|\mathrm{\Phi }{\prime }_{\mathit{m\sigma }\prime \prime }⟩⟨\xi {\prime }_{\mathit{m\sigma }}|\xi {\prime }_{\mathit{m\sigma }\prime \prime }⟩={\delta }_{\mathit{m\sigma },\mathit{m\sigma }\prime \prime }{\delta }_{\sigma ,\sigma \prime \prime },$$E7

$${\mathrm{\Phi }}_{\mathit{m\sigma }}\left(\mathit{r}\right)=\sum _{a,l}^{{\mathbf{G}}_{\mathbf{A}}}{C}_{a,l}^{\mathit{m\sigma }}{g}_l^a\left(\mathit{r}\right)=\sum _{a=1}^{N_A}\sum _{/}{C}_{a,l}^{\mathit{m\sigma }}{g}_l^a\left(\mathit{r}\right).$$E8

In Eq. (8), the AO basis functions, which may include polarization and diffuse functions, are assumed to be orthonormalized in each atom such as $⟨g_l^a|g_{l\prime }^a⟩={\delta }_{l,l\prime }$ but slightly overlap between the valence-electron orbitals of neighboring atoms, except for most AOs in the core levels. Hereafter, the single dashed quantities in Eqs. (4)–(7) will be used for the quantities in a thermal nonequilibrium state. The completeness Eq. (5) and the orthonormality Eqs. (6) and (7) are assumed to seamlessly hold even in such non-equilibrium states.

It is important to remind that the thermodynamic equilibrium state can be achieved in terms of the Rayleigh-Ritz variational principle applied to the one-particle reduced density matrix given by

$$\widehat{\mathrm{\Gamma }}\left[p^{\prime },{\mathrm{\Psi }}^{\prime }\right]\equiv \sum _{N\prime }^{\infty }\sum _{\tau }^{\mathbf{M}\prime }P{\prime }_{N\prime ,\tau }⌊\mathrm{\Psi }{\prime }_{N\prime ,\tau }⟩⟨\left.\mathrm{\Psi }{\prime }_{N\prime ,\tau }\right|,$$E9

$$=\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}⌊\mathrm{\Psi }{\prime }_{\tau ,n_{\tau }}⟩⟨\left.\mathrm{\Psi }{\prime }_{\tau ,n_{\tau }}\right|,$$E10

with respect to a set of the distribution probabilities, designated $\left\{P_{N\prime ,\tau }^{\prime };N^{\prime }=0,\dots ,\infty ;\tau =1,\dots ,M\right\}$ for fermions and bosons in Eq. (9) or $\left\{p{\prime }_{\tau ,}{n}_{\tau };\tau =1,\dots ,\mathrm{M};n_{\tau }=0,1\right\}$ for (NMO-transformed) fermions in Eq. (10). This nonequilibrium state will relax to the maximum entropy state keeping the normalization condition of Eq. (11) and the binary chemical potential (μ_α/μ_β) equilibrium conditions with the heat/particle reservoir leading to Eqs. (12) and (13):

$$\mathrm{Tr}\left[\widehat{\mathrm{\Gamma }}\right]=\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}=1,$$E11

$${\widehat{N}}_{\sigma }=\sum _{\tau \left(\sigma \right)}^{{}^{\sigma }\mathbf{M}\prime }{\widehat{a}}_{\tau \left(\sigma \right)}^{†}{\widehat{a}}_{\tau \left(\sigma \right)};\mathrm{\sigma }=\mathrm{\alpha },\mathrm{\beta }$$E12

$$\mathrm{Tr}\left[{\widehat{N}}_{\sigma }\widehat{\mathrm{\Gamma }}\right]=\sum _{\tau \left(\sigma \right)}^{{}^{\mathit{\sigma }}{\mathbf{M}}^{\prime }}\sum _{n_{\tau \left(\sigma \right)}=0}^1{n}_{\tau \left(\sigma \right)}p{\prime }_{\tau \left(\sigma \right),n_{\tau \left(\sigma \right)}}\equiv {\overline{N}}_{\sigma }\left[p^{\prime },{\mathrm{\Psi }}^{\prime }\right],$$E13

where ${\widehat{a}}_{\tau }$ and ${\widehat{a}}_{\tau }^{†}$are the annihilation and creation operators of an electron in the τth NMO-eigenstate, respectively. Subsequently, using the GC entropy

$$\mathrm{\Sigma }\left[p^{\prime },{\mathrm{\Psi }}^{\prime }\right]=-k_B\mathrm{Tr}\left(\widehat{\mathrm{\Gamma }}ln\widehat{\mathrm{\Gamma }}\right)=-k_B\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}lnp{\prime }_{\tau ,n_{\tau }},$$E14

and the second quantization expression of the Hamiltonian operator, that is

$$\widehat{H}={\widehat{H}}_0+{\widehat{H}}_1=\sum _{\tau }^{\mathbf{M}\prime }ϵ{\prime }_{\tau }{\widehat{a}}_{\tau }^{†}{\widehat{a}}_{\tau }+{\widehat{H}}_1,$$E15

which in general consists of the principal part ${\widehat{H}}_0\text{and}$ the secondary part ${\widehat{H}}_1$, responsible for the GC ensemble of mutually independent NMO fermions and a perturbational interaction between them, respectively, we obtain the grand potential in thermal nonequilibrium state:

$$\mathrm{\Omega }\left[p^{\prime },{\mathrm{\Psi }}^{\prime }\right]\equiv \mathrm{Tr}\left[\widehat{\mathrm{\Gamma }}\left(\theta ln\widehat{\mathrm{\Gamma }}+{\widehat{H}}_0-{\mu }_{\alpha }{\widehat{N}}_{\alpha }-{\mu }_{\beta }{\widehat{N}}_{\beta }\right)\right],$$E16

$$=\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}\left[\theta lnp{\prime }_{\tau ,n_{\tau }}+⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|\left(ϵ{\prime }_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right)n_{\tau }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩\right].$$E17

In Eqs. (16) and (17) should be noted that $\widehat{H}\mathit{\text{is replaced}}\mathit{by}\mathit{\text{the principal part}}{\widehat{H}}_0=\sum _{\tau }^{\mathbf{M}\prime }ϵ{\prime }_{\tau }{\widehat{a}}_{\tau }^{†}{\widehat{a}}_{\tau }.$ The variational equations subject to the normalization condition Eq. (11) are given by

$$\frac{\partial }{\partial p{\prime }_{\tau ,n_{\tau }}}\left(\mathrm{\Omega }\left[\left\{p{\prime }_{\tau ,n_{\tau }},\mathrm{\Psi }{\prime }_{\tau }\right\}\right]+\mathrm{\lambda }\left[\sum _{\tau }^{\mathbf{M}}\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}\mathrm{ー}1\right]\right)=0,$$E18

with λ being a Lagrange’s multiplier, satisfying

$$e^{-1-\lambda /\theta }=\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^1⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|ϵ{\prime }_{\tau }-{\mu }_{\alpha }{\overline{N}}_{\alpha }-{\mu }_{\beta }{\overline{N}}_{\beta }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩.$$E19

Thus, we obtain an intermediate solution with fixed {Ψ′_τ}:

$$p{\prime }_{\tau ,n_{\tau }}\stackrel{\mathit{\text{yields}}}{\to }{p}_{\tau ,n_{\tau }}^0\left({\mathrm{\Psi }}^{\prime }\right)=\frac{exp\left[⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|\left({\mu }_{\sigma \left(\tau \right)}-ϵ{\prime }_{\tau }\right)n_{\tau }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩/\theta \right]}{{\sum }_{n_{\tau =0}}^{1}exp\left[⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|\left({\mu }_{\sigma \left(\tau \right)}-ϵ{\prime }_{\tau }\right)n_{\tau }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩/\theta \right]}$$E20

$${\widehat{\mathrm{\Gamma }}}^0\left({\mathrm{\Psi }}^{\prime }\right)=\sum _{\tau }^{\mathbf{M}\prime }\sum _{n_{\tau }=0}^1{p}_{\tau ,n_{\tau }}^0\left({\mathrm{\Psi }}^{\prime }\right)⌊\mathrm{\Psi }{\prime }_{\tau ,n_{\tau }}⟩⟨\left.\mathrm{\Psi }{\prime }_{\tau ,n_{\tau }}\right|.$$E21

Then, the GC potential decreases by the non-negative quantity (the equality appears when p′ = p⁰), namely

$$\begin{array}{c}\mathrm{\Omega }\left[\left\{p{\prime }_{\tau ,n_{\tau }},\mathrm{\Psi }{\prime }_{\tau }\right\}\right]-\mathrm{\Omega }\left[\left\{p_{\tau ,n_{\tau }}^0,\mathrm{\Psi }{\prime }_{\tau }\right\}\right]\\ =\sum _{\tau }^{\mathbf{M}}\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}\left(\theta lnp{\prime }_{\tau ,n_{\tau }}+⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|\left(ϵ{\prime }_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right)n_{\tau }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩\right)\\ +\theta ln\sum _{\tau }^{\mathbf{M}}\sum _{n_{\tau }=0}^{1}exp\left[⟨\mathrm{\Psi }{\prime }_{\tau },n_{\tau }|\left({\mu }_{\sigma \left(\tau \right)}-ϵ{\prime }_{\tau }\right)n_{\tau }|\mathrm{\Psi }{\prime }_{\tau },n_{\tau }⟩/\theta \right]\\ =\theta \sum _{\tau }^{\mathbf{M}}\sum _{n_{\tau }=0}^{1}p{\prime }_{\tau ,n_{\tau }}\left[lnp{\prime }_{\tau ,n_{\tau }}-ln{p}_{\tau ,n_{\tau }}^0\left({\mathrm{\Psi }}^{\prime }\right)\right]\ge 0,\end{array}$$E22

which is notably induced only by the increase of the GC entropy Σ [p′,Ψ′] -Σ [p⁰,Ψ′]. This initial-guess state of the GC ensemble M′ = {Ψ′_τ, $ϵ{\prime }_{\tau }$} will relax to converge toward the self-consistent, orthonormal and complete eigenvectors/eigenvalues set M of Eqs. (23)–(27) by iteration technique, which will be described in Section 2.2.

$$\mathbf{M}={}^{\alpha }\mathbf{M}+{}^{\beta }\mathbf{M}\equiv \left\{{\mathrm{\Psi }}_{\mathrm{m\sigma }}\left(\mathbf{x}\right),{ϵ}_{m\mathrm{\sigma }},f_{\mathit{m\sigma }}\left({\mu }_{\sigma }\right);\mathit{m\sigma }\equiv \tau =1_{\mathrm{\alpha }},\dots ,{\mathrm{M}}_{\mathrm{\alpha }},1_{\mathrm{\beta }},\dots ,{\mathrm{M}}_{\mathrm{\beta }}\right\},$$E23

$$\sum _{m\mathrm{\sigma }}^{{}^{\sigma }\mathbf{M}}\left|{\mathrm{\Psi }}_{\mathit{m\sigma }}⟩\right.⟨\left.{\mathrm{\Psi }}_{\mathit{m\sigma }}\right|\cong \mathbf{1}\text{for}{\mathrm{M}}_{\mathrm{\sigma }}>>{\mathrm{N}}_{\mathrm{\sigma }};\mathrm{\sigma }=\mathrm{\alpha },\mathrm{\beta }.$$E24

$${\mathrm{\Psi }}_{\mathit{m\sigma }}\left(\mathit{x}\right)\equiv {\mathrm{\Phi }}_{\mathit{m\sigma }}\left(\mathit{r}\right){\xi }_{\mathit{m\sigma }}\left(\mathit{s}\right),$$E25

$$⟨{\mathrm{\Psi }}_{\mathit{m\sigma }}|{\mathrm{\Psi }}_{\mathit{m\sigma }\prime }⟩=⟨{\mathrm{\Phi }}_{\mathit{m\sigma }}|{\mathrm{\Phi }}_{\mathit{m\sigma }\prime }⟩⟨{\xi }_{\mathit{m\sigma }}|{\xi }_{\mathit{m\sigma }\prime }⟩={\delta }_{\mathit{m\sigma },\mathit{m\sigma }\prime \prime }{\delta }_{\sigma ,\sigma \prime \prime },$$E26

$$f_{\mathit{m\sigma }}\left({\mu }_{\sigma }\right)=p_{\mathit{m\sigma },1}={\left\{exp\left[\left({ϵ}_{\mathit{m\sigma }}-{\mu }_{\sigma }\right)/\theta \right]+1\right\}}^{-1}\equiv f\left({ϵ}_{\mathit{m\sigma }}-{\mu }_{\sigma }\right),$$E27

where it should be noted that only the populated distribution probability p⁰_mσ,1 of Eq. (20) converges to the Fermi-Dirac distribution probability $f_{\mathit{m\sigma }}\left({\mu }_{\sigma }\right)$ of Eq. (27).

Next, we need to introduce the first-order (for one-particle interactions) and the second order (for two-particle interactions) reduced electron density matrixes, Γ(x, x′) and $\mathrm{\Gamma }$₂(x₁x₂, x′₁x′₂), respectively, for the GC ensemble M in thermal equilibrium state, as given by

$$\mathrm{\Gamma }\left(\mathit{x},{\mathit{x}}^{\prime }\right)=\sum _{\tau }^{\mathbf{M}}{f}_{\tau }\left|{\mathrm{\Psi }}_{\tau }\left(\mathit{x}\right)⟩\right.⟨\left.{\mathrm{\Psi }}_{\tau }\left({\mathit{x}}^{\prime }\right)\right|,$$E28

$${\mathrm{\Gamma }}_2\left({\mathit{x}}_1{\mathit{x}}_2,{\mathit{x}}_1^{\prime }{\mathit{x}}_2^{\prime }\right)=\frac{1}{2}\left[\mathrm{\Gamma }\left({\mathit{x}}_1,\mathit{x}{\prime }_1\right)\mathrm{\Gamma }\left({\mathit{x}}_2,\mathit{x}{\prime }_2\right)-\mathrm{\Gamma }\left({\mathit{x}}_1,\mathit{x}{\prime }_2\right)\mathrm{\Gamma }\left({\mathit{x}}_2,\mathit{x}{\prime }_1\right)\right],$$E29

which will be used in Section 2.2.

2.2 The self-consistent field method in the UHFD approximation

Here, we derive the self-consistent field (SCF) method to generate such a realistic GC ensemble, M, as given by Eqs. (23)–(27), in which all NMO levels will be partially occupied with the Fermi-Dirac distribution probability f(${\mathit{ϵ}}_{\mathit{m\sigma }}\mathit{-}$μ_σ) constrained by the chemical potentials μ_σ defined by either Eq. (27) or the Gibbs free energy per a σ-spin electron for σ = (α, β), as given by Eq. (30), using the mean number of σ-spin NMO-fermions for σ = (α, β) in Eq. (31):

$$\mathrm{G}={\mu }_{\alpha }{\overline{N}}_{\alpha }+{\mu }_{\beta }{\overline{N}}_{\beta },$$E30

$${\overline{N}}_{\sigma }=\sum _{m\mathrm{\sigma }}^{{}^{\sigma }\mathbf{M}}f\left({ϵ}_{\mathit{m\sigma }}-{\mu }_{\sigma }\right).$$E31

At first, we define various spinless electron density matrixes for later usage:

$${\rho }_{\sigma }\left(\mathit{r},{\mathit{r}}^{\prime }\right)={\mathrm{Tr}}_{\mathit{s}}{\left[{\mathrm{\Gamma }}_{\sigma }\left(\mathit{x},{\mathit{x}}^{\prime }\right)\right]}_{\mathit{s}=\mathit{s}\prime }=\sum _{\mathit{m\sigma }}^{{}^{\sigma }\mathbf{M}}{f}_{\mathit{m\sigma }}\left({\mu }_{\sigma }\right){\mathrm{\Phi }}_{\mathit{m\sigma }}\left(\mathit{r}\right){\mathrm{\Phi }}_{\mathit{m\sigma }}^{\ast }\left({\mathit{r}}^{\prime }\right);\mathrm{\sigma }=\mathrm{\alpha },\mathrm{\beta },$$E32

$$\rho \left(\mathit{r},{\mathit{r}}^{\prime }\right)={\rho }_{\alpha }\left(\mathit{r},{\mathit{r}}^{\prime }\right)+{\rho }_{\beta }\left(\mathit{r},{\mathit{r}}^{\prime }\right),$$E33

$${\rho }_{\sigma }\left(\mathit{r}\right)\equiv {\rho }_{\sigma }\left(\mathit{r},\mathit{r}\right),\rho \left(\mathit{r}\right)={\rho }_{\alpha }\left(\mathit{r}\right)+{\rho }_{\beta }\left(\mathit{r}\right),$$E34

where Tr_s represents the trace on the spin coordinate s.

Next, we will prove that the internal energy $U\left[{\mathrm{\Gamma }}_{\mathrm{\alpha }},{\mathrm{\Gamma }}_{\mathrm{\beta }}\right]$ as a function of the reduced density matrix Γ = (Γ_α, Γ_β) can be decoupled into two parts, as seen in Eq. (36): (1) the principal part $U_{\mathit{\text{UHFD}}}^0\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]$ of Eq. (37) including the principal exchange-correlation energy functional -K_XC [ρ_α, ρ_β] defined by Eq. (40) and (2) a secondary part containing only a spin-dependent correlation energy functional $∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]$defined by Eq. (43). This ultimate decoupling scheme neglecting the secondary correlation term $∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}[\mathrm{\Gamma }$_α, $\mathrm{\Gamma }$_β] of Eq. (43) will be tentatively called “the Unrestricted Hartree-Fock-Dirac (UHFD) approximation”. because the Dirac’s spin permutation operator (σ_i is called spinor)

$$P_{12}^{\sigma }=\frac{1}{2}\left[1+\left({\mathit{\sigma }}_1,{\mathit{\sigma }}_2\right)\right]=\frac{1}{2}\left[1+4\left({\mathit{s}}_1,{\mathit{s}}_2\right)\right],$$E35

including an inner product of two of Pauli’s spin operators $\left({\mathrm{\sigma }}_1,{\mathrm{\sigma }}_2\right)=\left(2{\mathit{s}}_1,2{\mathit{s}}_2\right),$has played a decisive role in our discovery of this new decoupling scheme. This is indeed a revolutionary discovery a long way beyond the early Unrestricted Hartree (UH), Unrestricted Hartree-Fock (UHF) and Unrestricted Hartree-Fock-Slater (UHFS) approximations. This new UHFD decoupling scheme leads to a group of fundamental equations:

$$U\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]=U_{\mathit{\text{UHFD}}}^0\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]+∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right],$$E36

$$U_{\mathit{\text{UHFD}}}^0\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]\equiv \mathrm{T}\left[\rho \right]+V_{\mathit{ne}}\left[\rho \right]+J\left[\rho \right]-\frac{1}{2}{K}_{\mathit{XC}}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]+H_{\mathit{ES}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right],$$E37

$$V_{\mathit{ne}}\left[\rho \right]=\int d\mathit{r}\upsilon \left(\mathit{r},\mathrm{\epsilon }\right)\rho \left(\mathit{r}\right),$$E38

$$J\left[\rho \right]=\frac{1}{2}\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}\rho \left({\mathit{r}}_1\right)\rho \left({\mathit{r}}_2\right),$$E39

$$-K_{\mathit{XC}}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]=-\frac{1}{2}\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}\rho \left({\mathit{r}}_1,{\mathit{r}}_2\right)\rho \left({\mathit{r}}_2,{\mathit{r}}_1\right),$$E40

$$H_{\mathit{ES}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]=-2\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}Q\left({\mathit{r}}_1,{\mathit{r}}_2\right)Q\left({\mathit{r}}_2,{\mathit{r}}_1\right)\left({\mathit{s}}_1,{\mathit{s}}_2\right),$$E41

$$Q\left({\mathit{r}}_1,{\mathit{r}}_2\right)={\rho }_{\alpha }\left({\mathit{r}}_1,{\mathit{r}}_2\right)-{\rho }_{\beta }\left({\mathit{r}}_1,{\mathit{r}}_2\right),$$E42

$$\begin{array}{c}∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]=-2\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}[{\rho }_{\alpha }\left({\mathit{r}}_1,{\mathit{r}}_2\right){\rho }_{\beta }\left({\mathit{r}}_2,{\mathit{r}}_1\right)\\ +{\rho }_{\beta }\left({\mathit{r}}_1,{\mathit{r}}_2\right){\rho }_{\alpha }\left({\mathit{r}}_2,{\mathit{r}}_1\right)]\left({\mathit{s}}_1,{\mathit{s}}_2\right),\end{array}$$E43

$\text{where}$T [ρ] represents the expected value of the electron kinetic energy operator T, although this notation does not mean any explicit functional form of ρ, the other explicit energy functionals of ρ have usual meanings, and $H_{\mathit{ES}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]$ does the spin density coupling energy functional between two NMO-fermions, which is expected to contain H_ex. (CORREGENDUM: Please add a miss-dropped factor, 2, in Eq. (3.18) in [13], just like above Eq. (41)).

[Proof of Eqs. (36)–(42)] Substituting Eq. (28) into Eq. (29), we get the NMO expansion formula of the exchange-product matrix:

$$\begin{array}{c}\mathrm{\Gamma }\left({\mathit{x}}_1,\mathit{x}{\prime }_2\right)\mathrm{\Gamma }\left({\mathit{x}}_2,\mathit{x}{\prime }_1\right)=\sum _{\sigma ,\mathit{m\sigma }}^{\mathbf{M}}{f}_{\mathit{m\sigma }}\left({\mu }_{\sigma }\right){\mathrm{\Phi }}_{\mathit{m\sigma }}\left({\mathit{r}}_1\right){\mathrm{\Phi }}_{\mathit{m\sigma }}{\left(\mathit{r}{\prime }_2\right)}^{\ast }\\ \times \sum _{\sigma \prime ,\mathit{m\sigma }\prime }^{\mathbf{M}}{f}_{\mathit{m\sigma }\prime }\left({\mu }_{\sigma \prime }\right){\mathrm{\Phi }}_{\mathit{m\sigma }\prime }\left({\mathit{r}}_2\right){\mathrm{\Phi }}_{\mathit{m\sigma }}{\left(\mathit{r}{\prime }_1\right)}^{\ast }\left|{\xi }_{\mathit{m\sigma }}\left({\mathit{s}}_2\right){\xi }_{\mathit{m\sigma }\prime }\left({\mathit{s}}_1\right)⟩\right.⟨\left.{\xi }_{\mathit{m\sigma }}\left(\mathit{s}{\prime }_1\right){\xi }_{\mathit{m\sigma }\prime }\left(\mathit{s}{\prime }_2\right)\right|.\end{array}$$E44

Here, to restore the spin-pair wave function to the normal-order form, Dirac’s spin operator of Eq. (35) needs to be operated to the two-spin function:

$${\xi }_{\mathit{m\sigma }}\left({\mathit{s}}_2\right){\xi }_{\mathit{m\sigma }\prime }\left({\mathit{s}}_1\right)=P_{12}^{\sigma }{\xi }_{\mathit{m\sigma }}\left({\mathit{s}}_1\right){\xi }_{\mathit{m\sigma }\prime }\left({\mathit{s}}_2\right).$$E45

Then, using Eqs. (32), (33), (35), and (45), we can transform Eq. (44) into two different formulas:

$$\begin{array}{c}\mathrm{\Gamma }\left({\mathit{x}}_1,\mathit{x}{\prime }_2\right)\mathrm{\Gamma }\left({\mathit{x}}_2,\mathit{x}{\prime }_1\right)\\ =\rho \left({\mathit{r}}_1,\mathit{r}{\prime }_2\right)\rho \left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)\left[\frac{1}{2}\left(1+{\sigma }_{1z}{\sigma }_{2z}\right)+\frac{1}{4}\left({\sigma }_{1+}{\sigma }_{2-}+{\sigma }_{1-}{\sigma }_{2+}\right)\right],\end{array}$$E46

$$\begin{array}{c}=\frac{1}{2}\left[{\rho }_{\alpha }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\alpha }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)+{\rho }_{\beta }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\beta }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)\right]\\ +\frac{1}{2}\left[{\rho }_{\alpha }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\beta }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)+{\rho }_{\beta }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\alpha }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)\right]\\ +\frac{1}{2}Q\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right)Q\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)\left({\mathit{\sigma }}_1,{\mathit{\sigma }}_2\right)\\ +\frac{1}{2}\left[{\rho }_{\alpha }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\beta }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)+{\rho }_{\alpha }\left({\mathit{r}}_1,\mathit{r}{\prime }_2\right){\rho }_{\beta }\left({\mathit{r}}_2,\mathit{r}{\prime }_1\right)\right]\left({\mathit{\sigma }}_1,{\mathit{\sigma }}_2\right),\end{array}$$E47

with the use of the spin density matrix Q of Eq. (42) and the off-diagonal spinors ${\sigma }_{j\pm }$of Eq. (48):

$${\sigma }_{j\pm }={\sigma }_{\mathit{jx}}\pm {\mathit{i\sigma }}_{\mathit{jy}},j=1,2.$$E48

Apparently, there exist two decoupling schemes: (1) In Eq. (46) is decoupled a pair of off-diagonal spinor terms, leading to the UHF approximation, and (2) In Eq. (47) decoupled only the last term, leading to the UHFD approximation. In the UHFD approximation, we obtain the functional formula for the internal energy $U\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]$ in 3D-spin space decomposed into the principal part $U_{\mathit{\text{UHFD}}}^0\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]$ and the secondary part ΔE_UHFD^corr in Eq. (36). (QED)

Since all the DFT calculations can be made in the binary-spin Hilbert space, we must take the trace of H_ES [Γ_α, Γ_β] in Eq. (37) on the 3D-spin coordinates to obtain its functional of ρ(ρ_α, ρ_β), which is equal to a half of the principal exchange-correlation energy functional, that is

$$H_{\mathit{ES}}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]\equiv {\mathrm{Tr}}_{\mathit{S}}\left(H_{\mathit{ES}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]\right)=-\frac{1}{2}{K}_{\mathit{XC}}\left[{\rho }_{\alpha },{\rho }_{\beta }\right].$$E49

Using Eq. (49), we also obtain the principal internal energy functional of ρ:

$$U_{\mathit{\text{UHFD}}}^0\left[{\rho }_{\alpha },{\rho }_{\beta }\right]\equiv {\mathrm{Tr}}_{\mathit{S}}\left(U_{\mathit{\text{UHFD}}}^0\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]\right)=\mathrm{T}\left[\rho \right]+V_{\mathit{ne}}\left[\rho \right]+J\left[\rho \right]-K_{\mathit{XC}}\left[{\rho }_{\alpha },{\rho }_{\beta }\right],$$E50

which should be equated to the GC ensemble average of the principal part of the Hamiltonian operator in the second quantization representation in the thermal equilibrium state (see Eq. (15)), that is

$${\widehat{H}}_{\mathit{\text{UHFD}}}^0=\sum _{\tau }^{\mathbf{M}}{ϵ}_{\tau }{\widehat{a}}_{\tau }^{†}{\widehat{a}}_{\tau },$$E51

leading to

$$U_{\mathit{\text{UHFD}}}^0\left[{\rho }_{\alpha },{\rho }_{\beta }\right]\equiv {\mathrm{Tr}}_{\mathbf{M}}\left[\mathrm{\Gamma }{\widehat{H}}_0\right]=\int d\mathit{r}\sum _{\tau }^{\mathbf{M}}{ϵ}_{\tau }{f}_{\tau }\left({ϵ}_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right){\mathrm{\Phi }}_{\tau }^{\ast }\left(\mathit{r}\right){\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right).$$E52

Similarly, $U_{\mathit{\text{UHFD}}}^0$can be expanded as the GC ensemble average of a self-consistent effective Hamiltonian as given by

$$\begin{array}{c}{U}_{\mathit{\text{UHFD}}}^0\left[{\rho }_{\alpha },{\rho }_{\beta }\right]=\int d\mathit{r}\left\{\sum _{\tau }^{\mathbf{M}}{f}_{\tau }\left({ϵ}_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right){\mathrm{\Phi }}_{\tau }^{\ast }\left(\mathit{r}\right)\right.\\ \times \left.\left[-\frac{1}{2}{\nabla }^2{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right)+\int {\upsilon }_{\mathit{NMO}}\left(\mathit{r},{\mathit{r}}^{\prime }\right){\mathrm{\Phi }}_{\tau }\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }\right]\right\}\end{array},$$E53

with the use of the local and non-local NMO-based effective potential defined by

$${\upsilon }_{\mathit{NMO}}\left(\mathit{r},{\mathit{r}}^{\prime }\right)=\left[\upsilon \left(\mathit{r},\epsilon \right)+\frac{1}{2}\int d{\mathit{r}}^{″}\frac{\rho \left({\mathit{r}}^{″}\right)}{\left|\mathit{r}-{\mathit{r}}^{″}\right|}\right]\delta \left(\mathit{r}-{\mathit{r}}^{\prime }\right)-\frac{\rho \left(\mathit{r},{\mathit{r}}^{\prime }\right)}{2\left|\mathit{r}-{\mathit{r}}^{\prime }\right|},$$E54

$$\rho \left(\mathit{r}\right)=\sum _{\tau }^{\mathbf{M}}{f}_{\tau }\left({ϵ}_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right){\left|{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right)\right|}^2,$$E55

$$\rho \left(\mathit{r},{\mathit{r}}^{\prime }\right)=\sum _{\tau }^{\mathbf{M}}{f}_{\tau }\left({ϵ}_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right){\mathrm{\Phi }}_{\tau }^{\ast }\left({\mathit{r}}^{\prime }\right){\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right).$$E56

From equivalent Eqs. (52) and (53) leading to Eq. (57), we obtain a series of central Schrödinger equations, (58) by putting […]_τ = 0:

$$\begin{array}{c}{U}_{\mathit{\text{UHFD}}}^0\left[{\rho }_{\alpha },{\rho }_{\beta }\right]-{\mathrm{Tr}}_{\mathbf{M}}\left[\mathrm{\Gamma }{\widehat{H}}_0\right]=\int d\mathit{r}\left\{\sum _{\tau }^{\mathbf{M}}{f}_{\tau }\left({ϵ}_{\tau }-{\mu }_{\sigma \left(\tau \right)}\right){\mathrm{\Phi }}_{\tau }^{\ast }\left(\mathit{r}\right)\right.\\ \left.\left[-\frac{1}{2}{\nabla }^2{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right)+\int {\upsilon }_{\mathit{NMO}}\left(\mathit{r},{\mathit{r}}^{\prime }\right){\mathrm{\Phi }}_{\tau }\left({\mathit{r}}^{\prime }\right)d{\mathit{r}}^{\prime }-{ϵ}_{\tau }{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right)\right]\right\}=0\end{array}$$E57

$$\begin{array}{c}-\frac{1}{2}{\nabla }^2{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right)+\int d{\mathit{r}}^{\prime }{\upsilon }_{\mathit{NMO}}\left(\mathit{r},{\mathit{r}}^{\prime }\right){\mathrm{\Phi }}_{\tau }\left({\mathit{r}}^{\prime }\right)={ϵ}_{\tau }{\mathrm{\Phi }}_{\tau }\left(\mathit{r}\right),\\ \text{and}\mathrm{h}.\mathrm{c}.\text{for}\mathrm{all}\tau \in \mathbf{M}.\end{array}$$E58

This central (not variational!) solution will be stocked as the presumed GC ensemble M, in which the eigenvalues $\left\{{ϵ}_{\tau }\right\}$ are usually assumed to be rearranged from the minimum ${ϵ}_{1\sigma }$ to the maximum ${ϵ}_{\mathrm{M\sigma }}$ in the order of increasing energies, first for α-spin NMOs and second for β-spin NMOs, as

$$\tau =1_{\alpha },2_{\alpha },\cdots ,M_{\alpha },1_{\beta },2_{\beta },\cdots ,M_{\beta }=1,2,\cdots ,M\left(M=M_{\alpha }+M_{\beta }\right).$$E59

For simplicity, we assume that there exists no degeneracy in energy levels in the unrestricted large system without any structural symmetry.

On the other hand, the secondary correlation energy functional, $∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}[\mathrm{\Gamma }$_α, $\mathrm{\Gamma }$_β], defined by Eq. (43), represents the sole perturbation term in 3D-spin space. Taking the trace of it on the 3D-spin coordinates, we obtain the second quantization expression of it as follows

$$\begin{array}{c}{\widehat{H}}_{\mathit{\text{UHFD}}}^1={\mathrm{Tr}}_s\left(∆E_{\mathit{\text{UHFD}}}^{\mathit{\text{corr}}}\left[{\mathrm{\Gamma }}_{\alpha },{\mathrm{\Gamma }}_{\beta }\right]\right)\\ =-\frac{1}{2}\iint \frac{{d\mathit{r}}_1{d\mathit{r}}_2}{r_{12}}{\mathrm{Tr}}_{\mathit{s}}\left\{\left[{\rho }_{\alpha }\left({\mathit{r}}_1,{\mathit{r}}_2\right){\rho }_{\beta }\left({\mathit{r}}_2,{\mathit{r}}_1\right)+{\rho }_{\beta }\left({\mathit{r}}_1,{\mathit{r}}_2\right){\rho }_{\alpha }\left({\mathit{r}}_2,{\mathit{r}}_1\right)\right]\right.\left.\left({\mathit{\sigma }}_1,{\mathit{\sigma }}_2\right)\right\}\\ =-\sum _{\mathit{m\alpha }}^{{}^{\alpha }\mathbf{M}}\sum _{\mathit{m\beta }}^{{}^{\beta }\mathbf{M}}f\left({ϵ}_{\mathit{m\alpha }}-{\mu }_{\alpha }\right)f\left({ϵ}_{\mathit{m\beta }}-{\mu }_{\beta }\right)\iint \frac{{d\mathit{r}}_1{d\mathit{r}}_2}{r_{12}}\\ \times \left[{\mathrm{\Phi }}_{\mathit{m\alpha }}^{\ast }\left({\mathit{r}}_2\right){\mathrm{\Phi }}_{\mathit{m\alpha }}\left({\mathit{r}}_1\right){\mathrm{\Phi }}_{\mathit{m\beta }}^{\ast }\left({\mathit{r}}_1\right){\mathrm{\Phi }}_{\mathit{m\beta }}\left({\mathit{r}}_2\right)\right]\left({\sigma }_{+,\mathit{m\alpha },}{\sigma }_{-,\mathit{m\beta }}+{\sigma }_{-,\mathit{m\alpha }}{\sigma }_{+,\mathit{m\beta }}\right)\\ =-\sum _{\mathit{m\alpha }}^{{}^{\alpha }\mathbf{M}}\sum _{\mathit{m\beta }}^{{}^{\beta }\mathbf{M}}\left({\sigma }_{+,\mathit{m\alpha },}{\sigma }_{-,\mathit{m\beta }}+{\sigma }_{-,\mathit{m\alpha }}{\sigma }_{+,\mathit{m\beta }}\right)\iint \frac{{d\mathit{r}}_1{d\mathit{r}}_2}{r_{12}}\\ \times \left[{\mathrm{\Phi }}_{\mathit{m\alpha }}^{\ast }\left({\mathit{r}}_2\right){\mathrm{\Phi }}_{\mathit{m\beta }}^{\ast }\left({\mathit{r}}_1\right){\mathrm{\Phi }}_{\mathit{m\beta }}\left({\mathit{r}}_2\right){\mathrm{\Phi }}_{\mathit{m\alpha }}\left({\mathit{r}}_1\right)\right]{\widehat{a}}_{\mathit{m\beta }}^{†}{\widehat{a}}_{\mathit{m\beta }}{\widehat{a}}_{\mathit{m\sigma }}^{†}{\widehat{a}}_{\mathit{m\alpha }}.\end{array}$$E60

The first-order perturbation term vanishes owing to the nondiagonal spinors. However, the second-order perturbation correction can always induce a finite attractive force between any pair of NMO-fermions with antiparallel spins. The most interesting example would be a positive enhancement effect on the Cooper-pair superconductivity due to an additional attractive force between two conductive Bloch-electrons with antiparallel spins near Fermi level, as will be discussed in Section 2.4.

2.3 MGC-SDFT-UHFD method for polynuclear transition metal complexes

We next consider a variety of paramagnetic systems including plural n (≥ 2) spins, designated {S_i, i = 1,…,n}, which arise from transition metal (TM) cations, C/N/O-radicals, -C=C- bond radicals, and so on. These spins are quantum-mechanically interacting with each other via the exchange coupling constants (J_{i, j}) in the Heisenberg spin Hamiltonian defined by

$$H_{\mathit{ex}}=-2\sum _{i=1}^{n-1}\sum _{j=i+1}^n{J}_{i,j}\left({\mathit{S}}_i,{\mathit{S}}_j\right).$$E61

However, this H_ex model takes account only the pure spin operators {S_i} but does not contain any kind of polarized spins of ligand atoms, designated {s_iL}, “Why?” The most fundamentally important question is, “What is the origin of H_ex?” These questions have been recently solved by Kusunoki [13], as reviewed in this subsection. Let us investigate what kinds of spin-dependent physical processes are involved in the spin-density coupling energy functional H_ES [Γ_α,Γ_β] of Eqs. (41) and (42).

Since in the binary spin space appear 2ⁿ⁻¹ (n ≥ 2) mutually-independent up/down-iES arrangements (ESA), which represent one ferromagnetic and the other anti-ferromagnetic states, we must prepare a set of multiple grand canonical (MGC) ensembles, as given by

$$\mathbf{M}=\sum _{k=1}^{2^{n-1}}{\mathbf{M}}^{\left(k\right)},{\mathbf{M}}^{\left(k\right)}={}^{\alpha }{\mathbf{M}}^{\left(k\right)}+{}^{\beta }{\mathbf{M}}^{\left(k\right)};$$E62

$${{}^{\mathit{\sigma }}\mathbf{M}}^{\left(k\right)}\equiv \left\{{\mathrm{\Psi }}_{\mathit{m\sigma }}^{\left(k\right)}\left(\mathbf{x}\right),{ϵ}_{\mathit{m\sigma }}^{\left(k\right)},f_{\mathit{m\sigma }}^{\left(k\right)}\left({ϵ}_{\mathit{m\sigma }}^{\left(k\right)}-{\mu }_{\sigma }^{\left(k\right)}\right);\mathit{m\sigma }=1,\dots ,M_{\sigma }\right\};\mathrm{\sigma }=\mathrm{\alpha },\mathrm{\beta }.$$E63

Practically, we can calculate only a set of 2ⁿ⁻¹ principal internal energy functionals:

$$U_{\mathit{\text{UHFD}}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]=T^{\left(k\right)}\left[\rho \right]+V_{\mathit{ne}}^{\left(k\right)}\left[\rho \right]+J^{\left(k\right)}\left[\rho \right]-K_{\mathit{XC}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right];k=1,\dots ,2^{n-1}.$$E64

However, the origin of H_ex must be traced to a set of 2ⁿ⁻¹ equality relationships between the principal exchange-correlation energy functional, $-K_{\mathit{XC}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]$, and the projected value of the spin-dependent XC energy functional, i.e.

$$-K_{\mathit{XC}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]=2H_{\mathit{ES}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]=2{\mathrm{Tr}}_s\left(H_{\mathit{ES}}^{\left(k\right)}\left[{\mathrm{\Gamma }}_{\alpha }^{\left(k\right)},{\mathrm{\Gamma }}_{\beta }^{\left(k\right)}\right]\right);k=1,\dots ,2^{n-1}.$$E65

We note that the kth projected value of $\left({\mathit{S}}_i,{\mathit{S}}_j\right)$ onto the binary Hilbert space spanned by the kth GC ensemble M^(k) must depend not only on the kth principal exchange-correlation energy between iES and jES, but also on the polarized spins of bridging and non-bridging ligand atoms, iLaj (j ≠ i) and iLnb, respectively, via the conservation law of each projected spin number of n_i^(k) and n_j^(k), defined by

$$n_i^{\left(k\right)}\equiv 2S_i{\sigma }_{z,i}^{\left(k\right)};{\sigma }_{z,i}^{\left(k\right)}=\pm 1,k=1,\dots ,2^{n-1};i=1,\dots ,n.$$E66

$$\sum _{i=1}^n{S}_i{\sigma }_{z,i}^{\left(k\right)}\ge 0,k=1,\dots ,2^{n-1}.$$E67

Although the arrangement order of n ESA-codes $\left\{{\sigma }_{z,i}^{\left(k\right)}=1\text{or}-1,i=1-n\right\}$ leaves the choice of one’s best depending on the arrangement order of n TM-cation spins {S_i, −1 -- n}, a non-negative sum rule of Eq. (67) must be satisfied owing to the time-reversal symmetry.

For the ith transition metal (TM) cations, its non-bridging ligand atoms iLnb and its bridging ligand atoms iLaj between the ith and jth TM-cations, this wave-packet spin projection may be entrusted to the respective spin operator by itself, that is S_i, s_iLnb and s_iLaj, by imposing each projection equation acting on a spin-dependent AO in a LCAO-NMO wave function, without change of the F-D distribution for the former two and with change to its half distribution for the latter one:

$${\mathit{S}}_i{g}_l^a\left(\mathit{x}\right)={\delta }_{a,i}{\delta }_{l,3d_i}{\mathit{S}}_i{g}_{3d}^i\left(\mathit{x}\right),$$E68

$${\mathit{s}}_{\mathit{\text{iLnb}}}{g}_l^a\left(\mathit{x}\right)={\delta }_{a\in \mathit{\text{iLnb}}}{\delta }_{l,2p}{\mathit{s}}_{\mathit{\text{iLnb}}}{g}_{2p}^{\mathit{\text{iLnb}}}\left(\mathit{x}\right),$$E69

$${\mathit{s}}_{\mathit{\text{iLaj}}}{g}_l^a\left(\mathit{x}\right)={\delta }_{a\in \mathit{\text{iLaj}}}{\delta }_{l,2p}{\mathit{s}}_{\mathit{\text{iLaj}}}{g}_{2p}^{\mathit{\text{iLaj}}}\left(\mathit{x}\right);f_{\mathit{\text{iLajmσ}}}=\frac{f_{\mathit{im\sigma }}}{{\nu }_{\mathit{\text{iLaj}}}}=\frac{f_{\mathit{im\sigma }}}{2},$$E70

where we have introduced the share frequency ν_iLajmσ among n different subsets (in this case, it is 2), δ_{a ∈ A} = (1 for a ∈ A; 0 for otherwise) and δ_x,y is Kronecker’s δ.

Concomitantly, the kth GC ensemble M^(k) might be decomposed into n spin-dependent NMO-subsets associated with these elements, and the other spin-free subset as in Eq. (71), each ${\mathbf{M}}_i^{\left(k\right)}$ further decomposed into three components as in Eqs. (72) and (73), finally to define the ith ES in terms of two components in Eq. (74).

$${\mathbf{M}}^{\left(k\right)}=\sum _{i=1}^n{\mathbf{M}}_i^{\left(k\right)}+\mathbf{M}{\prime }^{\left(k\right)},$$E71

$${\mathbf{M}}_i^{\left(k\right)}={\mathbf{M}}_{\mathit{id}}^{\left(k\right)}+{\mathbf{M}}_{\mathit{iL}}^{\left(k\right)}\mathbf{+}{\mathbf{M}}_{\mathit{io}}^{\left(k\right)},$$E72

$${\mathbf{M}}_{\mathit{iL}}^{\left(k\right)}=\sum _{j\ne i}^n{\mathbf{M}}_{\mathit{\text{iLaj}}}^{\left(k\right)}+{\mathbf{M}}_{\mathit{\text{iLnb}}}^{\left(k\right)},$$E73

$${\mathbf{M}}_{\mathit{iES}}^{\left(k\right)}={\mathbf{M}}_{\mathit{id}}^{\left(k\right)}+\sum _{j\ne i}^n{\mathbf{M}}_{\mathit{\text{iLaj}}}^{\left(k\right)},$$E74

where the ith subset in Eq. (72) consists of the subset ${\mathbf{M}}_{\mathit{id}}^{\left(k\right)}$ associated with (3d, 4 s, 4p)-electron AO’s in the ith TM cation, $\text{the subset}{\mathbf{M}}_{\mathit{iL}}^{\left(k\right)}$ associated with 2p-valence electron AO’s in the iLth assembly of ligand atoms, and $\text{the subset}{\mathbf{M}}_{\mathit{io}}^{\left(k\right)}$ associated with other doubly-occupied core-shell AO’s in the ith TM cation. The iLth subset in Eq. (73) can be further decomposed into two kinds of ligand assembly: (1) a thermal equipartition half of the iLajth assembly of bridging ligand atoms between the ith and jth TM cations, which can mediate the antiferromagnetic super-exchange coupling, and (2) the iLnbth assembly of non-bridging ligand atoms around the ith TM cation, which can be either paramagnetically or diamagnetically polarized depending on the ligand C/N/O atomic structure and hence control the ith ES magnitude via the spin number (n_i) conservation law governing the Mulliken atomic spin densities {$M_a^{\left(k\right)}$} [14] of these magnetically-interacting atoms, finally to define the ith ES density $n_{\mathit{iEF}}^{\left(k\right)}$ so as to satisfy Eq. (76):

$$n_{\mathit{iEF}}^{\left(k\right)}\equiv M_{\mathit{id}}^{\left(k\right)}+\sum _{j\ne i}^n\sum _{\mathit{\text{iLaj}}}\frac{M_{\mathit{\text{iLaj}}}^{\left(k\right)}\vartheta \left(M_{\mathit{\text{iLaj}}}^{\left(k\right)}{\sigma }_{z,i}^{\left(k\right)}\right)}{{\nu }_{\mathit{\text{iLaj}}}^{\left(k\right)}}+{\delta }_{k>1}\sum _{j\ne i}^n\sum _{\mathit{\text{iLaj}}}{\sigma }_{z,i}^{\left(k\right)}\frac{M_{\mathit{\text{iLaj}}}^{\left(1\right)}\vartheta \left(-{\sigma }_{z,j}^{\left(1\right)}{\sigma }_{z,i}^{\left(1\right)}\right)}{{\nu }_{\mathit{\text{iLaj}}}^{\left(1\right)}};$$E75

$$n_{\mathit{\text{iLnb}}}^{\left(k\right)}\equiv \sum _{\mathit{\text{iLnb}}\prime }{M}_{\mathit{\text{iLnb}}\prime }^{\left(k\right)};n_{\mathit{iEF}}^{\left(k\right)}+n_{\mathit{\text{iLnb}}}^{\left(k\right)}=n_i.k=1\left(\mathrm{F}\right),2\left(\mathrm{AF}\right),\dots ,2^{n-1}\left(\mathrm{AF}\right).$$E76

Here, ϑ(sign1*sign2) is the Heaviside step function, and ${\nu }_{\mathit{\text{iLaj}}}^{\left(k\right)}$ is the frequency of the iLaj spin density shared by plural ESs with the same sign, which can be calculated as

$${\nu }_{\mathit{\text{iLaj}}}^{\left(k\right)}=1+\vartheta \left({\sigma }_{z,i}^{\left(k\right)}{\sigma }_j^{\left(k\right)}\right)\text{for}\mathrm{all}k’\mathrm{s},$$E77

(CORRIGENDAM: Eq. (5.6c) in [13] should be replaced by Eq. (77)).

Decomposition into these NMO-subsets allows us to provide a noble systematic and quantitative method to derive a set of the expected values of H_ex {<H_ex > ^(k); k = 1, …, 2ⁿ⁻¹} from the spin-dependent XC energy functional in Eq. (65), as follows:

The kth spin density matrix for Eq. (42) can be decomposed into

$$Q^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)=\sum _{\sigma =\alpha }^{\beta }{\left(-1\right)}^{\sigma }\sum _{\mathit{m\sigma }}^{{}^{\sigma }{\mathbf{M}}^{\left(\mathit{k}\right)}}{f}_{\mathit{m\sigma }}{\mathrm{\Phi }}_{\mathit{m\sigma }}\left(\mathit{r}\right){\mathrm{\Phi }}_{\mathit{m\sigma }}{\left(\mathit{r}\right)}^{\ast };{\left(-1\right)}^{\alpha /\beta }=\pm 1,$$E78

$$=\sum _{i=1}^n\left[Q_{\mathit{iES}}^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)+Q_{\mathit{\text{iLnb}}}^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)+Q_{\mathit{io}}^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)\right]+Q{\prime }^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right),$$E79

$$\begin{array}{c}=\sum _{i=1}^n\left[Q_{\mathit{iES}}^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)+Q_{\mathit{\text{iLnb}}}^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)\right]\\ +\sum _{i=1}^n\left[{\rho }_{\mathit{io},\alpha }^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)-{\rho }_{\mathit{io},\beta }^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)\right]+\rho {\prime }_{\alpha }^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right)-\rho {\prime }_{\beta }^{\left(k\right)}\left(\mathit{r},{\mathit{r}}^{\prime }\right),\end{array}$$E80

in which Eq. (80) indicates that the first term will contribute to both intra-atomic and interatomic spin-density coupling energies generated by open 3d-shell electrons, as given by

$$2H_{\mathit{ES}1}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]\equiv 2E_{\mathit{ES}1}^{\left(k\right)}=-\sum _{i=1}^n{J}_{i,i}^{\left(k\right)}{S}_i\left(S_i+1\right),$$E81

$$J_{i,i}^{\left(k\right)}\equiv \frac{1}{2n_i^2}\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}\left[{\rho }_{\mathit{iES}}^{\left(k\right)}\left({\mathit{r}}_1,{\mathit{r}}_2\right)+{\rho }_{\mathit{\text{iLnb}}}^{\left(k\right)}\left({\mathit{r}}_1,{\mathit{r}}_2\right)\right]\left[{\rho }_{\mathit{iES}}^{\left(k\right)}\left({\mathit{r}}_2,{\mathit{r}}_1\right)+{\sigma }_{\mathit{\text{iLnb}}}^{\left(k\right)}\left({\mathit{r}}_2,{\mathit{r}}_1\right)\right],$$E82

and

$$2H_{\mathit{ES}2}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]=2E_{\mathit{ES}2}^{\left(k\right)}={⟨H_{\mathit{ex}}⟩}^{\left(k\right)}=-2\sum _{i<j}^n{J}_{i,j}{⟨\left({\mathit{S}}_i,{\mathit{S}}_j\right)⟩}^{\left(k\right)},$$E83

$${⟨\left({\mathit{S}}_i,{\mathit{S}}_j\right)⟩}^{\left(k\right)}=\mathrm{Tr}\left[\left({\mathit{S}}_i,{\mathit{S}}_j\right){\mathrm{\Gamma }}_2^{\left(k\right)}\left({\mathit{x}}_1{\mathit{x}}_2,{\mathit{x}}_1{\mathit{x}}_2\right)\right],$$E84

$$=\sum _{M_i=-S_i}^{S_i}\sum _{M_j=-S_j}^{S_j}\mathrm{Tr}\left[⟨M_i,M_j|\left({\mathit{S}}_i,{\mathit{S}}_j\right)|M_i,M_j⟩\right.\times \left.⟨M_i,M_j|{\mathrm{\Gamma }}_{i,j}^{\left(k\right)}\left({\mathit{x}}_1{\mathit{x}}_2,{\mathit{x}}_1{\mathit{x}}_2\right)|M_i,M_j⟩\right],$$E85

with

$${\mathrm{\Gamma }}_{i,j}^{\left(k\right)}\left({\mathit{x}}_1{\mathit{x}}_2,{\mathit{x}}_1{\mathit{x}}_2\right)=\frac{2}{2}\left[{\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left({\mathit{x}}_1,{\mathit{x}}_1\right){\mathrm{\Gamma }}_{\mathit{jES}}^{\left(k\right)}\left({\mathit{x}}_2,{\mathit{x}}_2\right)-{\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left({\mathit{x}}_1,{\mathit{x}}_2\right){\mathrm{\Gamma }}_{\mathit{jES}}^{\left(k\right)}\left({\mathit{x}}_2,{\mathit{x}}_1\right)\right],$$E86

respectively, and the second and third terms to a major H_ex-less.

component in the principal XC energy, as given by

$$-K_{\mathit{XC}-H_{\mathit{ex}}}^{\left(k\right)}\left[{\rho }_{\alpha }^{\left(k\right)},{\rho }_{\beta }^{\left(k\right)}\right]=2H_{\mathit{ES}0}^{\left(k\right)}\left[{\rho }_{\mathit{ES}0,\alpha }^{\left(k\right)},{\rho }_{\mathit{ES}0,\beta }^{\left(k\right)}\right]=2E_{\mathit{ES}0}^{\left(k\right)},$$E87

$$\begin{array}{c}=-\frac{1}{2}\iint \frac{d{\mathit{r}}_{1}d{\mathit{r}}_2}{r_{12}}\left[{\rho }^{\left(k\right)}\left({\mathit{r}}_1,{\mathit{r}}_2\right)-\sum _{i=1}^n\left\{{\rho }_{\mathit{iES}}^{\left(k\right)}\left({\mathit{r}}_1,{\mathit{r}}_2\right)+{\rho }_{\mathit{\text{iLnb}}}^{\left(k\right)}\left({\mathit{r}}_1,{\mathit{r}}_2\right)\right\}\right]\\ \times \left[{\rho }^{\left(k\right)}\left({\mathit{r}}_2,{\mathit{r}}_1\right)-\sum _{i=1}^n\left\{{\rho }_{\mathit{iES}}^{\left(k\right)}\left({\mathit{r}}_2,{\mathit{r}}_1\right)+{\rho }_{\mathit{\text{iLnb}}}^{\left(k\right)}\left({\mathit{r}}_2,{\mathit{r}}_1\right)\right\}\right];\end{array}$$E88

To calculate Eq. (84), we need the total spin operator given by

$${\mathit{S}}_{\mathit{tot}}=\sum _{i=1}^n\left({\mathit{S}}_i+\sum _{\mathit{\text{iLaj}}}{\mathit{s}}_{\mathit{\text{iLaj}}}+\sum _{\mathit{\text{iLnb}}}{\mathit{s}}_{\mathit{\text{iLnb}}}\right),$$E89

The ith ES spin operator, S_i, is considered to turn around between up and down states, |S_i > and |-S_i>, in the binary Hilbert space spanned by two rules:

$$\sum _{M_i=-S_i}^{S_i}\left|M_i⟩\right.⟨\left.M_i\right|={\mathbf{1}}_i;i=1,\dots ,n,$$E90

and

$$⟨{\alpha }_i|M_i⟩=\frac{{\delta }_{M_i,S_i}}{\sqrt{\left|n_i\right|}},⟨\beta |M_i⟩=\frac{{\delta }_{M_i,-S_i}}{\sqrt{\left|n_i\right|}};$$E91

while the iLajth and iLnbth ES spin operators, s_iLaj and s_iLnb, are assumed to automatically respond to the up or down state of S_i. Then, the expected value of the z-component of S_tot in the kth ESA state is given by

$${⟨S_{\mathit{tot},z}⟩}^{\left(k\right)}=\sum _{i=1}^n\mathrm{Tr}\left\{S_{i,z}\left[{\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left(\mathit{x},\mathit{x}\right)+{\mathrm{\Gamma }}_{\mathit{\text{iLnb}}}^{\left(k\right)}\left(\mathit{x},\mathit{x}\right)\right]\right\},$$E92

$$=\sum _{i=1}^n{S}_i{\sigma }_{z,i}^{\left(k\right)}\left(P_{\mathit{iES}}^{\left(k\right)}+P_{\mathit{\text{iLnb}}}^{\left(k\right)}\right)=\sum _{i=1}^n{S}_i{\sigma }_{z,i}^{\left(k\right)},$$E93

$${\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left(\mathit{x},{\mathit{x}}^{\prime }\right)={\mathrm{\Gamma }}_{\mathit{id}}^{\left(k\right)}\left(\mathit{x},{\mathit{x}}^{\prime }\right)+\sum _{j\ne i}^n\sum _{\mathit{\text{iLaj}}\prime }{\mathrm{\Gamma }}_{\mathit{\text{iLaj}}\prime }^{\left(k\right)}\left(\mathit{x},{\mathit{x}}^{\prime }\right),$$E94

$$P_{\mathit{iES}}^{\left(k\right)}=\frac{n_{\mathit{iES}}^{\left(k\right)}}{n_i}=\frac{1}{n_i}\int d\mathit{r}{Q}_{\mathit{iES}}^{\left(k\right)}\left(\mathit{r},\mathit{r}\right),$$E95

$$P_{\mathit{\text{iLnb}}}^{\left(k\right)}=\frac{n_{\mathit{\text{iLnb}}}^{\left(k\right)}}{n_i}=\frac{1}{n_i}\int d\mathit{r}{Q}_{\mathit{\text{iLnb}}}^{\left(k\right)}\left(\mathit{r},\mathit{r}\right),$$E96

$$P_{\mathit{iES}}^{\left(k\right)}+P_{\mathit{\text{iLnb}}}^{\left(k\right)}=\left[n_{\mathit{iES}}^{\left(k\right)}+n_{\mathit{\text{iLnb}}}^{\left(k\right)}\right]/n_i=1.$$E97

It is important to remind that the last additional term in Eq. (75) in the antiferromagnet kAFth ESA state was absolutely required for a systematic better agreement with experimental J_i,j indicating that a large-positive $M_{\mathit{\text{iLaj}}}^{\left(1F\right)}$density in the 1Fth-ESA state can be divided into two half densities which will be reversed to be distributed with phase matching to an antiferromagnetic pair of $n_{\mathit{iES}}^{\left(\mathit{kAF}\right)}$ and $n_{\mathit{jES}}^{\left(\mathit{kAF}\right)}$ in the kAFth ESA state, as given by a programmatic equation [13],

$$M_{\mathit{iL}\mathrm{a}j}^{\left(k\right)}\Leftarrow \frac{1}{2}{\sigma }_{z,i}^{\left(k\right)}{M}_{\mathit{\text{iLaj}}}^{\left(1\right)}\vartheta \left(-{\sigma }_{z,i}^{\left(k\right)}{\sigma }_{z,j}^{\left(k\right)}\right)+M_{\mathit{\text{iLaj}}}^{\left(k\right)}.$$E98

Now, substituting Eq. (86) into Eq. (85) and taking the traces over spin-orbital coordinates with use of Eqs. (90), (91), (95) and (96), we obtain a noble formula

$${⟨\left({\mathit{S}}_i,{\mathit{S}}_j\right)⟩}^{\left(\mathrm{k}\right)}=n_{\mathit{iES}}^{\left(k\right)}{n}_{\mathit{jES}}^{\left(k\right)}\left(1-S_{\mathit{iES},\mathit{jES}}^{\left(k\right)}\right),$$E99

$$S_{\mathit{iES},\mathit{jES}}^{\left(k\right)}\equiv \frac{\mathrm{Tr}\left[{\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left({\mathit{x}}_1,{\mathit{x}}_2\right){\mathrm{\Gamma }}_{\mathit{jES}}^{\left(k\right)}\left({\mathit{x}}_2,{\mathit{x}}_1\right)\right]}{\mathrm{Tr}\left[{\mathrm{\Gamma }}_{\mathit{iES}}^{\left(k\right)}\left({\mathit{x}}_1,{\mathit{x}}_1\right)\right]\mathrm{Tr}\left[{\mathrm{\Gamma }}_{\mathit{jES}}^{\left(k\right)}\left({\mathit{x}}_1,{\mathit{x}}_1\right)\right]},k=1,\dots ,2^{n-1}.$$E100

Here, {$S_{\mathit{iES},\mathit{jES}}^{\left(k\right)}$} represents a set of the Exchange-Correlation vs. Classical Coulomb Density Overlap Integral (XC/CC-DOI) ratios. Although it appears almost impossible to directly calculate a set of 2ⁿ⁻²n(n-1) XC/CC-DOI ratios, we could find out a reasonable solution of Eq. (103) by imposing 2ⁿ⁻¹ equations to eliminate all the residues {${∆}^2{⟨{\mathit{S}}_{\mathit{tot}}^2⟩}^{\left(k\right)}$} from a set of the expected values of {${⟨{\mathit{S}}_{\mathit{tot}}^2⟩}^{\left(k\right)}$}, which are given by

$${⟨{\mathit{S}}_{\mathit{tot}}^2⟩}^{\left(k\right)}={⟨S_{\mathit{tot},z}⟩}^{\left(k\right)}\left[{⟨S_{\mathit{tot},z}⟩}^{\left(k\right)}+1\right]+\sum _{i=1}^n{S}_i\left(1-{\sigma }_{z,i}^{\left(k\right)}\right)+{∆}^2{⟨{\mathit{S}}_{\mathit{tot}}^2⟩}^{\left(k\right)},$$E101

$${∆}^2{⟨{\mathit{S}}_{\mathit{tot}}^2⟩}^{\left(k\right)}=\sum _{i=1}^n{S}_i^2\left[1-{\left(P_{\mathit{iES}}^{\left(k\right)}\right)}^2\right]-\frac{1}{2}\sum _{i<j}{n}_{\mathit{iES}}^{\left(k\right)}{n}_{\mathit{jES}}^{\left(k\right)}{S}_{\mathit{iES},\mathit{jES}}^{\left(k\right)}=0,$$E102

$$S_{\mathit{iES},\mathit{jES}}^{\left(k\right)}=\frac{4}{n\left(n-1\right)}\left\{\sum _{i\prime =1}^n{S}_{i\prime }^2\left[1-{\left(n_{i\prime \mathit{ES}}^{\left(k\right)}/2S_{i\prime }^{\left(k\right)}\right)}^2\right]\right\}{\left(n_{\mathit{iES}}^{\left(k\right)}{n}_{\mathit{jES}}^{\left(k\right)}\right)}^{-1};k=1,\dots ,2^{n-1},$$E103

Thus, we could derive the MGC-set of the internal energy functionals taking each different decomposition form from Eq. (64) involving the projected Heisenberg spin Hamiltonian:

$${⟨H_{\mathit{ex}}⟩}^{\left(k\right)}=-\frac{1}{2}\sum _{i<j}^n{J}_{i,j}\left(1-S_{\mathit{iES},\mathit{jES}}^{\left(k\right)}\right)n_{\mathit{iES}}^{\left(k\right)}{n}_{\mathit{jES}}^{\left(k\right)},k=1,\dots ,2^{n-1},$$E104

$$U_{\mathit{\text{UHFD}}-H_{\mathit{ex}}}^{\left(k\right)}\left[\rho \right]=T^{\left(k\right)}\left[\rho \right]+V_{\mathit{en}}^{\left(k\right)}\left[\rho \right]+J^{\left(k\right)}\left[\rho \right]-K_{\mathit{XC}-\mathit{Hex}}^{\left(k\right)}\left[\rho -\sum _{i=1}^n\left\{{\rho }_{\mathit{iES}}+{\rho }_{\mathit{\text{iLnb}}}\right\}\right],$$E105

$$U_{\mathit{\text{UHFD}}}^{\left(k\right)}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]=U_{\mathit{\text{UHFD}}-H_{\mathit{ex}}}^{\left(k\right)}\left[\rho \right]+{⟨H_{\mathit{ex}}⟩}^{\left(k\right)}.$$E106

Notably, one may expect that the H_ex-less internal energy function defined by Eq. (106) will become almost constant over 2ⁿ⁻¹ ESA states, owing to the sum of four different components with each having a weak k-dependency. If this is the case, one can utilize Eqs. (104) and (106) to determine the only unknown set of mean ES-exchange coupling constants, {J_i,j; i(j > i) = 1,…,n}, since n2ⁿ⁻¹ effective spin densities {n_iES^(k)}, 2ⁿ⁻²n(n-1) XC/CC-DOI ratios defined by Eq. (103) and (2ⁿ⁻¹-1) energy-difference equations (107) could be quantitatively calculated using UB3LYP/PBS(ε6)/lacvp** method:

$$∆U_{\mathit{\text{UHFD}}}^{\left(k\right)}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]\equiv U_{\mathit{\text{UHFD}}}^{\left(k\right)}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]-U_{\mathit{\text{UHFD}}}^{\left(1\right)}\left[{\rho }_{\alpha },{\rho }_{\beta }\right]\cong {⟨H_{\mathit{ex}}⟩}^{\left(k\right)}-{⟨H_{\mathit{ex}}⟩}^{\left(1\right)};k=2,\dots ,2n-1,$$E107

so that the transformed equations can be written in regular matric form:

$${\mathit{A}}^{\mathit{T}}\mathit{AX}={\mathit{A}}^{\mathit{T}}\mathit{B},$$E108

$$\mathit{A}=-\frac{1}{2}\left[\begin{array}{ccc}\mathrm{\Delta }\left[\left(1-S_{1\mathit{ES},2\mathit{ES}}^{\left(2\right)}\right)n_{1\mathit{ES}}^{\left(2\right)}{n}_{2\mathit{ES}}^{\left(2\right)}\right]& \cdots & \mathrm{\Delta }\left[\left(1-S_{\left(n-1\right)\mathit{ES},\mathit{nES}}^{\left(2\right)}\right)n_{\left(n-1\right)\mathit{ES}}^{\left(2\right)}{n}_{\mathit{nES}}^{\left(2\right)}\right]\\ ⋮& \ddots & ⋮\\ \mathrm{\Delta }\left[\left(1-S_{(n-19\mathit{ES},\mathit{nES}}^{\left(2^{n-1}\right)}\right)n_{\left(n-1\right)\mathit{ES}}^{\left(2^{n-1}\right)}{n}_{\mathit{nES}}^{\left(2^{n-1}\right)}\right]& \cdots & \mathrm{\Delta }\left[\left(1-S_{\left(n-1\right)\mathit{ES},\mathit{nES}}^{\left(2^{n-1}\right)}\right)n_{\left(n-1\right)\mathit{ES}}^{\left(2^{n-1}\right)}{n}_{\mathit{nES}}^{\left(2^{n-1}\right)}\right]\end{array}\right],$$E109

$$\mathit{B}=\left(B_2,\cdots ,B_{2^{n-1}}\right),B_k=\mathrm{\Delta }{U}_{\mathit{\text{UHFD}}}^{\left(k\right)}-\mathrm{\Delta }{U}_{\mathit{\text{UHFD}}-H_{\mathit{ex}}}^{\left(k\right)};k=2,\dots ,2^{n-1},$$E110

$$\mathit{X}=\left(X_1,\cdots ,X_{n\left(n-1\right)/2}\right)\equiv \left(X_{\mathit{ij}}\right);X_{\mathit{ij}}=J_{i,j},$$E111

where A^T is the transpose of A. Thus, we get a unique solution:

$$J_{i,j}={\left[{\left({\mathit{A}}^{\mathit{T}}\mathit{A}\right)}^{-1}{\mathit{A}}^{\mathit{T}}\mathit{B}\right]}_{\mathit{ij}}.$$E112

For n = 2, Eq. (112) reduces to

$$J_{1,2}\cong \frac{2{\mathrm{\Delta U}}_{\mathit{\text{UHFD}}}^{\left(2\right)}}{\left(1-S_{1\mathit{ES},2\mathit{ES}}^{\left(1\right)}\right)n_{1\mathit{ES}}^{\left(1\right)}{n}_{2\mathit{ES}}^{\left(1\right)}-\left(1-S_{1\mathit{ES},2\mathit{ES}}^{\left(2\right)}\right)n_{1\mathit{ES}}^{\left(2\right)}{n}_{2\mathit{ES}}^{\left(2\right)}}.$$E113

In Table 1, we show again the results of benchmark-test calculations of the ES-exchange coupling constants J_1,2, designated ($J_{1,2/\mathbf{a}}^{\left(m\right)}$, $J_{1,2/\mathbf{c}}^{\left(m\right)},\dots ,J_{1,2/\mathbf{k}}^{\left(m\right)}$), for 10 biomimetic binuclear Cu, Mn and Fe complexes, named (a, c, …, k), were made using 13 conventional mXC/PBS/lacvp** method (m = 4 ∼ 16) in place of the present MGC-SDFT-UHFD (≡0XC) method, which is unfortunately not yet implemented. These data sets were compared with the observed values, named $\left(J_{1,2/\mathbf{a}}^{\mathit{exp}},J_{1,2/\mathbf{c}}^{\mathit{exp}},\dots .,J_{1,2/\mathbf{k}}^{\mathit{exp}}\right)$, to show all the excellent quantitative agreements between the theoretical values $\left(J_{1,2/\mathbf{a}}^{\left(4\right)},J_{1,2/\mathbf{c}}^{\left(4\right)},\dots ,J_{1,2/\mathbf{k}}^{\left(4\right)}\right)$ and the experimental values mentioned above only by the standard B3LYP (≡4XC) method [13]. Here, we raise two possible explanations for the best performance by the B3LYP hybrid XC energy functional; (1) the best atomic structure of each TM-dimeric complex could be obtained by further geometry-optimization near the observed XRD structure by the B3LYP/PBS(ε)/lacvp** method [15, 16] with the dielectric constant ε of the solvent being chosen the best one from 5, 10, 20, and 40 [13]; (2) as an ideally-good balance between the exchange and correlation energy in the UHFD approximation is considered to be a key factor, B3LYP/PBS(ε)/lacvp** method may satisfy this condition most closely.

2.4 Two most stable isomers of the S₁⁽⁰⁾ state Mn₄CaO₅ clusters: Identified by two EPR signals

We have recently applied the UB3LYP/PBS(ε)/lacvp** method in place of the UHFD/ PBS(ε)/lacvp** method to all the water-splitting and oxygen-evolving reactions catalyzed by the Mn₄CaO₅ cluster in photosystem II (PSII). The electron-abstracting and proton-releasing reactions from the so-called oxygen-evolving complex (OEC) are considered to occur serially via five redox states, called Kok’s S_i-states (i = 0, 1, …, 4), where S₁ is the dark-stable state, and S₄ spontaneously decays to the initial S₀-state after releasing two protons and evolving dioxygen: the generalized reaction schemes are symbolically given by

$$S_1^{\left(0\right)}\stackrel{\mathit{-e}}{\to }{S}_2^{\left(+1\right)}\stackrel{\mathit{-e}}{\to }\stackrel{\mathit{-}{\mathit{H}}^{+}}{\to }{S}_3^{\left(+1\right)}\stackrel{\mathit{-e}}{\to }\stackrel{\mathit{-}2{\mathit{H}}^{+}}{\to }{S}_4^{\left(0\right)}\stackrel{\text{+2}{H}_{2}O}{\to }\stackrel{{\mathit{-O}}_2}{\to }{S}_0^{\left(0\right)}\stackrel{\mathit{-e}}{\to }\stackrel{{\mathit{-H}}^{+}}{\to }\left\{S_1^{\left(0\right)},\circlearrowleft \right\},$$E114

where the figure k in the superfix parentheses of S_i^(k) represents a formal charge of the ith OEC, −e above an arrow (→) indicate one electron transfer from OEC to P680⁽⁺⁾, an oxidized PSII reaction center intermittingly generated by every ∼10 μs light-pulse, −H⁺above an arrow (→) does a proton released into aqueous phase, and the symbols +(−) indicate to go out(in) of OEC, respectively. Among many controversial problems remained to be elucidated, we here take up the molecular structure of the S₁⁽⁰⁾-state Mn₄CaO₅ cluster, that is not yet established because the experimental data from XFEL, EPR and EXAFS spectroscopies appear to be apparently inconsistent if these are assumed to have been observed from the same S₁⁽⁰⁾-state. Although we can’nt exclude the possibility that the XFEL model [21] may reflect a photo-reduced S₀⁽⁻¹⁾ state of the S₁⁽⁰⁾-state Mn₄CaO₅ cluster, we have no reason to doubt the fact that two kinds of broad g4.8 and g12-multiline EPR signals were observed from the S₁⁽⁰⁾-state samples of cyanobacteria [18, 19] and spinach [20], respectively, which must have slightly different structures due to the different peripheral proteins between them. Indeed, we could prove that these EPR signals are attributable to two different structural isomers, named S_1A and S_1B in [23], which coexist with quasi-degenerate lowest energies in the respective S₁⁽⁰⁾-state. Two papers substantiating these ideas will be submitted for publication near future.

2.5 Superconductivity enhanced by the secondary correlation interaction in metals

It is well known that many materials become superconducting (S-phase) at lower temperatures than the critical temperature Tc where each system makes the transition from the normal metallic phase (N-phase). This phenomenon has been explained in terms of the Bardeen-Cooper-Schrieffer model [23, 24] combining the Fröhlich electron-lattice attractive interaction model [25] and the Bogoliubov Cooper-pairing model [26]. The highest T_c that had been achieved on 2015 is the sulfur hydride system at 203 K at high pressure (155 GP), identified from the observed magnetization vs. θ/k_B stepping-curve [26]. The observed H/D isotope effect on the down-shift & down-size of this curve appears to be consistent with the BCS model. Drozdov et al. raised three conditions required for such much higher-T_c than those of normal metals: (1) higher-frequency-phonon, (2) stronger electron-phonon coupling, and (3) a higher-density of Cooper pairing states [27]. At least the former two conditions could in principle be fulfilled for metallic and covalent compounds dominated by hydrogen. But notably the BCS model contains a serious deficiency that it is based on the free electron model but not on the Bloch-electron model depending on any approximation of self-consistent exchange-correlation potential, so that it is forced to take account only the screened Coulomb repulsive force between conductive electrons, as given by a Fourie transform,

$$\underset{\kappa \to 0^{+}}{lim}\frac{e^2}{r}{e}^{-\mathit{\kappa r}}=\int V_{\mathit{sc}}\left(\mathit{q}\right)e^{i\mathit{q}\mathbf{\bullet }\mathit{r}}\mathrm{d}\mathbf{q}\mathbf{;}{V}_{\mathit{sc}}\left(\mathit{q}\right)\mathbf{=}\underset{\kappa \to 0^{+}}{lim}\frac{4\pi e^2}{q^2+{\kappa }^2},$$E115

where κ is called “Thomas-Fermi wave number” and the limit of κ → 0⁺ implies a bare-Coulomb interaction. In order to treat such a many-body effect for Bloch-electrons near the Femi-level μ_F, we should adopt the grand potential Ω as the more appropriate thermodynamic free energy than the internal energy U, as given by

$$\mathrm{\Omega }=-\mathit{pV}=U-{\mu }_F\overline{N}-\theta \mathrm{\Sigma }.$$E116

Although here we assume that the decrease of entropy Σ upon the N-to-S phase transition may be relatively much smaller than the decrease of $U-{\mu }_F\overline{N}$ at least at low temperatures. This choice of Ω is also consistent with the fact that the higher pressure appears to be directly correlated with higher-T_c superconductors [27].

In contrast to the conventional idea of repulsive Coulomb force, the principal exchange-correlation energy Hamiltonian ${\widehat{H}}_{\mathit{\text{UHFD}}}^0$ of Eq. (57) is already incorporated in the present GC-UHFD-SDFT theory to define a binary set of Bloch eigenstates occupied with the Fermi-Dirac distribution f(${ϵ}_{\mathit{k}}-{\mu }_F$), designated ${}_B{}^{\sigma }\mathbf{M}$(σ=α, β), and there remains only the secondary correlation energy part of Eq. (60) to be regarded as giving rise to the attractive exchange force between opposite-spin itinerant Bloch-electrons. Hence, neglecting the entropic term, we need only to treat a small part of the grand potential of Eqs. (116) alone, that can be expanded by the perturbation theory:

$${}_B{\mathrm{\Omega }}_{\mathit{\text{UHFD}}}={}_B{\mathrm{\Omega }}_{\mathit{\text{UHFD}}}^0+{}_B{\mathrm{\Omega }}_{\mathit{\text{UHFD}}}^1+{}_B{\mathrm{\Omega }}_{\mathit{\text{UHFD}}}^2+\cdots ,$$E117

$${}_B{\widehat{\mathrm{\Omega }}}_{\mathit{\text{UHFD}}}^0\cong {{}_B\widehat{H}}_{\mathit{\text{UHFD}}}^0-{\mu }_F\sum _{\mathit{k}\sigma }^{{}_B{}^{\sigma }\mathbf{M}}{a}_{\mathit{k}\sigma }^{†}{a}_{\mathit{k}\sigma }=\sum _{\sigma =\alpha ,\beta }\sum _{\mathit{k}\sigma }^{{}_B{}^{\sigma }\mathbf{M}}\left[{ϵ}_{\mathit{k}\sigma }-{\mu }_F\right]a_{\mathit{k}\sigma }^{†}{a}_{\mathit{k}\sigma },$$E118

$${{}_B\widehat{\mathrm{\Omega }}}_{\mathit{\text{UHFD}}}^1=\sum _{\sigma =\alpha ,\beta }\sum _{\mathit{k\sigma },\mathit{k}\prime \mathit{\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}{V}_{\mathit{\text{UHFD}}}^{\mathit{k\sigma }\mathbf{,}\mathit{k}\prime \overline{\sigma }}{V}_{\mathit{k}\mathit{\sigma }\mathbf{,}\mathit{k}\prime \overline{\sigma }}{a}_{\mathit{k}\mathit{\sigma }}^{†}{a}_{\mathit{k}\prime \overline{\mathit{\sigma }}}^{†}{a}_{\mathit{k}\prime \overline{\sigma }}{a}_{\mathit{k}\mathit{\sigma }},$$E119

$$\begin{array}{c}{V}_{\mathit{\text{UHFD}}}^{\mathit{k\sigma }\mathbf{,}\mathit{k}\prime \overline{\sigma }}\equiv -e^2\iint \frac{{d\mathit{r}}_1{d\mathit{r}}_2}{2r_{12}}{\mathrm{\Phi }}_{\mathit{k}\sigma }^{\ast }\left({\mathit{r}}_{\mathbf{2}}\right){\mathrm{\Phi }}_{\mathit{k}\prime \overline{\mathit{\sigma }}}^{\ast }\left({\mathit{r}}_{\mathbf{1}}\right){\mathrm{\Phi }}_{\mathit{k}\prime \overline{\mathit{\sigma }}}\left({\mathit{r}}_2\right){\mathrm{\Phi }}_{\mathit{k}\sigma }\left({\mathit{r}}_1\right)\\ \times \left({\sigma }_{+,\mathit{k}\sigma ,}{\sigma }_{-,\mathit{k}\prime \overline{\sigma }}+{\sigma }_{-,\mathit{k}\sigma }{\sigma }_{+,\mathit{k}\prime \overline{\sigma }}\right)\end{array},$$E120

$$\approx -\underset{\kappa \to 0^{+}}{lim}\frac{4\pi e^{2}V}{{\left|\mathit{k}-{\mathit{k}}^{\prime }\right|}^2+{\kappa }^2}\left({\sigma }_{+,\mathit{k}\sigma ,}{\sigma }_{-,\mathit{k}\prime \overline{\sigma }}+{\sigma }_{-,\mathit{k}\sigma }{\sigma }_{+,\mathit{k}\prime \overline{\sigma }}\right),$$E121

where we put $\overline{\alpha }=\beta$ and $\overline{\beta }=\alpha ,$and used the plane-wave approximation for the conductive Bloch-wave functions: ${\mathrm{\Phi }}_{\mathit{k}\sigma }\left({\mathit{r}}_2\right)\approx exp[i\mathit{k}•$r₂], ${\mathrm{\Phi }}_{\mathit{k}\sigma }^{\ast }\left({\mathit{r}}_{\mathbf{2}}\right)\approx exp[-i\mathit{k}•$r₂] et c. to make the double integrations on the coordinate vectors r₁ and r₂ after transformed into R₁₂ = (r₁ + r₂)/2 and r₁₂ = r₂-r₁, together with introducing the screened-Coulomb damping factor exp.(−κr₁₂). Note that the volume of the system V appears from the integral on the center-of-mass coordinate R₁₂.

Notably, this spin-dependent first-order perturbation of Eq. (120) is off-diagonal in the binary Hilbert space, so that it can’nt contribute to the renormalized eigenstates (S-state) through any odd-number order of perturbation term. Then, the predominant contribution could arise from the second order perturbation term as given by

$$\begin{array}{c}{{}_B\mathrm{\Omega }}_{\mathit{\text{UHFD}}}^2=\sum _{\sigma =\alpha ,\beta }\sum _{{\mathit{k}}_{1\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}\sum _{\mathit{k}{\prime }_{1\overline{\sigma }}}^{{}_B{}^{\overline{\sigma }}\mathbf{M}}\sum _{{\mathit{k}}_{2\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}\sum _{\mathit{k}{\prime }_{2\overline{\sigma }}}^{{}_B{}^{\overline{\sigma }}\mathbf{M}}\frac{V_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{1}},\mathit{k}{\prime }_{\mathbf{1}}}{\left[V_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{2}},\mathit{k}{\prime }_{\mathbf{2}}}\right]}^{\ast }}{{ϵ}_{{\mathit{k}}_{1\sigma }}+{ϵ}_{\mathit{k}{\prime }_{1\overline{\sigma }}}-{ϵ}_{{\mathit{k}}_{2\sigma }}-{ϵ}_{\mathit{k}{\prime }_{\mathbf{2}\overline{\mathit{\sigma }}}}}\\ \times f\left({ϵ}_{{\mathit{k}}_{1\sigma }}-\mu \right)f\left({ϵ}_{\mathit{k}{\prime }_{1\overline{\sigma }}}-\mu \right)\left[1-f\left({ϵ}_{{\mathit{k}}_{2\sigma }}-\mu \right)\right]\left[1-f\left({ϵ}_{\mathit{k}{\prime }_{2\overline{\sigma }}}-\mu \right)\right],\end{array}$$E122

$$V_{\mathit{\text{UHFD}}}^{\mathit{k},\mathit{k}\prime }=-\underset{\kappa \to 0^{+}}{lim}\frac{4\pi e^{2}V}{{\left|\mathit{k}-{\mathit{k}}^{\prime }\right|}^2+{\kappa }^2}.$$E123

Significantly in the second-quantization representation this term may be transformed into

$$\begin{array}{c}{{}_B\widehat{\mathrm{\Omega }}}_{\mathit{\text{UHFD}}}^2=\sum _{\sigma =\alpha ,\beta }\sum _{{\mathit{k}}_{1\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}\sum _{\mathit{k}{\prime }_{1\overline{\sigma }}}^{{}_B{}^{\overline{\sigma }}\mathbf{M}}\sum _{{\mathit{k}}_{2\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}\sum _{\mathit{k}{\prime }_{2\overline{\sigma }}}^{{}_B{}^{\overline{\sigma }}\mathbf{M}}{V}_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{1}},\mathit{k}{\prime }_{\mathbf{1}}\mathbf{;}{\mathit{k}}_{\mathbf{2}},\mathit{k}{\prime }_{\mathbf{2}}}\\ \times {\widehat{a}}_{\mathit{k}{\prime }_{2\overline{\sigma }}}{\widehat{a}}_{\mathit{k}{\prime }_{2\overline{\sigma }}}^{†}{\widehat{a}}_{{\mathit{k}}_{2\sigma }}{\widehat{a}}_{{\mathit{k}}_{2\sigma }}^{†}{\widehat{a}}_{\mathit{k}{\prime }_{1\overline{\sigma }}}^{†}{\widehat{a}}_{\mathit{k}{\prime }_{1\overline{\sigma }}}{\widehat{a}}_{{\mathit{k}}_{1\sigma }}^{†}{\widehat{a}}_{{\mathit{k}}_{1\sigma }},\end{array}$$E124

$$V_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{1}},\mathit{k}{\prime }_{\mathbf{1}}\mathbf{;}{\mathit{k}}_{\mathbf{2}},\mathit{k}{\prime }_{\mathbf{2}}}=\underset{\kappa \to 0^{+}}{lim}\frac{{\left(4\pi e^{2}V\right)}^2}{\left({\left|{\mathit{k}}_1-\mathit{k}{\prime }_1\right|}^2+{\kappa }^2\right)\left({\left|{\mathit{k}}_2-\mathit{k}{\prime }_2\right|}^2+{\kappa }^2\right)\left({ϵ}_{{\mathit{k}}_1}+{ϵ}_{\mathit{k}{\prime }_1}-{ϵ}_{{\mathit{k}}_2}-{ϵ}_{\mathit{k}{\prime }_{\mathbf{2}}}\right)}<0.$$E125

The first point to notice is that the second-order perturbation Eq. (122) might be too complicated but could generate the attractive interaction between two Cooper-pair particles if it be approximated by the appropriate form (simply putting k′_iβ = −k_iβ; i = 1, 2 and multiplying twice the state number in each spherical-shell volume, 4πk_F³ω_D/μ_F, as given by N(k_F) ≈ 4πk_F³ω_D/μ_F (2π)³V = k_F³ω_D/2π²Vμ_F):

$${{}_B\widehat{\mathrm{\Omega }}}_{\mathit{\text{UHFD}}}^2\approx \sum _{\sigma =\alpha ,\beta }\sum _{{\mathit{k}}_{1\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}\sum _{{\mathit{k}}_{2\sigma }}^{{}_B{}^{\sigma }\mathbf{M}}{V}_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{1}},\mathbf{-}{\mathit{k}}_{\mathbf{1}}\mathbf{;}{\mathit{k}}_{\mathbf{2}},\mathbf{-}{\mathit{k}}_{\mathbf{2}}}{\widehat{a}}_{\mathbf{-}{\mathit{k}}_{2\overline{\sigma }}}{\widehat{a}}_{-{\mathit{k}}_{2\overline{\sigma }}}^{†}{\widehat{a}}_{{\mathit{k}}_{2\sigma }}{\widehat{a}}_{{\mathit{k}}_{2\sigma }}^{†}{\widehat{a}}_{\mathbf{-}{\mathit{k}}_{1\overline{\sigma }}}^{†}{\widehat{a}}_{\mathbf{-}{\mathit{k}}_{1\overline{\sigma }}}{\widehat{a}}_{{\mathit{k}}_{1\sigma }}^{†}{\widehat{a}}_{{\mathit{k}}_{1\sigma }},$$E126

$$V_{\mathit{\text{UHFD}}}^{{\mathit{k}}_{\mathbf{1}},\mathbf{-}{\mathit{k}}_{\mathbf{1}}\mathbf{;}{\mathit{k}}_{\mathbf{2}},\mathbf{-}{\mathit{k}}_{\mathbf{2}}}\approx \underset{\kappa \to 0^{+}}{lim}\frac{{\left(4\pi e^2\right)}^2{\left(k_F^3{\omega }_D/2{\pi }^2{\mu }_F\right)}^2}{2\left(4k_1^2+{\kappa }^2\right)\left(4k_2^2+{\kappa }^2\right)\left({ϵ}_{{\mathit{k}}_1}-{ϵ}_{{\mathit{k}}_2}\right)}<0;$$E127

which is an attractive potential under the BCS restrictions:

$$-{\omega }_D<{ϵ}_{{\mathit{k}}_1}-{\mu }_F<0<{ϵ}_{{\mathit{k}}_2}-{\mu }_F<{\omega }_D,$$E128

where ω_D is the Debye frequency, and note that the matrix element Eq. (126) does not contain V and it is proportional to ω_D²/(${ϵ}_{{\mathit{k}}_2}-{ϵ}_{{\mathit{k}}_1}$), which diverges as $\left({ϵ}_{{\mathit{k}}_2}-{ϵ}_{{\mathit{k}}_1}\right)\to 0$ and hence may not be approximated as a constant. Examination of this singularity problem must be postponed in future, because of the page limitation.

Up to the present stage, however, we find out that in the principal GC-SDFT-UHFD method the remained secondary correlation interaction between Bloch-electrons near the Fermi-surface could generate an additional attractive force to promote the Cooper-pairing superconductivity by increasing not only the concentration of Cooper-pair particles but also the energy gap at the Fermi level.

3. Conclusion

In this chapter, we have reviewed the MGC-SDFT-UHFD method proposed in [13] in order to advance beyond the conventional KS-DFT-UHF method. We need more clearly to explain why the KS-formalism must be regarded as incomplete, because it is a kind of double standard or hybrid theory based the quantum-mechanical rule in closed system and the thermodynamic rule in open system, as clearly seen from their use of two distinct variation-principal equations. This inconsistent theory results in two problematic notions, (1) “eternally-unknown correlation energy functional” including a separated part of kinetic energy, and (2) a set of mutually interacting LCAO-MO quasi-particles.

Here, we have widely proposed a thermodynamic alternative to derive the principal internal energy functional, which has been required to define the self-consistent one-body potential in the Schrödinger equation yielding the ultimate ground and excited states, further which have been required multiple grand canonical ensembles to properly describe all kinds of spin-dependent systems, like the paramagnetic properties of the water-splitting Mn₄CaO₅-cluster in photosystem II. This one-body quasi-particle world picture has been completed by our two revolutionary discoveries of the principal exchange-correlation energy functional, that is, a non-local exchange-correlation interaction, and a complete set of self-consistent LCAO-NMOs, which extensively span all the energy levels below dissociation limit (called the work function W) with the Fermi-Dirac distribution.

Significantly, we have presented in Sections 2.3 and 2.4 two experimental evidences directly supporting the quantitative and systematic aspects of the MGC-SDFT-UHFD method, and in Section 2.5 one more evidence indirectly supporting this UHFD decoupling scheme retaining the only secondary correlation energy functional, which spin-dependent interaction between Bloch-electrons can promote Cooper-pairings of Bloch-electrons near Fermi-level in superconductor, provided that their eigen states might be exactly determined by the MGC-SDFT-UHFD method under the crystalline periodic conditions. This implies that the Bloch-electrons near the Fermi surface are unstable in the normal phase and hence tend to make the phase transition to the superconducting phase. Further, this provides an additional mechanism for the high-temperature superconductivity. It is further emphasized that the MGC-SDFT-UHFD/PSB(ε)/lacvp** method can help meet the demand for an eagerly awaited, first principle, quantitative, and practical method to elucidate the enzymatic function of paramagnetic Mn₄CaO_x clusters in a series of water-splitting and oxygen-evolving reactions in PSII. Moreover, the present method have very high potential to be able to extend the application fields to the optical excited states, the van-der Waals interactions between fragments in the molecular system and the high-temperature superconductor.

Acknowledgments

The authors would like to acknowledge the infra-structure support for the Joint Research Programs in Graduate School of Science and Technology, Meiji University, Japan.