We are IntechOpen, the world's leading publisher of Open Access books Built by scientists, for scientists

6,800+

Open access books available

183,000+

International authors and editors

195M+ Downloads

156 Countries delivered to Our authors are among the

Top 1% most cited scientists

12.2%

Contributors from top 500 universities

Selection of our books indexed in the Book Citation Index in Web of Science™ Core Collection (BKCI)

### Interested in publishing with us? Contact book.department@intechopen.com

Numbers displayed above are based on latest data collected. For more information visit www.intechopen.com

## Meet the editors

Dr. Sajjad Haider is an associate professor in the Department of Chemical Engineering, King Saud University, Saudi Arabia. He received an MSc in 1999 and an MPhil in 2004 from the Institute of Chemical Sciences, University of Peshawar, KPK, Pakistan, and a Ph.D. in 2009 from the Department of Polymer Science and Engineering, Kyungpook National University, Taegu, South Korea. His research work focuses on the devel-

opment of scaffolds for tissue regeneration, biopolymer composites, polymer hydrogels, fabrication of electrospun nanofibers and metal nanoparticles, and evaluation of their potential application in biomedicine and removal of hazardous materials from aqueous media. He has published more than 165 research articles in reputed journals and conferences. He has edited nine books and authored fifteen book chapters. He also holds two US patents, which are also registered in Europe and Saudi patent offices. He has been invited several times to speak at international workshops and conferences.

Dr. Adnan Haider is an associate professor in the Department of Biological Sciences, Division of Nanomedicine, National University of Medical Sciences (NUMS) Rawalpindi, Pakistan. He holds a Ph.D. and post-doctorate in Polymer Science from South Korea. His research work focuses on the development of scaffolds for tissue regeneration, biopolymer composites, polymer hydrogels, drug delivery systems, fabrication of electrospun

nanofibers and metal nanoparticles, and evaluation of their potential application in biomedicine and removal of hazardous materials from aqueous media. Dr. Haider has 88 publications to his credit, including journal articles and book chapters. He is the editor of several books in his field of specialization. He is also the managing editor of Life and Science. He has been invited to speak at several international workshops and conferences.

Dr. Salah Ud-Din Khan is an associate professor at the Sustainable Energy Technologies (SET) Center, King Saud University (KSU), Saudi Arabia. Before joining KSU, he was an assistant professor, visiting scientist, researcher, and lecturer at various institutes and universities. His areas of research include clean energy sources, energy storage, radiological protection, and computational techniques for radiation shielding materials. Dr.

Khan teaches and supervises postgraduate students in renewable energy technologies at SET Center. He has published more than 250 research articles in reputed journals and conferences. He has also authored three books and two book chapters. He has one US patent to his credit, which is also registered in Europe and Saudi patent offices. He has been invited as a speaker several times to international workshops and conferences.

### Contents


Preface

Density functional theory (DFT) is a course that was registered as a "curriculum course" at the University of Munich in the 1990s, mainly covering topics such as time- and orbital-dependent functions. Over the years, DFT has proven to be a very useful theory in various chemical applications. This has prompted us to present to interested readers the most important breakthroughs and the current state of knowledge on the subject. Since the development of DFT as a quantum mechanical theory based on electron density, numerous efforts have been made to redefine important concepts already known in organic chemistry and make them useful for semiquantitative analysis. The book begins with the principles and theory from the ground states to the excited states of photoelectron emission. The detailed principles and theories are not limited to metal surfaces but also include nanoscale particles. Meanwhile, the DFT model (both time-dependent and time-independent) is widely used and has the advantage that it directly accounts for electron correlation in its formalism. In addition, DFT can be programmed for computer simulations to determine eigenvalues and eigenvectors. It provides enough information to determine the electronic structure and the total energy of the particles. An important discovery of Kohn-Sham density functional theory (KS-DFT), originally formulated half a century ago for Hamiltonians in isolated many-electron systems, was the driving force for materials science and led to a wide range of specialized codes for the prediction of molecular and crystalline properties. The first important improvement came from the discovery of the most important exchange-correlation energy functional, which led to the scheme provisionally referred to as the UHFD scheme. This scheme is the ultimate approximation beyond UHF and all other hybrid exchange-correlation energy functionals, suggesting that the search for a more powerful hybrid exchange-correlation energy functional may no longer be necessary to escape the chaotic state. The book also discusses the physicochemical properties of the drug ellipticine obtained by the DFT-B3LYP/6-311 G (d, p) method using the Gaussian computational package. HOMO (highest occupied molecular orbital), LUMO (lowest unoccupied molecular orbital), MEP (minimum energy paths) surface area, and chemical reactivity descriptors can be used computationally for

**Sajjad Haider, Ph.D.**

 College of Engineering, King Saud University, Riyadh, Saudi Arabia

Department of Chemical Engineering,

practical applications.

## Preface

Density functional theory (DFT) is a course that was registered as a "curriculum course" at the University of Munich in the 1990s, mainly covering topics such as time- and orbital-dependent functions. Over the years, DFT has proven to be a very useful theory in various chemical applications. This has prompted us to present to interested readers the most important breakthroughs and the current state of knowledge on the subject. Since the development of DFT as a quantum mechanical theory based on electron density, numerous efforts have been made to redefine important concepts already known in organic chemistry and make them useful for semiquantitative analysis. The book begins with the principles and theory from the ground states to the excited states of photoelectron emission. The detailed principles and theories are not limited to metal surfaces but also include nanoscale particles. Meanwhile, the DFT model (both time-dependent and time-independent) is widely used and has the advantage that it directly accounts for electron correlation in its formalism. In addition, DFT can be programmed for computer simulations to determine eigenvalues and eigenvectors. It provides enough information to determine the electronic structure and the total energy of the particles. An important discovery of Kohn-Sham density functional theory (KS-DFT), originally formulated half a century ago for Hamiltonians in isolated many-electron systems, was the driving force for materials science and led to a wide range of specialized codes for the prediction of molecular and crystalline properties. The first important improvement came from the discovery of the most important exchange-correlation energy functional, which led to the scheme provisionally referred to as the UHFD scheme. This scheme is the ultimate approximation beyond UHF and all other hybrid exchange-correlation energy functionals, suggesting that the search for a more powerful hybrid exchange-correlation energy functional may no longer be necessary to escape the chaotic state. The book also discusses the physicochemical properties of the drug ellipticine obtained by the DFT-B3LYP/6-311 G (d, p) method using the Gaussian computational package. HOMO (highest occupied molecular orbital), LUMO (lowest unoccupied molecular orbital), MEP (minimum energy paths) surface area, and chemical reactivity descriptors can be used computationally for practical applications.

> **Sajjad Haider, Ph.D.** Department of Chemical Engineering, College of Engineering, King Saud University, Riyadh, Saudi Arabia

#### **Adnan Haider**

Section 1

DFT Theory and Practice

 Department of Biological Sciences, Nanomedicine Section, PWD Campus, National University of Medical Sciences, Rawalpindi, Punjab, Pakistan

#### **Salah Ud-Din Khan, Ph.D.**

 Sustainable Energy Technologies (SET) Center, College of Engineering, King Saud University, Riyadh, Saudi Arabia Section 1

## DFT Theory and Practice

#### **Chapter 1**

## The Use of DFT-Based *ab-initio* Technique to Determine the Stability Difference in B2 Ti-PGM Compounds

*Ramogohlo Diale, Duduzile Nkomo, Bongani Ngobe and Maje Phasha*

#### **Abstract**

In this chapter, the density functional theory (DFT) based first-principles approach is used to predict the underlying lattice properties associated with the phase transformation and stability of B2 phase in titanium-platinum group metal (Ti-PGM) compounds. This *ab- initio* technique provides a good platform to accurately explore phase stability variation between the successful Ti-PGM shape memory alloys (SMAs) (Ti50M50, M = Rh, Pd, Ir, Pt) and other B2 Ti-PGM compounds that do not show any shape memory effect (SME), such as Ti50Os50 and Ti50Ru50. The B2 TiFe, TiNi and TiAu have also been considered in this chapter in order to draw similarities and differences. Amongst the predicted results, the heat of formation was calculated to determine the thermodynamic stability, whereas the total densities of states were used to evaluate the electronic stability of these compounds. Insights on the mechanical stability of the B2 crystals were derived from the calculated elastic constants. Mechanical instability was revealed in some compounds, indicative of a possible phase transition responsible for the intrinsic shape memory character. Although an attempt to correlate this mechanical instability with imaginary frequencies established from the phonon dispersion curves is made, the correlation is not yet conclusive due to some discrepancies observed in TiNi.

**Keywords:** B2 Ti-PGM, shape memory alloys, DFT, thermodynamic stability, electronic stability, mechanical stability, lattice dynamic stability

#### **1. Introduction**

An interest in PGM-containing SMAs continues to rise due to their unique functional properties desirable for use in various high-temperature applications such as aerospace, automotive, power plants and chemical industries [1]. Ti-PGM intermetallic compounds such as TiPt [2, 3] and TiPd [4] have gained attention as promising candidates for development of high-temperature shape memory alloys (HTSMAs). This is so because they exhibit a martensitic transformation (MT) from cubic

CsCl-type B2 to orthorhombic AuCd-type B19 phase at temperatures higher (above 373 K) than that of a well-known commercial TiNi. TiPt has attracted significant interest in various industries due to its high MT from B2 to B19 at approximately 1273 K [2], whilst the MT for TiPd has been observed at 823 K [5, 6]. Similarly, for TiIr, its B2 phase undergoes MT at 2023 K to monoclinic phase [7], whereas the MT for TiRh binary system occurs at 1118 K from B2 to tetragonal L10 [8–10]. Similarly, the MT of TiAu from B2 to B19 occurs at 900 K [11].

On the other hand, some of the Ti-PGM compounds that have been investigated were previously found to have a highly stable B2 phase to room temperature, with no sign of martensitic transformation occurring, thus no shape memory behavior was observed. Amongst those are TiOs [12] and TiRu [13–16], which showed similar behavior to TiFe [17, 18]. So far, it is not known, at least from experimental studies, if these compounds may undergo martensitic transformation at cryogenic temperatures. However, such useful information can be generated using DFT computational tools. Moreover, from scientific point of view, it is very important to establish factors that suppress MT.

It is clear that the incorporation of PGMs in Ti influences both the crystal and electronic structure and consequently, the stability of austenite phase and formation of martensite phases. These phases are thoroughly studied as they are key in the design and functionality of SMAs for specific applications [19]. Most researchers have used computational calculations to gain some insight into the underlying shape memory properties of TiPt [20], TiPd [21] and TiRh [10]. Computational studies usually include structural, elastic and electronic properties. First principle calculations have been proven to be a useful and reliable tool for studying ground-state properties of these compounds [22–25].

In addition to structural, elastic and electronic properties, there has been an increase in the use of lattice vibration properties to predict the dynamic stability for a particular phase by calculating phonon dispersions. For example, Haskins et al. [26] recently used phonon dispersion calculations to resolve discrepancies associated with determination of ground-state phases for TiNi [26, 27]. Also, phonon dispersion calculations have been used by several researchers in order to determine the lattice vibration properties and predict the phase transformation path for various alloy systems [10, 20, 27].

Although DFT predictions may not agree perfectly with experimental observations, the accuracy of predicted results can still be scrutinized against the available experimental data in order to validate the accuracy of DFT calculations. In predicting phase stability, there is a link between the predicted thermodynamic, electronic and mechanical and lattice dynamic stability, which in turn can be used to predict the expected phases at a particular temperature. Some researchers [26, 28, 29] often report these stabilities separately and although that is acceptable, we observed that these stabilities can be used in connection to gain clear understanding of phase transformations for the investigated Ti-PGM compounds.

Furthermore, in order to improve shape memory properties or induce MT in TiOs and TiRu, factors influencing the shape memory behavior must be understood and this can be done by studying the predicted mechanical, electronic and dynamic stabilities. Thus, the main aim of this chapter is to demonstrate the versatile use of *ab-initio* methods based on DFT to track the shape memory behavior of these wellknown Ti-PGM compounds by evaluating their ground-state stability of the hightemperature austenite phase (B2) with reference to mechanical, electronic and dynamical stability. This will assist the reader in identifying the underlying

fundamental factors that drive the existence of shape memory behavior and further give an insight towards the design of new and improvement of existing SMAs.

#### **2. Computational methodology**

#### **2.1 Density functional theory (DFT)**

Computational modeling techniques offer an alternative way of investigating the properties of materials using computers, whereby the simulator builds a model of a real system and explores its behavior. One of the computational techniques that have received immense attention over the last few decades is the first-principles calculations, also known as *ab-initio* calculations. This interest is attributed to significant usefulness and insightfulness of its calculated data in the materials design.

First-principles methods are based on DFT formalism in which properties of materials, that is, the values of the fundamental constants and the atomic numbers of the atoms present can be calculated using the Schrödinger equation. Due to its improved accuracy achieved over the past years in predicting properties of real solids, an *ab-initio* approach is adopted in the current work to predict the ground-state and structural properties of several B2 intermetallic compounds at equiatomic compositions. Computing total energies of any system is a necessary starting point for firstprinciple calculations.

DFT is a quantum-mechanical method used for calculating ground-state properties of condensed matter systems without dealing directly with many electron states [30]. It was first formulated by Hohenberg and Kohn in 1964 [31] and then secondly developed by Kohn and Sham in 1965 [32]. DFT has helped in the development of independent-particle methods that take into account the particle's correlations and interactions. Hohenberg-Kohn demonstrated the first theorem that the ground-state properties of a many-electron system are uniquely determined by an electron density that depends on only three spatial coordinates [30].

$$E = E\left[n\left(\overrightarrow{r'}\right)\right] \tag{1}$$

Where *E* is the total energy and *n* is the density. Within the Kohn-Sham scheme [30], consideration of an interacting electron gas moving in an external potential *ve*ð Þ*r* , as a variational principle leads to the effective single-electron Schrödinger equations,

$$\left\{-\nabla^2 + v\left([n]; \vec{r}\right)\right\}\mu\_j\left(\vec{r}\right) = \mathsf{E}\_j\mu\_j\left(\vec{r}\right) \tag{2}$$

Kohn-Sham electrons in an effective potential, *veff* for a system of non-interacting, is solved as follows:

$$\mathbf{V}\_{\text{eff}}\left(\overrightarrow{r}\right) = \mathbf{v}\left(\overrightarrow{r}\right) + \int \frac{\mathbf{n}\left(\overrightarrow{r}'\right)}{|\overrightarrow{r} - \overrightarrow{r}'|} d\overrightarrow{r}' + \frac{\delta E\_{\text{xc}}\left[\mathbf{n}\left(\overrightarrow{r}\right)\right]}{\delta n\left(\overrightarrow{r}\right)}\tag{3}$$

Where *v r*!� � is the external potential and *Exc n r*! h i � � is the exchange-correlation density functional [30].

#### **2.2 Approximations to exchange-correlation functional**

The two main types of exchange-correlation functionals used in DFT are the local density approximation (LDA) [33] and the generalized gradient approximation (GGA) [34], which have been discussed in the sub-sections below.

#### *2.2.1 Local density approximation*

The local density approximation (LDA) is an approximation in which the exchange-correlation (XC) energy functional in density functional theory (DFT) depends upon the value of the electronic density at each point in space. It was first discovered by Kohn and Sham, which is expressed in the Eq. (4):

$$E\_{\rm xc}^{LDA}[n(r)] = \int dr n(r) \varepsilon\_{\rm xc}(n(r)) \tag{4}$$

Where *εXC*ð Þ *n* is the exchange-correlation energy per electron in a uniform electron gas of density *n* [33]. The uniform electron gas represents a group of systems of interacting electrons with an arbitrary spatially constant density *n,* which acts as a parameter. The local density approximation quantity is known for the limit of high density and can be calculated accurately at densities of interest by the use of Monte Carlo methods. LDA has been proven to give accurate results for many atomic, molecular and crystalline interacting electron systems, even though in these systems, the density of electrons is not slowly varied.

#### *2.2.2 Generalized gradient approximation*

The GGA is known to be a semi-local approximation, which means that there is no use of local density *n r*ð Þ value but its gradient ∇*n r*ð Þ. Perdew and Wang developed generalized gradient approximation (GGA), which is based on a real-space cut-off of the spurious long-range components for the second-order gradient expansion for the exchange-correlation hole [34]. GGA improves total energies, atomization energies, energy barriers and also the difference in structural energies. GGA takes the form:

$$E\_{\rm xr}^{GGA}\left[n\left(\overrightarrow{r}\right)\right] = \int e\_{\rm xr}^{GGA}\left(n\left(\overrightarrow{r}\right), |\nabla n\left(\overrightarrow{r}\right)|\right) n\left(\overrightarrow{r}\right) d\overrightarrow{r} \tag{5}$$

There are several GGA-based functionals, that is, the PBE [35], PBEsol [36], RPBE [37], BLYP [38] and AM05 [39]. Other known GGA-based functionals are meta-GGA [40], hyper-GGA and generalized random phase approximation.

In this chapter, the GGA-PBE [35] functional was used to optimize the Ti50M50 (M = PGMs, Ni, Fe, Au) systems as it provides accurate parameters for these materials.

#### **2.3 Computational code and implementation**

#### *2.3.1 CASTEP code*

In this book chapter, the plane-wave Cambridge Serial Total Energy Package (CASTEP) [41, 42] code was used to investigate the properties of the B2 structures.

#### *The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

CASTEP is a module embedded within the Materials Studio software package. It is a first-principle quantum mechanical code based on DFT formalism, used for performing electronic structure calculations. It can be used to simulate a wide range of materials, including crystalline solids, surfaces, molecules, liquids and amorphous materials; the properties of any material that can be thought of as an assembly of nuclei and electrons can be calculated with the only limitation being the finite speed and memory of the computers being used. This approach to simulation is extremely ambitious given that the aim is to use no experimental data but to rely purely on quantum mechanics.

Aiming to calculate any physical property of the system from first principles, the basic quantity is the total energy from which many other quantities are derived. For example, the derivative of total energy concerning atomic positions results in the forces and the derivative concerning cell parameters gives stresses. To do this, the total-energy code, on CASTEP code, performs a variational solution to the Kohn-Sham equations by using a density mixing scheme to minimize the total energy and also conjugate gradients to relax the ions under the influence of the Hellmann-Feynman forces. CASTEP uses fast Fourier transforms (FFTs) to provide an efficient way of transforming various entities (wave functions and potentials) from real to reciprocal space and back, as well as to reduce the computational cost and memory requirement for operating with the Hamiltonian on the electronic wave functions, a plane-wave basis for the expansion of the wave functions. These are then used to perform full geometry optimizations [41, 42].

A summary of the methodology for electronic structure calculations as implemented in CASTEP is as follows: a set of one-electron Schrödinger (Kohn-Sham) equations are solved using the plane-wave pseudopotential approach. The wave functions are expanded in a plane wave basis set defined by the use of periodic boundary conditions and Bloch's Theorem. The electron-ion potential is described employing *ab initio* pseudopotentials within both norm-conserving and ultrasoft formulations [41, 42].

Direct energy minimisation schemes are used to obtain self-consistent electronic wave functions and corresponding charge density. In particular, the conjugate gradient and density mixing schemes are implemented. Also, the robust electron ensemble DFT approach can be used for systems with partial occupancies, in particular, metals [41, 42].

#### *2.3.2 Implementation*

**Figure 1** illustrates the equiatomic B2 crystal geometry used to carry out all the calculations reported in this chapter.

The resulting geometry-optimized crystal structure was used to carry out all the calculations, including structural, thermodynamic and elastic properties, of all considered compounds. Only valence electrons were considered through the use of ultrasoft pseudopotentials [41, 42]. All of our calculations were performed with pseudo-potentials in generalized gradient approximation (GGA) [32] refined by Perdew, Burke and Ernzerhof (PBE) [35]. Before any calculation could be performed, a convergence test was also conducted in this code to determine the suitable cut-off energy and k-point mesh parameter for systems. A plane wave cut-off energy of 500 eV was found to be sufficient enough to converge the total energy of the systems. The Brillouin zone (BZ) sampling was performed using the k-point mesh of 13x13x13 according to the Monkhorst–Pack method [43]. A full geometry optimization was performed to determine the ground-state parameters for the binary systems.

**Figure 1.**

*Schematic representation of B2 crystal structure of Ti50M50 (M = Ru, Rh, Pd, Os, Ir, Pt, Ni, Fe, Au) intermetallic compounds.*

To obtain the stable structure of Ti50M50 with minimum total energy, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization scheme [44] was performed in geometry optimization. The maximum ionic Hellmann-Feynman force was given to be 0.03 eV/Å. Furthermore, the elastic constants of a solid were calculated by an efficient strain–stress method through a linear least-square fit of first-principles calculation results. The maximum stress set is below 0.05 GPa. The phonons dispersion curves were also calculated.

#### **2.4 Theoretical background on calculated properties**

#### *2.4.1 Thermodynamic stability*

The heat of formation (Δ*Hf* ) is the enthalpy change when one mole of a compound is formed from the constituent elements in their stable states and is essential in determining the structural stabilities of the different crystal structures. The heat of formation is estimated by the following expression:

$$
\Delta H\_f = E\_C - \sum\_i \mathbf{x}\_i E\_i \tag{6}
$$

where *EC* is the calculated total energy of the compound and *Ei* is the calculated total energy of the constituent element in the compound. For a structure to be stable, the heat of formation must have the lowest negative value (Δ*Hf* < 0). The heat of formation is used to determine the stability trend of the B2 systems.

#### *2.4.2 Electronic stability*

The electronic stability is determined by the density of states (DOSs), which refer to the occupancy and density of the electronic states in a crystalline solid. It is

*The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

described by a function, g (*E*), as the number of electrons per unit volume and energy with electron energies near Fermi-level *EF*. In the case of states with DOS being zero implies no state has been occupied (empty orbital). In general, a DOS is an average of all available spaces multiplied by the number of states occupied by the system. The local density of states (LDOSs) is a measure of variation due to distortion of the original system. LDOS can locally be non-zero if the DOS of an undisturbed system is zero due to the presence of local potential. In that case, the DOSs are the total number of states that are available in the system within the plane-wave framework of DFT. It is possible to calculate each orbital's contribution (partial DOS) to determine which orbitals are occupied or involved in bonding. The electronic behavior of a material is determined by the location of *EF* within the DOS. Alloys'stability can be predicted using the DOS.

In the case of partial density of states (PDOSs), the states are attributed to the basic functions and then to the atoms constituting the unit cell. DOS is then calculated as the sum of atomic contributions. The DOS is calculated by using the following expression:

$$m(\boldsymbol{\varepsilon}) = 2 \sum\_{n,k} \delta(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}\_n^k) = \frac{2}{V\_{BZ}} \sum\_n \int \delta(\boldsymbol{\varepsilon} - \boldsymbol{\varepsilon}\_n^k) d\boldsymbol{k} \tag{7}$$

where *δ* is the Dirac delta function and the *k* is integral extends over the BZ. The number of the electrons in the unit cell is given by:

$$\int\_{-\infty}^{\underline{\ell}} \pi(\varepsilon)d\varepsilon.$$

#### *2.4.3 Mechanical stability*

There are various criteria established to deduce the mechanical stability of crystals for different lattice crystals. Accuracy in determining the elasticity of a compound is vital in understanding its mechanical stability and elastic properties. The elastic constants depend on the type of lattice *i.e.* for the cubic, there are three (c11, c12, c44) independent elastic constants [20, 45]. For example, applying two types of strains ð*ε*1*ε*4Þ to the cubic system gives stresses relating to three elastic coefficients, this is a useful method for obtaining elastic constants. The mechanical stability condition for the cubic system as outlined in Ref. [45] is given as follows:

$$c\_{44} > 0; \ c\_{11} > c\_{12} \text{ and } c\_{11} + 2c\_{12} > 0\tag{8}$$

According to Born-Huang's lattice dynamical theory [46, 47], the stability criterion for the elastic constants must be completely satisfied for the structure to be stable. The positive *C*<sup>0</sup> = (1/2(c11–c12) > 0) indicates the mechanical stability of the crystal, otherwise, it is unstable.

#### *2.4.4 Lattice dynamic stability*

A phonon dispersion curve along a high symmetry direction is calculated by using interplanar force constants [45], as every plane perpendicular to this direction is

displaced within an elongated supercell. Generally, lattice dynamics are analyzed with the *ab-initio* evaluation of forces on all atoms resulting from finite displacements of few atoms within otherwise perfect crystals. It is usually necessary to construct supercells of the appropriate size to ensure that interactions of the perturbation with all its translational symmetry equivalent copies are small. The techniques for selecting suitable supercells and atomic displacements, assembling force constant matrices from the calculated forces and calculating phonon dispersion relations *via* Fourier transform are well documented.

Using one of the 230 crystallographic space groups, phonon constructs a crystal structure, calculates the Hellmann-Feynman force constant, builds the dynamical matrix, diagonalizes it and calculates phonon dispersion relations [45].

In phonon dispersion calculations, polarization vectors and irreducible representations (Gamma points) of phonon modes are found, and the total and partial phonon densities are calculated. It plots the internal energy, free energy, entropy, heat capacity and tensor of mean square displacements (Debey-Waller factor). Phonons calculate the dynamical structure factor of coherent inelastic neutron scattering and incoherent doubly differential scattering in single crystal and polycrystalline systems [45].

The properties of phonons can be determined using a harmonic approximation with one fundamental quantity, the force constants matrix [45]:

$$D\_{\mu\nu}(R - R') = \frac{\delta^2 E}{\partial u\_{\mu}(R)\partial u\_{\nu}}|\_{u=0} \tag{9}$$

Where *u* represents the displacement of a given atom and *E* is the total energy in the harmonic approximation. This matrix of force constants can also be represented in reciprocal space and is known as a dynamic matrix:

$$D\_{\mu\nu}(q) = \frac{1}{N\_R} \sum\_{R} D\_{\mu\nu}(R) \exp(-iqR) \tag{10}$$

Classical equations of motion can be written as eigenvalue problems with each atomic displacement in the form of plane waves:

$$\mu(R,t) = \varepsilon e^{(qR - o(q)t)} \tag{11}$$

Where ε is the 3N-dimensional eigenvector of the eigenvalue problem is:

$$M o^2(q) \varepsilon = D(q) \varepsilon \tag{12}$$

The ω on the wave vector is well-known as the phonon dispersion. A guide to the basic theory of phonons has been described in detail by Born and Huang (1954) and Ashcroft and Mermin (1976) [48, 49]. In this chapter, we use the CASTEP code [41, 42] to calculate phonon dispersion curves and their density of states (PHDOS). It is well documented that compounds that radiate real (only positive) vibrational modes along high symmetry directions in the Brillouin zone are considered to be vibrationally stable with no possibility to undergo a martensitic phase transition at lower temperatures [10, 13]. On the other hand, compounds that radiate both positive and negative vibrational modes turn out to be vibrationally unstable with high chances to undergo a martensitic phase transition at lower temperatures, a primary feature of alloys with shape memory effect [20, 23, 50].

#### **3. Results and discussion**

#### **3.1 Structural and thermodynamic stability**

The stability of the investigated B2 compounds is first discussed on the basis of heat of formation. **Table 1** presents the ground-state lattice parameters and the heat of formation of the investigated B2 phases. The reported lattice parameters were obtained from the geometrically relaxed structures, whereby their respective volumes and unit cells were allowed to change to obtain their ground state.

The lattice parameters were found to be in good agreement with those reported previously by other authors [10, 12, 13, 19, 22, 54, 55]. The lattice parameters of the investigated Ti50M50 compounds increase with the position of the M atom along the groups or the periods of the periodic table of elements. This is in line with their respective densities and volume changes observed, as presented in **Table 1**.

Furthermore, **Table 1** also presents the thermodynamic stability of the investigated B2 compounds that were determined by calculating their respective heat of formation using Eq. (6). Heat of formation provide primary insight into the existence of the phase. All the investigated B2 phases reported in this research work were found to be negative (*ΔHF* < 0), an indication that they are all thermodynamically stable. Amongst the investigated B2 compounds, the heat of formation value of Ti50Ir50 was found to be the most thermodynamically stable indicated by the lowest value of 0.913 eV/atom. The reported heat of formation results were found to be in accordance with other data in literature [51, 52, 55–58].

#### **3.2 Electronic stability**

**Figure 2** shows the total density of states (tDOS) of all the investigated B2 compounds reported in this chapter. It is noted that the DOS values *g* (*E*) of the investigated compounds were found to be non-zero across the Fermi level (*EF*) indicating that all the investigated B2 compounds were mainly characterized by metallic bonds.


**Table 1.**

*Lattice parameters and heat of formation of the investigated binary B2 compounds.*

**Figure 2.** *The total density of states (tDOS) of the nine investigated austenite compounds.*

Generally, the density of states'spectra of CsCl-type B2 compounds consists of two peaks [25] that are separated by a pseudogap (deep valley) at the Fermi level (*EF*). With the lower energy peak representing the anti-bonding region and the higher energy peak representing the bonding region. Therefore, the position of *EF* at the pseudogap provides insight into the phase stability of a compound at the ground state. Such that, if *EF* is found to fall at the centre of the pseudogap, that signifies phase stability and if *EF* falls at the peak or shoulder of the bonding region. **Figure 2** shows that Ti50Fe50, Ti50Ru50 and Ti50Os50 will retain their stable high-temperature austenite phase as their ground-state as their *EF* cuts at the centre of the deep valley, while the other investigated B2 compounds (Ti50M50, M = Ni, Rh, Pd, Ir, Au and Pt) will certainly undergo a phase transition from the high-temperature B2 to unstable phases at 0 K because their pseudogap shifted towards the anti-bonding region and enabled phase transition.

The aforesaid provides insight information about the phase transition of the investigated compounds from high-temperature austenite (B2) to low-temperature martensite phase, a gauge for shape memory effect observed on alloys with shape memory properties.

#### **3.3 Dynamic phase stability**

The phonon-dispersion curves are used to gain insight into the underlying lattice vibrations that may influence the ability of the crystal to transform to martensite phase on cooling. The information on phonons is very useful for accounting variety of properties and behaviors of crystalline materials, such as thermal properties, phase transition, and superconductivity [59].

As detailed in Section 2.4.4, B2 compounds that show only positive frequencies remain stable with no prospect of undergoing martensitic phase transition, while those that show both positive and negative frequencies are prone to become unstable and undergo martensitic phase transition at lower temperatures.

*The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

**Figure 3** represents the sets of phonon-dispersion curves of the investigated compounds. It can be shown that Ti50Fe50, Ti50Ru50 and Ti50Os50 consist of only the real positive vibrational modes, while the rest of the investigated compounds were found to consist of both positive and negative vibrational modes. This observed behavior agrees well with the DOS results reported in the previous section as well as the results reported by other authors [20, 23, 50].

**Figure 3.** *Phonon-dispersion curves of the investigated B2 phases plotted along selected Brillouin zone directions.*

Furthermore, **Figure 3** gives the calculated B2 phonons density of states (PHDOS). It is evident that the vibrational stable B2 compounds (Ti50M50, M = Fe, Ru and Os) consist of a valence gap between the positive and negative frequency modes. Such observed behavior renders the aforementioned compounds to be brittle-metal-compounds, a kind of semi-conductive material properties. While those austenite phases that were found to be vibrationally unstable (Ti50M50, M = Ni, Rh, Pd, Ir, Pt and Au) show a clear interaction of valence and conduction electrons with a non-zero DOS.

#### **3.4 Mechanical stability**

Elastic constants are the primary output parameters of most DFT codes, and they determine the response of the crystal to external forces. They are crucial in predicting the mechanical properties of the crystal structure and their subsequent modulus of elasticity.

**Table 2** presents the calculated elastic constants of the investigated alloys reported in this chapter. Furthermore, the trend of the elastic constants plotted against different B2 compounds is shown in **Figure 4**.

Based on the mechanical stability criterion for cubic crystals as outlined in Section 2.4.3, as well as very high value of C<sup>0</sup> , Ti50Fe50, Ti50Ru50 and Ti50Os50 completely


**Table 2.**

*The calculated elastic constants (Cij) of the investigated Ti50M50 compounds.*

**Figure 4.** *The calculated elastic constants (GPa) of the investigated compounds.*

*The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

satisfy the criterion, rendering the austenite phase to be mechanically stable in these compounds, thus unlikely to undergo any phase transition at lower temperatures.

Although a similar stability criterion was observed (*C*<sup>0</sup> *> 0*) for Ti50Ni50, which is a well-known SMA [14, 60], it should be noted that its *C44* was found to be larger than *C*0 , resulting in anisotropic factor *A* greater than 1. This is another important elastic property factor worth attention in design of SMAs and can be assessed on the relationship between the tetragonal shear modulus (*C*<sup>0</sup> ) and its corresponding monoclinic shear constant (*C44*). It is reported that the mechanical instability of a hightemperature phase decreases when nearing the martensitic transition temperature during heating [61], this is the case where *C*<sup>0</sup> becomes close to zero. Furthermore, if *C*<sup>0</sup> < 0, the corresponding B2 phase is unstable at 0 K, signaling a great potential for this high-temperature phase to undergo a martensitic phase transition on cooling, yielding SME character. Since this MT occurs at much higher temperatures, such B2 alloys are of interest for high-temperature structural applications. **Table 2** shows that Ti50Os50 and Ti50Ir50 are the most and least mechanically stable amongst the investigated B2 compounds, and the same results are graphically presented as shown in **Figure 4**.

#### **4. Conclusions**

The thermodynamic, electronic, lattice dynamic and mechanical stability of the investigated B2 compounds have been calculated using a first-principles approach and are reported in this chapter. The calculated heat of formation, elastic constants, densities of states and vibrational properties have provided a much-needed clarity in determining the differences in stabilities of various Ti-PGM B2 compounds. This was made possible by choosing the correct model, resulting in reliable data that compares well with results reported in literature.

From all the calculated properties, it was shown that the higher-temperature austenite phase of Ti50Fe50, Ti50Ru50 and Ti50Os50 remains stable with no prospect of martensitic phase transition that is associated with materials with shape memory effect.

This was evident by the electronic and dynamic stability of the investigated compounds. Furthermore, their corresponding tDOS spectra were found to coincide with the pseudogap at the *EF*, rendering electronic stability with no phase transition. This was further substantiated by phonon curves, which possessed only positive vibrational frequency indicating that they are not likely to undergo a phase transition. The predicted mechanical stability of the Ti50Fe50, Ti50Ru50 and Ti50Os50 B2 compounds was found to be very high, strongly signaling absence of any possible phase transformation.

On the other hand, using the electronic and elastic properties, this work has shown that other B2 compounds (Ti50M50, M = Ni, Rh, Pd, Ir, Pt and Au) considered in this study are unstable at 0 K, thus predicting the possibility to undergo phase transition to martensite phase on cooling, resulting in shape memory effect character if the transformation is diffusionless.

#### **Acknowledgements**

This book chapter is published with the permission of MINTEK. The authors would like to thank the Advanced Metals Initiative (AMI) of the Department of Science and Innovation (DSI) for financial support. The gratitude is also extended to the Centre for High-Performance Computing (CHPC) in Cape Town for allowing us to carry out the calculations using their remote computing resources.

### **Author details**

Ramogohlo Diale<sup>1</sup> , Duduzile Nkomo1,2\*, Bongani Ngobe1,3 and Maje Phasha<sup>1</sup>

1 Advanced Materials Division, Mintek, Johannesburg, South Africa

2 Department of Materials Science and Metallurgical Engineering, University of Pretoria, Pretoria, Hatfield, South Africa

3 School of Physics, University of the Witwatersrand, Johannesburg, South Africa

\*Address all correspondence to: dudunk@mintek.co.za

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

#### **References**

[1] Jani JM, Leary M, Subi A, Gibson MA. A review of shape memory alloy research, applications and opportunities. Materials and Design. 2014;**56**:1078-1113

[2] Biggs T, Cortie MB, Witcomb MJ, Cornish LA. Martensitic transformations, microstructure, and mechanical workability of TiPt. Metallurgical and Materials Transactions A. 2001;**32**:1881-1886

[3] Yamabe-Mitarai Y, Hara T, Miura S, Hosoda H. Shape memory effect and pseudoelasticity of TiPt. Intermetallics. 2010;**18**:2275-2280

[4] Otsuka K, Oda K, Ueno Y, Piao M, Ueki T, Horikawa H. The shape memory effect in a Ti50Pd50 Alloy. Scripta Metallurgica et Materialia. 1993;**29**: 1355-1358

[5] Solomon VC, Nishida M. Martensitic transformation in Ti-Rich Ti–Pd shape memory alloys. Materials Transactions. 2002;**43**:908-915

[6] Hisada S, Matsuda M, Yamabe-Mitarai Y. Shape change and crystal orientation of B19 martensite in equiatomic TiPd alloy by isobaric test. Meta. 2020;**10**:3754

[7] Yamabe-Mitarai Y, Hara T, Phasha MJ, Ngoepe PE, Chikwanda HK. Phase transformation and crystal structure of IrTi. Intermetallics. 2012;**31**: 26-33

[8] Yi SS, Chen BH, Franzen HF. Phase transitions in RhTi. Journal of the Less-Common Metals. 1988;**143**:243

[9] Szajek A. The electronic and superconducting properties of ordered Ti-Rh alloys. Journal of Physics. Condensed Matter. 1991;**3**:1089

[10] Si L, Jiang Z-Y, Zhou B, Chen W-Z. First-principles study of martensitic phase transformation of TiRh alloy. Physica B. 2012;**407**:347-351

[11] Donkersloot HC, Van Vucht JHN. Martensitic transformations in Gold– Titanium, Palladium–Titanium, and Platinum–Titanium alloys near the equiatomic composition. Journal of the Less-Common Metals. 1970;**20**: 83-91

[12] Murray JL. Ti Os-Ti (Osmium-Titanium) system. Bulletin of Alloy Phase Diagrams. 1982;**3**:212-215

[13] Kong Z, Duan Y, Peng M, Qu D, Bao L. Theoretical predictions of thermodynamic and electronic properties of TiM (M=Fe, Ru and Os). Physica B: Condensed Matter. 2019;**573**: 13-21

[14] Xing W, Chen XQ, Li D, Li Y, Fu CL, Meschel SV, et al. First-principles studies of structural stabilities and enthalpies of formation of refractory intermetallics: TM and TM3 (T=Ti, Zr, Hf; M=Ru, Rh, Pd, Os, Ir, Pt). Intermetallics. 2012;**28**:16-24

[15] Murray JL. The Ru-Ti (Ruthenium-Titanium) system. Bulletin of Alloy Phase Diagrams. 1982;**3**:216-221

[16] Jain E, Pagare G, Chouhan SS, Sanyal PS. Electronic structure, phase stability and elastic properties of ruthenium based four intermetallic compounds: Ab-initio study. Intermetallics. 2014;**59**:79-84

[17] Cacciamani G, De Keyzer J, Ferro R, Klotz UE, Lacaze J, Wollants P. Critical evaluation of the Fe–Ni, Fe–Ti and Fe–Ni–Ti alloy systems. Intermetallics. 2006;**14**:1312-1325

[18] Murray JL. The Fe-Ti (Iron-Titanium) system. Bulletin of Alloy Phase Diagrams. 1981;**2**:320-334

[19] Wadood A, Yamabe-Mitarai Y. Recent research and developments related to near-equiatomic titaniumplatinum alloys for high-temperature applications. Platinum Metals Review. 2014;**58**:61-67

[20] Mahlangu R, Phasha MJ, Chauke HR, Ngoepe PE. Structural, elastic and electronic properties of equiatomic PtTi as potential hightemperature shape memory alloy. Intermetallics. 2013;**33**:27-32

[21] Chen XQ, Fu CL, Morris JR. The electronic, elastic, and structural properties of Ti–Pd intermetallics and associated hydrides from first principles calculations. Intermetallics. 2010;**18**: 998-1006

[22] Bao W, Liu D, Duan Y, Peng M. First-principles predictions of anisotropies in elasticity and sound velocities of CsCl-type refractory intermetallics: TiTM, ZrTM and HfTM (TM = Fe, Ru, Os). Philosophical Magazine. 2019;**99**:2681-2702

[23] Baloyi ME, Modiba R, Ngoepe PE, Chauke HR. Structural, elastic and electronic properties of titanium-based shape memory alloys. In: The Proceedings of South African Institute of Physics 2018. South Africa: SAIP; 2018. pp. 8-13. ISBN: 978-0-620-85406-1

[24] Hart GLW, Curtarolo S, Massalski TB, Levy O. Comprehensive search for new phases and compounds in binary alloy systems based on platinum-group metals, using a computational first-principles approach. Physical Review X. 2013;**3**: 041035

[25] Zhang JM, Guo GY. Electronic structure and phase stability of three series of B2 Ti-transition-metal compounds. Journal of Physics: Condensed Matter. 1995;**7**:6001-6017

[26] Haskins JB, Lawson JW. Ab initio simulations of phase stability and martensitic transitions in NiTi. Physical Review B. 2016;**94**(21):1-10

[27] Haskins JB, Hessam M, Honrao SJ, Sandoval LA, Lawson JW. Lowtemperature mechanical instabilities govern high-temperature thermodynamics in the austenite phase of shape memory alloy constituents: Ab Initio simulations of NiTi, NiZr, NiHf, PdTi, and PtTi. Acta Materialia. 2021; **212**:1-13

[28] Vishnu KG, Alejandro S. Phase stability and transformations in NiTi from density functional theory calculations. Acta Materialia. 2010;**58**(3): 745-752

[29] Abdoshahi N, Dehghani M, Hatzenbich L, Spoerk-Erdely P, Ruban AV, Musi M, et al. Structural stability and mechanical properties of TiAl+Mo alloys: A comprehensive ab initio study. Acta Materialia. 2021;**221**: 1-14

[30] Mattson AE, Schultz PA, Desjarlais MP, Mattsson TR, Leung K. Designing meaningful density functional theory calculations in materials science-a primer. Modelling and Simulation in Materials Science and Engineering. 2005; **13**:1-32

[31] Hohenberg P, Kohn W. Introduction to density functional theory. Physical Review B. 1964;**136**:864

[32] Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Physical Review A. 1965;**140**:1133

*The Use of DFT-Based* ab-initio *Technique to Determine the Stability… DOI: http://dx.doi.org/10.5772/intechopen.112933*

[33] Hedin L, Lundqvist BI. Explicit local exchange-correlation potentials. Journal of Physics C. 1971;**4**:2064-2082

[34] Perdew JP, Wang Y. Accurate and simple analytic representation of the electron-gas correlation energy. Physical Review B. 1992;**45**:13244-13249

[35] Perdew JP, Burke K, Ernzerhof M. Generalized gradient approximation made simple. Physical Review Letters. 1996;**77**:3865-3868

[36] Perdew JP, Ruzsinszky A, Csonka GI, Vy-drov OA, Scuseria GE, Constantin LA, et al. Restoring the density-gradient expansion forexchange in solids and surfaces. Physical Review Letters. 2008;**100**:136406

[37] Hammer B, Hansen LB, Nørskov JK. Improved adsorption energetics within density-functionaltheory using revised Perdew-Burke-Ernzerhof functionals. Physical Review B. 1999;**59**:7413

[38] Gill PMW, Johnson BG, Pople JA, Frisch MJ. The performance of the Becke-Lee-Yang-Parr (B-LYP) density functional theory with various ba-sis sets. Chemical Physics Letters. 1992;**197**: 499

[39] Armiento R, Mattsson AE. Functional designed toinclude surface effects in self-consistent density functional theory. Physical Review B. 2005;**72**:085108

[40] Tao J, Perdew JP, Staroverov VN, Scuseria GE. Climbing the density functional ladder: nonempirical metageneralized gradient approximation designed for molecules and solids. Physical Review Letters. 2003;**91**:146401

[41] Segall MD, Philip JDL, Probert MJ, Pickard CJ, Hasnip PJ, Clark SJ, et al. First-principles simulation: Ideas,

illustrations and the CASTEP code. Journal of Physics D: Condensed Matter. 2002;**14**:2717-2744

[42] Clark SJ, Segall MD, Pickard CJ, Hasnip PJ, Probert MI, Refson K, et al. First principles methods using CASTEP. Zeitschrift für Kristallographie. 2005; **220**:567-570

[43] Monkhorst HJ, Pack JD. Special points for Brillouin-zone integrations. Physical Review B. 1996;**13**:5188-1592

[44] Fischer T, Almlof J. General methods for geometry and wave function optimization. Journal of Physical Chemistry. 1992;**96**:9768-9774

[45] Mehl MJ, Klein BM. First-principles calculation of elastic properties. In: Fleischer JW, editor. Intermetallic Compounds – Principles and Practice. Vol. 1. London: Principles, John Wiley and Sons; 1994. pp. 195-210

[46] Gomisa O, Manjónb FJ, Rodríguez-Hernándezc P, Muñozc A. Elastic and thermodynamic properties of α-Bi2O3 at high pressures: Study of mechanical and dynamical stability. Journal of Physics and Chemistry of Solids. 2019;**124**: 111-120

[47] Wagner MF-X, Windl W. Lattice stability, elastic constants and macroscopic moduliof NiTi martensites from first principles. Acta Materialia. 2008;**56**:6232-6245

[48] Born M, Huang K. Dynamical Theory of Crystal Lattices. Oxford: Oxford University Press; 1954

[49] Ashcroft NW, Mermin ND. Solid State Physics. Philadelphia: Saunders College; 1976

[50] Diale RG, Modiba R, Ngoepe PE, Chauke HR. The effect of Ru on

Ti50Pd50 high temperature shape memory alloy: A first-principles study. MSR Advances. 2019;**4**:2419-2429

[51] Semenova EL. Alloys with the shape memory effect in systems of d-transition metals containing metals of the platinum group. Powder Metallurgy and Metal Ceramics. 1997;**36**:394-404

[52] Gachon JC, Selhaoui N, Aba B, Hertz J. Comparison between measured and predicted enthalpies of formation. Journal of Phase Equilibria. 1992;**13**: 506-511

[53] Mizuno M, Akari H, Shirai Y. Compositional dependence of structures of NiTi martensite from first principles. Acta Materialia. 2015;**95**: 184-191

[54] Philip TV, Beck PA. CsCI-type ordered structures in binary alloys. Journal of Metals. 1957;**1957**: 1269-1271

[55] Jung J, Ghoshand G, Olson GB. A comparative study of precipitation behavior of Heusler phase (Ni2TiAl) from B2-TiNi in Ni–Ti–Al and Ni–Ti–Al– X (X = Hf, Pd, Pt, Zr) alloys. Acta Materialia. 2003;**51**:6341-6357

[56] John R, Ruben H. Theoretical investigations of Ti-based binary shape memory alloys. Materials Sciences and Applications. 2011;**2**:1355-1366

[57] Topor L, Kleppa OJ. Thermochemistry of the intermetallic compounds RuTi, RuZr, and RuHf. Metallurgical Transactions A. 1988;**9A**: 1061-1066

[58] Guo Q, Kleppa OJ. Standard enthalpies of formation of some alloys formed between group IV elements and group VIII elements, determined by high-temperature direct synthesis

calorimetry II. Alloys of (Ti, Zr, Hf) with (Co, Ni). Journal of Alloys and Compounds. 1998;**269**:181-186

[59] Togo A, Tanaka I. First principles phonon calculations in materials science. Scripta Materialia. 2015;**108**:1-5

[60] Buehler WJ, Gilfrich JW, Wiley RC. Effect of low temperature phase changes on the mechanical properties of alloys near composition TiNi. Journal of Applied Physics. 1963;**34**:1475

[61] Sedlák P, Janovská M, Bodnárová L, Heczko O, Seiner H. Softening of shear elastic coefficients in shape memory alloys near the martensitic transition: A study by laser-based resonant ultrasound spectroscopy. Meta. 2020;**10**:1-13

#### **Chapter 2**

## DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission

*Brahim Ait Hammou, Abdelhamid El Kaaouachi, El Hassan Mounir, Hamza Mabchour, Abdellatif El Oujdi, Adil Echchelh, Said Dlimi and Driss Ennajih*

#### **Abstract**

The Density-Functional Theory (DFT) is a reformulation of the quantum study of a correlated N-body system into a simpler system with independent equations being solved iteratively. The DFT considers only ground states of the systems. The extension to the time-dependent case of this theory is the Time-Dependent Density-Functional Theory (TDDFT) that also takes into account the excited states of the system. These calculations are very interesting in photonics areas. In fact, the interaction between electrons and light in the vicinity of solid surfaces and nanostructures is important as pathway to integrate photonics and electronics. The capability to couple light and electrons in purposefully designed device depends on the capability of creating such devices and the understanding of the underlying science.

**Keywords:** DFT, TDDFT, solid surfaces, ground states, excited states

#### **1. Introduction**

In the fields of photonics, researchers use the dielectric function and the electronic density in their calculations and investigations [1–8] to determine the physical and optical characteristics of materials such as noble metals (Au and Ag). In fact, the determination of the electronic structure of a material can give us all the physical and optical information about it. DFT and TDDFT are used to determine the electronic structure of materials such as noble metals in their ground and excited states. The purpose of these calculations is to determine the electron density *ρ r* ! . Indeed, it informs us about all the physical properties of the studied system. We simplify this very complex calculation by using several approximations and theorems. We propose at the end a flowchart to carry out a computer simulation allowing to highlight the values of the electron density *ρ r* ! .

#### **2. DFT calculation steps**

The DFT is based on wave Schrodinger equation

$$\hat{\mathcal{H}}\_T \mathcal{W}\left(\{\vec{\mathbb{R}}\_l\}, \{\vec{\mathbb{Z}\_l^\*}\}\_{l}\right) = i \frac{\partial}{\partial t} \mathcal{W}\left(\{\vec{\mathbb{R}}\_l\}, \{\vec{r\_i}\}\_{l}\right) \tag{1}$$

We consider that the system contains N cores and M electrons.

Η^ *<sup>T</sup>*: is the total Hamiltonian representing N cores and M electrons.

*<sup>ψ</sup> <sup>R</sup>* ! *I* � �, *ri* � � ! ,*<sup>t</sup>* � �: is the wave function representing N cores and M electrons.

*R* ! *I* n oand *ri* � � ! : represent respectively the set of nuclear and electronic coordinates.

In the case of stationary processes, the Schrodinger becomes:

$$\hat{\mathcal{H}}\_T \mathcal{\upmu} \left( \{ \vec{\kappa}\_l \}, \{ \overrightarrow{r\_i}^\cdot \} \right) = i \mathbf{E} \mathcal{\upmu} \left( \{ \vec{\kappa}\_l \}, \{ \overrightarrow{r\_i}^\cdot \} \right) \tag{2}$$

where E represent the energy of the system described by *<sup>ψ</sup> <sup>R</sup>* ! *I* � �, *ri* � � � � !

$$
\hat{\mathbf{H}}\_T = \hat{\mathbf{T}}\_T + \hat{\mathbf{V}}\_T \tag{3}
$$

where T^*<sup>T</sup>* is the total kinetic energy of the system operator, and *V*^ *<sup>T</sup>* is the operator describing all Colombian interaction.

Η^ *<sup>T</sup>* can be written as:

$$
\hat{\mathbf{H}}\_T = \hat{\mathbf{T}}\_n + \hat{\mathbf{T}}\_\epsilon + \hat{V}\_{n-\epsilon} + \hat{V}\_{\epsilon-\epsilon} + \hat{V}\_{n-n} \tag{4}
$$

where

<sup>T</sup>^*<sup>n</sup>* <sup>¼</sup> �ℏ<sup>2</sup> 2 P *I* ∇2*R* ! *I Mn* : is the kinetic energy of N cores with mass *Mn*. <sup>T</sup>^*<sup>n</sup>* <sup>¼</sup> �ℏ<sup>2</sup> 2 P *i* ∇<sup>2</sup> *r* !*i me* : is the kinetic energy of M electrons with mass *me*. *<sup>V</sup>*^ *<sup>n</sup>*�*<sup>e</sup>* <sup>¼</sup> <sup>1</sup> 4*πε*<sup>0</sup> P *i*,*j e*<sup>2</sup>*Zi R* ! *<sup>i</sup>*� *r* !*j* � � � � : is the core-electron attractive Colombian interaction. *<sup>V</sup>*^ *<sup>e</sup>*�*<sup>e</sup>* <sup>¼</sup> <sup>1</sup> 8*πε*<sup>0</sup> P *<sup>i</sup>*6¼*<sup>j</sup> <sup>e</sup>*<sup>2</sup> *r* ! *<sup>i</sup>*� *r* ! j j*<sup>j</sup>* : is the electron-electron repulsive Colombian interaction. *<sup>V</sup>*^ *<sup>n</sup>*�*<sup>n</sup>* <sup>¼</sup> <sup>1</sup> 8*πε*<sup>0</sup> P *i*6¼*j e*<sup>2</sup>*ZiZj R* ! *<sup>i</sup>*�*R* ! *j* � � � � : is the core-core repulsive Colombian interaction.

Eq. (4) takes into account N cores and M electrons. Some simplifications must be made:


*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

#### **2.1 Adiabatic born-Oppenheimer approximation**

It offers the opportunity to treat separately the cores and the electrons. So we have

$$
\Psi\_R\left(\stackrel{\rightarrow}{R},\stackrel{\rightarrow}{r}\right) = \Phi\_R\left(\stackrel{\rightarrow}{R}\right)\Psi\_R\left(\stackrel{\rightarrow}{r}\right) \tag{5}
$$

where

Φ*<sup>R</sup> R* � �! : is the wave function describing the cores. Ψ*<sup>R</sup> r* !� �: is the wave function describing the electrons. In this approximation, the interaction electron-phonon is neglected. Using Eqs. (3) and (4), the new Hamiltonian becomes

$$
\hat{\mathbf{H}}\_T = \hat{\mathbf{T}}\_\epsilon + \hat{\mathbf{V}}\_{n-\epsilon} + \hat{\mathbf{V}}\_{\epsilon-\epsilon} + \left(\hat{\mathbf{V}}\_{n-n} = \hat{\mathbf{V}}\_{\text{ext}} = \mathbf{C}\mathbf{st}\right) \tag{6}
$$

so

$$\dot{\mathbf{H}}\_{T} = \frac{-\hbar^{2}}{2} \sum\_{i} \frac{\nabla^{2} \overrightarrow{r}\_{i}}{m\_{e}} + \frac{1}{4\pi\varepsilon\_{0}} \sum\_{i,j} \frac{e^{2} Z\_{i}}{\left| \overrightarrow{R}\_{i} - \overrightarrow{r}\_{j} \right|} + \frac{1}{8\pi\varepsilon\_{0}} \sum\_{i \neq j} \frac{e^{2}}{\left| \overrightarrow{r}\_{i} - \overrightarrow{r}\_{j} \right|} \tag{7}$$
 
$$+ \left( \frac{1}{8\pi\varepsilon\_{0}} \sum\_{i \neq j} \frac{e^{2} Z\_{i} Z\_{j}}{\left| \overrightarrow{R}\_{i} - \overrightarrow{R}\_{j} \right|} = \text{Cst} \right)$$

The nuclear kinetic energy term independent from electrons is canceled (equal to zero). The attractive potential energy (electron-core) term is removed and the repulsive potential energy (core-core) becomes a constant simply evaluated for determined geometry. Using the Hohenberg and Kohn Theorems (HK theorems) [9], the formulation of the Schrödinger equation can be now based on the electron density *ρ r* !� �. This is due to the two HK theorems. In fact, the first KH theorem indicated that the total energy of the system in the ground state is a single universal function of the electronic density*.*

$$E = E\left[\rho\left(\overrightarrow{r}\right)\right] \tag{8}$$

HK gives

$$E\left[\rho\left(\overrightarrow{r}\right)\right] = F\_{HK}\left[\rho\left(\overrightarrow{r}\right)\right] + \int \hat{V}\_{\text{ext}}\left(\overrightarrow{r}\right)\rho\left(\overrightarrow{r}\right)d\overrightarrow{r},\tag{9}$$

where *FHK ρ r* ! h i � � represent the HK universal functional and *<sup>V</sup>* ! *ext r* !� � represent the external potential. The second HK theorem says that it is an analogue variational principle to the originally proposed in the approach of Hartree-Fock for the functional of the wave function (*<sup>∂</sup>E*½ � *<sup>ψ</sup> <sup>∂</sup><sup>ψ</sup>* ¼ 0), but this time applied to the electronic density functional.

$$\left. \frac{\partial E\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} \right|\_{\rho\_0\left(\overrightarrow{r}\right)} = 0 \tag{10}$$

$$G\left[\rho\left(\overrightarrow{r}\right)\right] = \int \rho\left(\overrightarrow{r}\right)d\overrightarrow{r} - \mathbf{M} \tag{11}$$

$$A\left[\rho\left(\overrightarrow{r}\right)\right] = E\left[\rho\left(\overrightarrow{r}\right)\right] - \mu G\left[\rho\left(\overrightarrow{r}\right)\right] \tag{12}$$

$$
\partial A \left[ \rho \left( \overrightarrow{r} \right) \right] = \int \frac{\partial A \left[ \rho \left( \overrightarrow{r} \right) \right]}{\partial \rho \left( \overrightarrow{r} \right)} d\rho dr = \mathbf{0} \tag{13}
$$

$$
\partial \left\{ E \left[ \rho \left( \overrightarrow{r} \right) \right] - \mu \left[ \left[ \rho \left( \overrightarrow{r} \right) d\overrightarrow{r} - M \right] \right] \right\} = \mathbf{0}
$$

$$\frac{\partial \mathcal{A}\left[\rho\left(\overrightarrow{\boldsymbol{r}}\right)\right]}{\partial \rho\left(\overrightarrow{\boldsymbol{r}}\right)} = \frac{\partial}{\partial \rho\left(\overrightarrow{\boldsymbol{r}}\right)} \left\{ E\left[\rho\left(\overrightarrow{\boldsymbol{r}}\right)\right] - \mu\left[\left[\rho\left(\overrightarrow{\boldsymbol{r}}\right)d\overrightarrow{\boldsymbol{r}} - \mathcal{M}\right]\right] \right\} = \frac{\partial E\left[\rho\left(\overrightarrow{\boldsymbol{r}}\right)\right]}{\partial \rho\left(\overrightarrow{\boldsymbol{r}}\right)} - \mu \frac{\partial}{\partial \rho\left(\overrightarrow{\boldsymbol{r}}\right)} \left[\left[\rho\left(\overrightarrow{\boldsymbol{r}}\right)d\overrightarrow{\boldsymbol{r}}\right] \right]$$

$$\frac{\partial \mathcal{A}\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} = \frac{\partial \mathcal{E}\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} - \mu\tag{14}$$

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

By replacing this expression in the equation of *∂A ρ r* ! h i � � , we obtained

$$
\partial A \left[ \rho \left( \overrightarrow{r} \right) \right] = \int \left[ \frac{\partial E \left[ \rho \left( \overrightarrow{r} \right) \right]}{\partial \rho \left( \overrightarrow{r} \right)} - \mu \right] \partial \rho dr = 0
$$

$$
\int \frac{\partial E \left[ \rho \left( \overrightarrow{r} \right) \right]}{\partial \rho \left( \overrightarrow{r} \right)} \partial \rho dr = \int \mu \partial \rho dr
\tag{15}
$$

$$
\int \frac{\partial E \left[ \rho \left( \overrightarrow{r} \right) \right]}{\partial \rho \left( \overrightarrow{r} \right)} \partial \rho dr = \mu
$$

Using Eq. (9) and calculating the functional derivative of *E ρ r* ! h i � � we obtained

$$\frac{\partial E\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} = \hat{V}\_{\text{ext}}\left(\overrightarrow{r}\right) + \frac{\partial F\_{HK}\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)}\tag{16}$$

By replacing Eq. (16) in Eq. (15) we obtained

$$\mu = \frac{\partial E\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} = \hat{\mathcal{V}}\_{\text{ext}}\left(\overrightarrow{r}\right) + \frac{\partial F\_{HK}\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)}\tag{17}$$

This equation constitutes the fundamental DFT formalism.

It only remains to determine the expression of the unknown function *FHK ρ r* ! h i � � in Eq. (17), that's why we are brought to use the Kohn-Sham approximations.

#### **2.2 Kohn-Sham (KS) approximations (equations)**

This step [10] consists of two approximations to transform the theorems of Hohenberg-Kohn into a practical workable from. However, the real system studied is redefined as a fictitious fermions system without interaction and with the same electronic density *ρ r* !� � characterizing the real system to being up the terms of inter-electronic as corrections of the other terms. Also, single particle orbitals are induced to treat the kinetic energy term of electrons more precisely than under the Thomas-Fermi theory.

#### *2.2.1 First approximation of KS*

The transformation of the real interactive system as a fictitious non-interactive system. Considering the first HK theorem, the functional *E ρ r* ! h i � � can be written as indicated in Eq. (9):

$$E\left[\rho\left(\overrightarrow{r}\right)\right] = F\_{HK}\left[\rho\left(\overrightarrow{r}\right)\right] + \int \hat{V}\_{\text{ext}}\left(\overrightarrow{r}\right)\rho\left(\overrightarrow{r}\right)d\overrightarrow{r} \tag{18}$$

The functional *FHK ρ r* ! h i � � is independent from *<sup>V</sup>*^ *ext <sup>r</sup>* !� � (external potential), and it is valid and applies regardless of the system studied. It contains a component of kinetic energy of electron *T*^*<sup>e</sup> ρ r* ! h i � � , and another component corresponding to the mutual coulomb interaction between electrons *<sup>V</sup>*^ *<sup>e</sup>*�*<sup>e</sup> <sup>ρ</sup> <sup>r</sup>* ! h i � � . The minimization of this functional with the constraint of preservation of the number of particles M:

Ð *ρ r* !� �*d r*! <sup>¼</sup> *<sup>M</sup>*gives us directly the total energy of the system and the charge density of ground state from which all other physical properties can be extracted:

$$F\_{HK}\left[\rho\left(\overrightarrow{r}\right)\right] = T\_S\left[\rho\left(\overrightarrow{r}\right)\right] + E\_H\left[\rho\left(\overrightarrow{r}\right)\right] + E\_{XC}\left[\rho\left(\overrightarrow{r}\right)\right] + V\_{ext}\left[\rho\left(\overrightarrow{r}\right)\right] \tag{19}$$

Where

*TS ρ r* ! h i � � : is the kinetic energy of non-interactive electron gas.

*EH ρ r* ! h i � � <sup>¼</sup> <sup>1</sup> 2 ÐÐ *<sup>ρ</sup> <sup>r</sup>* ! ð Þ*<sup>ρ</sup> <sup>r</sup>* !0 � � *r* !� *<sup>r</sup>* !0 � � � � *d r*!*d r*!0 : is the chemical coulomb interaction between electrons (Hartree term).

*EXC ρ r* ! h i � � : the additional functional describing interactions inter-electronic not obtained from non-interactive system.

*Vext ρ r* ! h i � � <sup>¼</sup> <sup>Ð</sup> *V*^ *ext r* !� �*<sup>ρ</sup> <sup>r</sup>* !� �*d r*!: is the external potential.

The functional *EXC ρ r* ! h i � � is called "exchange-correlation energy," this term contain all differences between fictitious non-interactive system and the real interactive system.

$$E\_{\rm XC}\left[\rho\left(\vec{r}\right)\right] = \left(T\left[\rho\left(\vec{r}\right)\right] - T\_S\left[\rho\left(\vec{r}\right)\right]\right) + \left(V\_{\epsilon-\epsilon}\left[\rho\left(\vec{r}\right)\right] - V\_H\left[\rho\left(\vec{r}\right)\right]\right) \tag{20}$$

where:

*T ρ r* ! h i � � : is the real kinetic energy (of real interactive system).

*TS ρ r* ! h i � � : is the non-interactive fermions system energy of KS.

*T ρ r* ! h i � � � *TS <sup>ρ</sup> <sup>r</sup>* ! h i � � : this difference is low and is neglected.

*EXC ρ r* ! h i � � traduced only the difference between Colombian energy of the real system *Ve*�*<sup>e</sup> ρ r* ! h i � � and the Colombian energy of non-interactive fermions system *VH ρ r* ! h i � � of KS.

$$E\_{\rm XC} \left[ \rho \left( \overrightarrow{r} \right) \right] \cong \mathbf{V}\_{e-e} \left[ \rho \left( \overrightarrow{r} \right) \right] - \mathbf{V}\_{H} \left[ \rho \left( \overrightarrow{r} \right) \right] \tag{21}$$

#### *2.2.2 Second approximation of KS*

This second approximation is based on the formulation of kinetic energy using an orbital approach. The exact formulation of kinetic energy *T* for the ground state systems is given by:

$$T = \sum\_{i}^{M} n\_{i\left<\rho\_{i}^{\*}\right|} \frac{-\hbar^{2}\nabla^{2}\overrightarrow{r}}{2m\_{\varepsilon}} |\rho\_{i}\rangle\tag{22}$$

$$T\_s \left[ \rho \left( \overrightarrow{r} \right) \right] = \sum\_{i}^{M} n\_{i \left< \rho\_i^\* \right|} \frac{-\hbar^2 \nabla^2 \overrightarrow{r}}{2m\_e} \left| \rho\_i \right> \tag{23}$$

$$\mu = \hat{V}\_{\text{eff}} \left( \overrightarrow{r} \right) + \frac{\partial T\_s \left[ \rho \left( \overrightarrow{r} \right) \right]}{\partial \rho \left( \overrightarrow{r} \right)} \tag{24}$$

$$
\hat{\mathbf{V}}\_{\text{eff}}\left(\overrightarrow{\mathbf{r}}\right) = \hat{\mathbf{V}}\_{\text{eff}}\left[\rho\left(\overrightarrow{\mathbf{r}}\right)\right] = \hat{\mathbf{V}}\_{\text{ext}}\left(\overrightarrow{\mathbf{r}}\right) + \frac{\partial E\_{\text{H}}\left[\rho\left(\overrightarrow{\mathbf{r}}\right)\right]}{\partial\rho\left(\overrightarrow{\mathbf{r}}\right)} + \frac{\partial E\_{\text{XC}}\left[\rho\left(\overrightarrow{\mathbf{r}}\right)\right]}{\partial\rho\left(\overrightarrow{\mathbf{r}}\right)}
$$

$$
\hat{\mathbf{V}}\_{\text{eff}}\left(\overrightarrow{\mathbf{r}}\right) = \hat{\mathbf{V}}\_{\text{ext}}\left(\overrightarrow{\mathbf{r}}\right) + \int \frac{\rho\left(\overrightarrow{\mathbf{r}}\right)}{\left|\overrightarrow{\mathbf{r}} - \overrightarrow{\mathbf{r}}\right|} d\overrightarrow{\mathbf{r}}\frac{\mathbf{r}}{\left|\overrightarrow{\mathbf{r}}\right|} + \hat{\mathbf{V}}\_{\text{xc}}\left(\overrightarrow{\mathbf{r}}\right) \tag{25}
$$

$$\hat{\mathcal{V}}\_{\text{xc}}\left(\overrightarrow{r}\right) = \frac{\partial \mathcal{E}\_{\text{XC}}\left[\rho\left(\overrightarrow{r}\right)\right]}{\partial \rho\left(\overrightarrow{r}\right)} \text{ and } \rho\left(\overrightarrow{r}\right) = \sum\_{i=1}^{M} \left|\rho\_{i}\left(\overrightarrow{r}\right)\right|^{2} \tag{26}$$

$$\left[\frac{-\hbar^2\nabla^2\vec{r}}{2m\_e} + \hat{V}\_{\text{eff}}\left(\vec{r}\right)\right] \left|\rho\_i\left(\vec{r}\right)\right> = \varepsilon\_i \left|\rho\_i\left(\vec{r}\right)\right>\tag{27}$$

The KS equations (Eqs. (25) and (27)) must be resolved with self-coherent method by starting with an initial electronic density. *V*^ *eff r* !� � is obtained for which the Schrodinger equation of Kohn-Sham (Eq. 27) is resolved, and a new electronic density is then calculated. From this new electronic density, a new *V*^ *eff r* !� � is determined. This process is repeated with self-coherent method up to that convergence is reached. The new electronic density obtained must be very close to the previous one corresponding to the criterion of convergence fixed before.

At this step, we only need to find the expression of *EXC ρ r* ! h i � � .Some exchangecorrelation functional can be considered using different functional approximation families like: Local Density Approximation (LDA) [11], Generalized Gradient Approximation (CGA) [12], Meta Generalized Gradient Approximation (Meta CGA) [13], and Hybrid Functional [14].


The order of accuracy is increasing from top to bottom. In the formalism of exchange-correlation functional, *EXC ρ r* ! h i � � is presented like an interaction between the electronic density *ρ r* !� � and density energy depending on *<sup>ρ</sup> <sup>r</sup>* !� �: *εXC ρ r* ! h i � � with

$$E\_{\rm XC} \left[ \rho \left( \overrightarrow{r} \right) \right] = \int \varepsilon\_{\rm XC} \left[ \rho \left( \overrightarrow{r} \right) \right] \rho \left( \overrightarrow{r} \right) d\overrightarrow{r} \tag{28}$$

*εXC ρ r* ! h i � � is considered as a summation of the contribution of exchange and correlation, (with) *εXC ρ r* ! h i � � <sup>¼</sup> *<sup>ε</sup><sup>X</sup> <sup>ρ</sup> <sup>r</sup>* ! h i � � <sup>þ</sup> *<sup>ε</sup><sup>C</sup> <sup>ρ</sup> <sup>r</sup>* ! h i � � and

$$E\_{\rm XC}\left[\rho\left(\vec{r}\right)\right] = E\_{\rm X}\left[\rho\left(\vec{r}\right)\right] + E\_{\rm C}\left[\rho\left(\vec{r}\right)\right] = \int \varepsilon\_{\rm X}\left[\rho\left(\vec{r}\right)\right] \rho\left(\vec{r}\right) d\vec{r} + \int \varepsilon\_{\rm C}\left[\rho\left(\vec{r}\right)\right] \rho\left(\vec{r}\right) d\vec{r} \tag{29}$$

#### **3.1 Local Density Approximation (LDA)**

In this work we will only use the LDA approximation in our calculations. In this approximation, the electronic density can be treated locally in the form of uniform electrons of gas.

In other words, this approach is to perform the following two hypotheses:


The fundamental hypothesis contained in the LDA formalism is to consider that the contribution of *EXC ρ r* ! h i � � to the total energy of the system can be added by cumulated method of each portion of the non-uniform gas as it was locally uniform.

$$E\_{\rm XC}^{LDA}\left[\rho\left(\overrightarrow{r}\right)\right] = \int \epsilon\_{\rm XC}^{LDA}\left[\rho\left(\overrightarrow{r}\right)\right] \rho\left(\overrightarrow{r}\right) d\overrightarrow{r} \tag{30}$$

Where *εLDA XC ρ r* ! h i � � represent the exchange-correlation energy by electron in a system with electron in mutual correlation with uniform density *ρ r* !� �.

Using *εLDA XC ρ r* ! h i � � , the exchange-correlation potential *<sup>V</sup>LDA XC r* !� � can be obtained like:

$$V\_{\rm XC}^{LDA}\left(\overrightarrow{r}\right) = \frac{\partial \left(\rho\left(\overrightarrow{r}\right)\varepsilon\_{\rm XC}^{LDA}\left[\rho\left(\overrightarrow{r}\right)\right]\right)}{\partial \rho\left(\overrightarrow{r}\right)}\tag{31}$$

There are several forms for the term of exchange and correlation of a homogeneous electron gas, among others those of Kohn and Sham [10], Wigner [15], Perdew and Wang [16], and Hedin and Lundqvis [17], Ceperly and Alder [18]. The last one is the most used.

#### **3.2 Approximation of Ceperley and Alder**

Ceperley and Alder [18] used the exchange energy of uniform electrons gas given by Dirac formula as:

$$E\_X\left[\rho\left(\overrightarrow{r}\right)\right] = -C\_X\int \rho^\sharp \left(\overrightarrow{r}\right)d\overrightarrow{r}\tag{32}$$

$$
\epsilon\_X^{LDA} \left[ \rho \left( \overrightarrow{r} \right) \right] = -\mathcal{C}\_X \int \rho^\sharp \left( \overrightarrow{r} \right) d\overrightarrow{r} \tag{33}
$$

with *CX* <sup>¼</sup> <sup>3</sup> 4 3 *π* � �<sup>1</sup> 3 , *εLDA <sup>X</sup> ρ r* ! h i � � can be expressed like *εLDA <sup>X</sup> ρ r* ! h i � � ¼ �0*:*738694*<sup>ρ</sup>* 1 <sup>3</sup> *r* !� �.

The correlation energy *εLDA <sup>C</sup> ρ r* ! h i � � is derived on the second order of Moller-Plesset perturbation theory [19–21]:

$$
\epsilon\_C^{LDA} \left[ \rho \left( \overrightarrow{r} \right) \right] = aLn \left( \mathbf{1} + \frac{b}{r\_s} + \frac{b}{r\_s^2} \right) \tag{34}
$$

With *<sup>a</sup>* <sup>¼</sup> *Ln*2�<sup>1</sup> <sup>2</sup>*π*<sup>2</sup> ¼ �0*:*01556111 and *b* ¼ 20*:*4562557.

*rs*is the Wigner-Seitz density parameter, in atomic unit we have *rs* <sup>¼</sup> <sup>4</sup>*πρ <sup>r</sup>* ! ð Þ 3 � ��1*=*<sup>3</sup> .

Considering the correlation term *εLDA <sup>C</sup> ρ r* ! h i � � , no explicit analytical expression is known. Several different parameterizations have been proposed since the early 1970's, and the most accurate results are based on quantum Monte Carlo simulations of Ceperley and Alder [18]. The most approximation commonly used today is that of Zunger [22] and who parameterized the Ceperley and Alder correlation functional for the spin-polarized electron gas and non-spin-polarized electron gas by the equation:

$$\varepsilon\_C^{pZ} \left[ \rho \left( \overrightarrow{r} \right) \right] = \begin{cases} A \times Lnr\_s + B + C \times r\_s \times Lnr\_s + D \times r\_s, r\_s \le 1 \\ \gamma/(1 + \beta\_1 \times r\_s + \beta\_2 \times r\_s), r\_s < 1 \end{cases} \tag{35}$$

For *rs* <sup>≤</sup>1, *<sup>A</sup><sup>U</sup>* <sup>¼</sup> <sup>0</sup>*:*0311, *BU* ¼ �0*:*048, *CU* <sup>¼</sup> <sup>0</sup>*:*002 and *DU* ¼ �0*:*0116. For *rs* <1, *γ* ¼ �0*:*1423, *β*<sup>1</sup> ¼ 1*:*0529, and *β*<sup>2</sup> ¼ 0*:*3334.

The improvements of the LDA approach must consider the gas of electrons in its real form (non-uniform and non-local), GGA, meta-GGA, and hybrid functional allow gradually approach to make in consideration their two effects.

#### **4. Resolution of the KS Schrödinger equation**

Combining between Eqs. (25) and (27) and using the expression *ρ r* !� � <sup>¼</sup> <sup>P</sup> *M i*¼1 *φ<sup>i</sup> r* ! � � � � � � � � 2 , where M is the number of electrons, Eq. (27) becomes:

$$\left[\frac{-\hbar^{2}\nabla^{2}\vec{r}}{2m\_{e}}+\hat{V}\_{ext}\left(\vec{r}\right)+\left[\frac{\rho\left(\vec{r}'\right)}{\left|\vec{r}-\vec{r}'\right|}d\vec{r}'+\hat{V}\_{xc}\left(\vec{r}\right)\right]\Big|\rho\_{i}\left(\vec{r}\right)\right>=\varepsilon\_{i}\left|\rho\_{i}\left(\vec{r}\right)\right>\tag{36}$$

The only unknown value in this expression is *V*^ *xc r* !� � (exchange-correlation potential).

Using LDA approximation of *V*^ *xc r* !� � (Eq. (30)),Where *<sup>ε</sup>LDA xc ρ r* ! h i � � <sup>¼</sup> *εLDA <sup>x</sup> ρ r* ! h i � � <sup>þ</sup> *<sup>ε</sup>LDA <sup>c</sup> ρ r* ! h i � � with *εLDA xc ρ r* ! h i � � : The exchange-correlation energy of one electron in a system of electron on mutual interaction with electronic density *ρ r* !� �. *εLDA <sup>x</sup> ρ r* ! h i � � : The exchange energy of one electron. *εLDA <sup>c</sup> ρ r* ! h i � � : The correlation energy of one electron. Eq. (36) becomes:

$$\left[\frac{-\hbar^{2}\nabla^{2}\vec{r}}{2m\_{\epsilon}}+\hat{V}\_{\text{ext}}\left(\vec{r}\right)+\left[\frac{\rho\left(\vec{r}'\right)}{|\vec{r}-\vec{r}'|}d\vec{r}'+\frac{\partial\left(\rho\left(\vec{r}\right)\epsilon\_{\text{xc}}^{LDA}\left[\rho\left(\vec{r}\right)\right]\right)}{\partial\rho\left(\vec{r}\right)}\right]\middle|\rho\_{i}\left(\vec{r}\right)\right>=\epsilon\_{i}\left|\rho\_{i}\left(\vec{r}\right)\right>\tag{37}$$

$$\begin{aligned} & \left[ \frac{-\hbar^2 \nabla^2 \vec{r}}{2m\_e} + \dot{V}\_{ext} \left( \vec{r} \right) + \int \frac{\rho \left( \vec{r'} \right)}{|\vec{r} - \vec{r'}|} d\vec{r'} + \frac{\partial \left( \rho \left( \vec{r'} \right) \left( \epsilon\_x^{LDA} \left[ \rho \left( \vec{r} \right) \right] \right) + \epsilon\_c^{LDA} \left[ \rho \left( \vec{r'} \right) \right] \right)}{\partial \rho \left( \vec{r} \right)} \right] \phi\_i \left( \vec{r'} \right) \end{aligned} \tag{20}$$
 
$$\begin{aligned} &= \varepsilon\_i \left| \rho\_i \left( \vec{r'} \right) \right> \end{aligned} \tag{20}$$

$$\begin{aligned} &\left[-\frac{\hbar^2 \nabla^2 \vec{r}}{2m\_\varepsilon} + \hat{V}\_{\rm ext} \left(\vec{r}\right) + \int \frac{\rho \left(\vec{r}'\right)}{\left|\vec{r} - \vec{r}'\right|} d\vec{r}' + \frac{\partial \left(\rho \left(\vec{r}\right) \left(\epsilon\_x^{LDA} \left[\rho \left(\vec{r}\right)\right]\right)\right)}{\partial \rho \left(\vec{r}\right)}\right) \\ &+ \frac{\partial \left(\rho \left(\vec{r}\right) \left(\epsilon\_x^{LDA} \left[\rho \left(\vec{r}\right)\right]\right)\right)}{\partial \rho \left(\vec{r}\right)}\right) \Big|\rho\_i \left(\vec{r}\right)\rangle \\ &= e\_i \left|\rho\_i \left(\vec{r}\right)\right| \end{aligned} \tag{39}$$

#### **4.1 Self-coherent technical to resolve KS Schrödinger equation**

The idea is that do not directly resolve Eq. (39), but to write previously the *ϕ<sup>m</sup> r* !� � in a finite basis function *ϕ<sup>b</sup> <sup>p</sup> r* !� � as:

$$\phi\_m(\overrightarrow{r'}) = \sum\_{p=1}^p \mathbf{C}\_p^m \phi\_p^b(\overrightarrow{r'}) \tag{42}$$

when *m* = *n.k* ! , *k* ! : is the wave vector belonging the first Brillouin zone in the case of crystal lattice.

The resolution of Eq. (38) consists to determine the coefficients *C<sup>m</sup> <sup>p</sup>* necessary to express *ϕ<sup>m</sup> r* !� � in a given basis *<sup>ϕ</sup><sup>b</sup> <sup>p</sup> r* !� �.

We need to search a basis to get as close as possible to *ϕ<sup>m</sup> r* !� � with *<sup>p</sup>* having a finite value.

Eq. (37) becomes:

$$
\begin{bmatrix}
\cdots & & & & \cdots \\
\cdots & \langle \phi\_i^b \middle| \hat{H}\_{\text{KS}} \middle| \phi\_j^b \rangle - \epsilon\_m \langle \phi\_i^b \middle| \phi\_j^b \rangle & \cdots \\
& & & \cdots \\
\cdots & & & \cdots \\
\end{bmatrix}
\begin{bmatrix} C\_1^m \\ \cdots \\ C\_p^m \end{bmatrix} = \begin{bmatrix} \mathbf{0} \\ \dots \\ \mathbf{0} \end{bmatrix}, \tag{43}
$$

in which one can identify the matrix Elements of Hamiltonian of single particle and the elements of the recovery matrix, *Hij* � *<sup>ε</sup>mSij* � �*C<sup>m</sup> <sup>p</sup>* <sup>¼</sup> 0, where *<sup>H</sup>*^ *ij* <sup>¼</sup>

*ϕb i* � � �*H*^ *KS <sup>ϕ</sup><sup>b</sup> j* � � � E and *Sij* <sup>¼</sup> *<sup>ϕ</sup><sup>b</sup> i* � *ϕ<sup>b</sup> j* � � � E respectively represent the Hamiltonian matrix and the recovery matrix.

For a solid, these equations need to be resolved for each vector *k* in the Brillouin zone. This secular equation system is linear with the energy. This system constitutes a problem of determination of the proper values *ε<sup>m</sup>* and proper functions *ϕ<sup>k</sup> <sup>i</sup> r* !� � that much we know within the Hartree theory and is commonly solved from standard numerical methods. The diagonalization of the Hamiltonian matrix provides *p* proper values clean and *p* sets values of coefficients that express each *p* proper functions in a given basis.

More *p* is big, more the proper functions are precise, but the matrix diagonalization time is also particularly high.

#### **4.2 Self-consistent cycle**

Eq. (39) must be resolved in an iterative way in self-consistent cycle procedure. The procedure starts by the definition of a density of departure corresponding to a determined geometry core. Generally, the initial density is constituted from a superposition of atomic densities: *<sup>ρ</sup>in* <sup>¼</sup> *<sup>ρ</sup>crystal* <sup>¼</sup> <sup>P</sup> *atρat*.

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

When the elements of the Hamiltonian matrix and recovery matrix were calculated, proper vectors and proper functions are determined from the diagonalization of the matrix *Hij* � *<sup>ε</sup>mSij* � �*Cm <sup>p</sup>* ¼ 0. Following the principle of Aufbau, the orbitals are filled, a new density is determined:

$$
\overline{\rho^{\rm out}\left(\overrightarrow{r}\right) = \sum\_{out} \left| \phi\_i \left(\overrightarrow{r}\right) \right|^2}.
$$

This step concludes the first cycle. At this stage of the process, acceleration of convergence is generally used to generate a new density realized from a mixture between this output density *<sup>ρ</sup>out*ð Þ*<sup>r</sup>* and density of entry of this cycle and there *<sup>ρ</sup>in*ð Þ*<sup>r</sup>* . One of the simplest procedures concerning this mixture can be formulated as: *ρin <sup>i</sup>*þ<sup>1</sup>ð Þ¼ *<sup>r</sup>* ð Þ <sup>1</sup> � *<sup>α</sup> <sup>ρ</sup>in <sup>i</sup>* ð Þþ *<sup>r</sup> αρout <sup>i</sup>* ð Þ*r* with 0≺ *α*≺1. *α* is the mixture parameter an *i* corresponds to iterative cycle number. Density of the new entry *ρin <sup>i</sup>*þ<sup>1</sup>ð Þ*r* is then introduced into a second self-consistent cycle. This process is repeated for an iterative manner until it has the convergence criterion (difference between *<sup>ρ</sup>out*ð Þ*<sup>r</sup>* and *<sup>ρ</sup>in*ð Þ*<sup>r</sup>* ) initially fixed is reached. When convergence is reached, the energy of the ground state of the system is reached (**Figure 1**).

Once we obtain *ρout r* !� � for the first cycle, we inject it like a new value of *<sup>ρ</sup>in <sup>r</sup>* !� � and we repeat this iterative operation until it is that the fixed precision at the beginning is reached.

**Figure 1.**

*Schematic representation of the self-consistent cycle.*

\$\text{Precision} = \$\rho^{out}\left(\overrightarrow{r}\right) - \rho^{in}\left(\overrightarrow{r}\right)\$

$$\rho^{in}\left(\overrightarrow{r}\right) = \rho^{orbital}\left(\overrightarrow{r}\right) = \sum\_{\text{at}} \rho^{\text{at}}\left(\overrightarrow{r}\right),$$

We can use *ϕ<sup>n</sup> <sup>k</sup> r* !� � <sup>¼</sup> *<sup>e</sup>ik* ! *:r* ! *un <sup>k</sup> r* !� �, where *<sup>e</sup>i k* ! *:r* ! is a plane wave and *u<sup>n</sup> <sup>k</sup> r* !� � a function processing Lattice periodicity (Bloch theorem *u<sup>n</sup> <sup>k</sup> r* ! þ *k* � �! <sup>¼</sup> *<sup>u</sup><sup>n</sup> <sup>k</sup> r* !� �).

#### **5. Time dependent density functional theory (TDDFT) calculations**

The Eigen values of Kohn-Sham does not match corresponding to the required energy for excited electrons in a correlated electron system with N electrons, because the Kohn-Sham [10] approach replaces the interacting electron problem by an independent electron problem. There is a way to study the potential and the wave function to determine the excitation energies. It is the use of the TDDFT theory which takes into account the time parameter. The excitation energy of a system can then be obtained from the response function to extreme electron density perturbations.

#### **5.1 Rung-Gross theorems**

TDDFT is based on the Rung-Gross theorem [23, 24] that is the analog of the timedependent of Hohenberg-Kohn theorems [9]. The theorems are cited below:

#### *5.1.1 First theorem*

The time depends on the external potential *V*^ *ext r* !, *t* � � is determined by the timedependent electron density *ρ r* !, *t* � � to a nearly additive function for an initial state *<sup>ψ</sup>* <sup>a</sup> state*ψ*ð Þ *t* ¼ 0 .

#### *5.1.2 Second theorem*

By difference with the ground state, the variational principle, which stated that there is a minimum associated with the total energy does not exist for the timedependent systems because the energy is not a conserved quantity. A similar amount of energy, which is applied on the stationary principle, is defined as:

$$\mathbf{A}[\boldsymbol{\psi}] = \int\_{t\_0}^{t\_1} dt \langle \boldsymbol{\psi}(t) | \boldsymbol{i} \frac{\partial}{\partial t} - \hat{H}(t) | \boldsymbol{\psi}(t) \rangle \tag{44}$$

where Α is called the action, and *ψ*ð Þ*t* is a function of the time-dependent poly electronic wave function.

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

#### **5.2 The time-dependent approach of Kohn-Sham**

The time-dependent approach [10] of Kohn-Sham shows an equation with partial derivatives for the effective system such as

$$i\frac{\partial}{\partial t}\rho\_i\left(\overrightarrow{r},t\right) = H\rho\_i\left(\overrightarrow{r},t\right) = \left[-\frac{1}{2}\nabla^2 + \hat{V}\_{\text{eff}}\left(\rho\left(\overrightarrow{r},t\right)\right)\right]\rho\_i\left(\overrightarrow{r},t\right) \tag{45}$$

$$\rho\left(\overrightarrow{r},t\right) = \sum\_{i=1}^{N} \left| \rho\_i \left(\overrightarrow{r},t\right) \right|^2 \tag{46}$$

With *φ<sup>i</sup> r* !, *t* � � is a single electron time-dependent wave function of the Kohn-Sham. *V*^ *eff ρ r* !, *t* � � � � defined in Eq. (46) contains in addition to the external field, a Hartree potential, and the exchange-correlation potential such as:

$$
\hat{\mathbf{V}}\_{\text{eff}}\left(\rho\left(\overrightarrow{r},t\right)\right) = \hat{\mathbf{V}}\_{\text{ext}}\left(\overrightarrow{r},t\right) + \frac{\partial E\_H\left(\overrightarrow{r},t\right)}{\partial \rho\left(\overrightarrow{r},t\right)} + \hat{\mathbf{V}}\_{\text{xc}}\left(\rho\left(\overrightarrow{r},t\right)\right) \tag{47}
$$

The exchange-correlation term *V*^ *xc ρ r* !, *t* � � � � is unknown. There is a major difficulty to evaluate it because it does not depend only on the density *ρ r* !, *t* � �, but also depends on the prior density at time *t* in every point in space. The evolution of *V*^ *xc ρ r* !, *t* � � � � therefore requires the introduction of the single fundamental approximation of TDDFT.

#### **6. Linear response theory**

The TDDFT approach allows one hand to study the perturbation of the system at time *t*0, secondly to propagate this disturbance for a time *t*> *t*0. The study of the evolution in the propagation of this disturbance leads to the production of the absorption spectrum. This theory [25, 26] is used in the space of frequencies ω rather than in the temporal space. The operation is performed by means of the Fourier transform.

Considering, the continuous function *f* of the time variable. The Fourier transform of *f* gives:

$$\text{TF}(f): a \mapsto f(a) = \frac{1}{\sqrt{2\pi}} \int\_{-\infty}^{+\infty} f(t)e^{-i\alpha t}dt\tag{48}$$

Considering a system in the ground state of electron density *ρ*<sup>0</sup> which is subject at *<sup>t</sup>* <sup>¼</sup> *<sup>t</sup>*<sup>0</sup> to a low disturbance of the external potential *<sup>∂</sup>V*^ *ext* the electron density obtained following this disturbance is, at first order given by:

$$
\rho\left(\overrightarrow{r}, w\right) = \rho\_0\left(\overrightarrow{r}\right) + \delta\rho\left(\overrightarrow{r}, w\right) \tag{49}
$$

In the interacting electrons system, the external potential is also written as the sum of the external potential calculated in the ground state, supplemented with a disruptive potential:

$$
\hat{\mathcal{V}}\_{\text{ext}}\left(\overrightarrow{r},\,\boldsymbol{\alpha}\right) = \hat{\mathcal{V}}\_{\text{ext}}\left(\overrightarrow{r}\right) + \delta\hat{\mathcal{V}}\_{\text{ext}}\left(\overrightarrow{r},\,\boldsymbol{\alpha}\right) \tag{50}
$$

The perturbation *δρ* depends only on the potential *δVext* and takes the expression:

$$
\delta\rho\left(\overrightarrow{r},\alpha\right) = \int\_{r'} \chi\left(\overrightarrow{r},\overrightarrow{r'},\alpha\right) \delta\hat{V}\_{\text{ext}}\left(\overrightarrow{r'},\alpha\right) d\overrightarrow{r'}\tag{51}
$$

*χ r* !, *r* !0 ,*ω* � �: represents the response function of the non-interacting electrons system. The evolution of *χ* is complex, but it can be simplified by using the Kohn-Sham approach that redefines *δρ r* !, *ω* � � according to the equation:

$$
\delta\rho\left(\overrightarrow{r},\,\omega\right) = \int\_{r'} \chi\_{\text{KS}}\left(\overrightarrow{r},\,\overrightarrow{r'},\,\omega\right) \delta\hat{V}\_{\text{eff}}\left(\overrightarrow{r'},\,\omega\right) d\overrightarrow{r'}\tag{52}
$$

where *χKS r* !, *r* !0 ,*ω* � �, according to Kohn-Sham response is easily calculated using the equation:

$$\chi\_{\rm KS}\left(\overrightarrow{r},\overrightarrow{r}',\alpha\right) = \lim\_{\eta \to 0^{+}} \sum\_{p=1}^{\infty} \sum\_{q=1}^{\infty} \left(f\_p - f\_q\right) \frac{\rho\_p^{\*\ast}\left(\overrightarrow{r}\right)\rho\_q\left(\overrightarrow{r}\right)\rho\_q^{\*\ast}\left(\overrightarrow{r}'\right)\rho\_p\left(\overrightarrow{r}'\right)}{\alpha - \alpha\_{p\to q} + i\eta} \tag{53}$$

This expression (52) represents the various excitations between the occupied Kohn-Sham orbital's *φ<sup>p</sup>* and unoccupied *φq*, for a total number N of orbital's approaching infinity. The occupation of the orbital's number *p* and *q* are respectively denoted *f <sup>p</sup>* and *f <sup>q</sup>*. The frequency *ω<sup>p</sup>*!*<sup>q</sup>* ¼ *ε<sup>p</sup>* � *εq*; where *ε<sup>p</sup>* and *ε<sup>q</sup>* are the eigenvalues respectively associated with wave functions *φ<sup>p</sup>* and *φq*. The infinitesimal positive number *η* is introduced to reflect the causality of the system response.

Disruption of the effective potential *δV*^ *eff* (Kohn-Sham potential) can be written as the sum of three terms:


Such:

$$
\delta \hat{V}\_{\rm eff} \left( \overrightarrow{r}, \alpha \right) = \delta \hat{V}\_{\rm ext} \left( \overrightarrow{r}, \alpha \right) + \int\_{r'} \frac{\delta \rho \left( \overrightarrow{r'}, \alpha \right)}{\left| \overrightarrow{r} - \overrightarrow{r'} \right|} d\overrightarrow{r'} + \delta \hat{V}\_{\rm xc} \left( \overrightarrow{r}, \alpha \right) \tag{54}
$$

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

Following a Fourier transformation, it is customary to write disturbance of the exchange-correlation potential as:

$$
\delta \hat{V}\_{\text{xc}} \left( \overrightarrow{r}, \alpha \right) = T F \left( \int\_{r'} \left[ f\_{\text{xc}}(r, t, r', t') \delta \rho(r', t') dt' dr' \right] \right) \tag{55}
$$

Where

$$f\_{\infty}(r, t, r', t') = \frac{\delta \hat{V}\_{\infty}\left(\overrightarrow{r}, a\right)}{\delta \rho(r', t')},\tag{56}$$

*f xc r*, *t*,*r*<sup>0</sup> , *t* <sup>0</sup> ð Þ is the expression of exchange and correlation core.

In this step of the calculation, the development of the disturbance of the electron density (Eq. 51) through the expression of the effective potential disturbance (Eq. 52) gives:

$$\delta\rho\left(\overrightarrow{r},\,\boldsymbol{w}\right) = \int\_{\boldsymbol{r}'} \boldsymbol{\chi\_{\rm KS}}\left(\overrightarrow{r},\,\overrightarrow{r}',\,\boldsymbol{w}\right) \delta\boldsymbol{\mathcal{V}}\_{\rm ext}\left(\overrightarrow{r}',\,\boldsymbol{w}\right) d\overrightarrow{r}' + \int\_{\boldsymbol{r}\_1} \int\_{\boldsymbol{r}\_2} \boldsymbol{\chi\_{\rm KS}}\left(\overrightarrow{r},\,\overrightarrow{r}\_1,\,\boldsymbol{w}\right) \boldsymbol{f}\_{\rm Hxc}\left(\overrightarrow{r}\_1,\,\overrightarrow{r}\_2,\,\boldsymbol{w}\right) \delta\rho\left(\overrightarrow{r}\_2,\,\boldsymbol{w}\right) d\boldsymbol{r}\_2 d\boldsymbol{r}\_1 \tag{57}$$

Where

$$f\_{H\text{xc}}\left(\overrightarrow{r},\overrightarrow{r}',o\right) = \frac{1}{\left|\overrightarrow{r}-\overrightarrow{r}'\right|} + f\_{\text{xc}}\left(\overrightarrow{r},\overrightarrow{r}',o\right) \tag{58}$$

*f Hxc r* !, *r* !0 , *ω* � � is the expression of the Hartree exchange-correlation core.

Expression of the disturbance of the electron density is known both in the interacting electron system (Eq. 57) and non-interacting electrons system (Eq. 51). The equality of these two expressions leads to that of the response function of the interacting electron system:

$$\chi\left(\overrightarrow{r},\overrightarrow{r}',o\right) = \chi\_{\text{KS}}\left(\overrightarrow{r},\overrightarrow{r}',o\right) + \int\_{r\_1}\int\_{r\_2}\chi\left(\overrightarrow{r},\overrightarrow{r}\_1,o\right)f\_{\text{Hxc}}\left(\overrightarrow{r}\_1,\overrightarrow{r}\_2,o\right)\chi\_{\text{KS}}\left(\overrightarrow{r}\_2,\overrightarrow{r}',o\right)dr\_2dr\_1\tag{59}$$

Equation of the susceptibility *χ* is important because this linear response has poles exactly to the energies of the electronic transitions of the system.


#### **7. Resolution of time-dependent problem**

#### **7.1 Adiabatic approximation**

No analytical expression exists for the time-dependent exchange correlation potential. So, we work using the adiabatic approximation [26] by ignoring the dependence on previous electron densities of these two terms and use only the instantaneous electron density. Furthermore, when the time-dependent potential changes slowly, thus adiabatically, the system remains in its instantaneous ground state. The exchange-correlation potential developed in the DFT to describe the ground state of the system can then be translated to time-dependent systems. The potential and the adiabatic exchange-correlation core then take local forms in time given to equations (Eq. 60) and (Eq. 59). It is possible to note that in this approximation, the exchangecorrelation core is independent of the frequency of the disturbance potential applied to the system.

$$
\hat{\mathcal{V}}\_{\text{xc}}^{\text{adi}}\left(\overrightarrow{r},t\right) = \hat{\mathcal{V}}\_{\text{xc}}^{DFT}\left[\rho\left(\overrightarrow{r}\right)\right]\_{\rho\left(\overrightarrow{r}\right) = \rho\left(\overrightarrow{r},t\right)}\tag{60}
$$

$$f\_{\rm xc}^{\rm adia}\left(\overrightarrow{r},t,r',t'\right) = \frac{\delta \hat{V}\_{\rm xc}\left(\overrightarrow{r}\right)}{\delta \rho\left(\overrightarrow{r}'\right)}\delta(t-t')\tag{61}$$

The simplicity of this approximation allows to make after Fourier transform, the exchange-correlation core *f adia xc* independent of the frequency *ω*. Therefore, *f Hxc* also becomes independent of the frequency as:

$$\hat{f}^{\text{adia}}\_{H\text{xx}}\left(\overrightarrow{r},\overrightarrow{r}',\alpha\right) = f^{\text{adia}}\_{H\text{xx}}\left(\overrightarrow{r},\overrightarrow{r}'\right) = \frac{\textbf{1}}{\left|\overrightarrow{r}-\overrightarrow{r}'\right|} + f^{\text{adia}}\_{\text{xx}}\left(\overrightarrow{r},\overrightarrow{r}'\right) \tag{62}$$

In the adiabatic approximation, using functional developed for DFT retains weaknesses of these functional. For example, functional LDA and most functional GGA falsely represent the decay of exchange-correlation potential in neutral finite systems.

#### **7.2 Equation with Eigenvalues**

When the exchange-correlation core was independent of the frequency (in the adiabatic approximation case), then search for the poles of the response factor by solving a system of equations to the eigenvalues [22, 26] of a matrix equation as the specific form:

$$\sum\_{p'=1}^{\infty} \sum\_{q'=1}^{\infty} \mathbf{M}\_{pq, p'q'} \left[ f\_{H\mathbf{x}c}^{adiai} \right] \mathbf{X}\_{p'q'} = \mathfrak{Q}\_{pq}^2 \mathbf{X}\_{pq} \tag{63}$$

where *p*, *p*<sup>0</sup> et *q*, *q*<sup>0</sup> respectively represent indices of non-occupied and occupied orbitals for the orbital's base with size tends to infinity. *Mpq*,*p*<sup>0</sup> *<sup>q</sup>*<sup>0</sup> is an element of the operator matrix writing such that:

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

$$\mathcal{M}\_{pq,p'q'}\left[\mathcal{F}\_{\text{Hxc}}^{\text{data}}\right] = \mathcal{D}\mathcal{W}\_{pq,p'q'}\left[\mathcal{F}\_{\text{Hxc}}^{\text{data}}\right]\sqrt{\alpha\_{p\to q}\alpha\_{p'\to q'}} + \alpha\_{p\to q}^2\delta(pq)(p'q'),\tag{64}$$

The matrix element *ωpq*,*p*<sup>0</sup> *<sup>q</sup>*<sup>0</sup> is function as the Hartree exchange-correlation core, itself-independent of the frequency as:

$$\rho\_{pq,p'q'}\left[\mathcal{f}\_{H\text{xc}}^{\text{ada}}\right] = \int\_{r}\int\_{r'} \rho\_p^\*\left(\vec{r'}\right)\rho\_{p'}^\*\left(\vec{r'}\right)\mathcal{f}\_{H\text{xc}}^{\text{ada}}\left(\vec{r},\vec{r'}\right)\rho\_q\left(\vec{r'}\right)\rho\_{q'}\left(\vec{r'}\right)d\vec{r}d\vec{r'}\tag{65}$$

Where *Xpq* et *Ωpq* are respectively the vectors and eigenvalues.

#### **8. Conclusion**

Our main objective is the determination of the electronic structure, so we have proposed the model of the Functional Density Theory dependent and independent of time which is widely used; it has the advantage of taking into account electron correlation directly in its formalism. To test the relevance of these theoretical calculations and the limit of validity of all these approximations used during the development of this theoretical model, a computer simulation model is necessary. Indeed, it will allow us to obtain the eigenvalues and eigenvectors. They provide enough information to access to the electronic structure and total energy. This theoretical investigation may allow to work toward the realistic modeling of the electronic structure of spherical and cylindrical nanoparticles, to provide effective potentials and orbitals that we can employ to calculate photoemission spectra, and may allow to improve the modeling of the ground-state electronic structure of metal surfaces within the framework of the streaked photoemission calculations, and to develop a theoretical framework for dielectric response of surfaces and nanostructures to external IR fields in order to model plasmonic electric fields [27–30] enhancement near nanostructures systems.

#### **Acknowledgements**

We are grateful to Professor Uwe Thumm who hosted us for 3 months in his James R. Macdonald laboratory at the Kansas State University in the USA, and who offered us an opportunity to collaborate on this subject, as part of the Fulbright Grant Merit Award.

### **Author details**

Brahim Ait Hammou<sup>1</sup> , Abdelhamid El Kaaouachi<sup>1</sup> \*, El Hassan Mounir<sup>2</sup> , Hamza Mabchour<sup>2</sup> , Abdellatif El Oujdi<sup>2</sup> , Adil Echchelh<sup>2</sup> , Said Dlimi<sup>3</sup> and Driss Ennajih<sup>2</sup>

1 Faculty of Sciences of Agadir, Materials and Physicochemistry of the Atmosphere and Climate Group, Agadir, Morocco

2 Faculty of Sciences Ibn Tofail, Laboratory of Energetic Engineering and Materials, Kenitra, Morocco

3 Faculty of Sciences, Laboratory of Electronics, Instrumentation and Energetic, Physics Department, Chouaib Doukkali, El Jadida, Morocco

\*Address all correspondence to: kaaouachi21@yahoo.fr

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*DFT and TDDFT Calculations of Ground and Excited States of Photoelectron Emission DOI: http://dx.doi.org/10.5772/intechopen.111611*

#### **References**

[1] Ambrosio MJ, Thumm U. Energyresolved attosecond interferometric photoemission from Ag(111) and Au (111) surfaces. Physical Review A. 2018; **97**:043431. DOI: 10.1103/PhysRevA.97. 043431

[2] Ambrosio MJ, Thumm U. Electronic structure effects in spatiotemporally resolved photoemission interferograms of copper surfaces. Physical Review A. 2017;**96**:051403(R)

[3] Leone SR, McCurdy CW, Burgdörfer J, Cederbaum LS, Chang Z, Dudovich N, et al. What will it take to observe processes in 'real time'? Nat. Photon. 2014;**8**:162-166. DOI: 10.1038/ nphoton.2014.48

[4] Calegari F, Sansone G, Stagira S, Vozzi C, Nisoli M. Advances in attosecond science. J. Phys. B. 2016;**49**: 062001. DOI: 10.1088/0953-4075/49/6/ 062001

[5] Goulielmakis E, Schultze M, Hofstetter M, Yakovlev VS, Gagnon J, Uiberacker M, et al. Single-cycle nonlinear optics. Science. 2008;**320**: 1614-1617. DOI: 10.1126/science.1157846

[6] Ambrosio MJ, Thumm U. Comparative time-resolved photoemission from Cu(100) and Cu (111) surfaces. Surfaces Phys. Rev. A. 2016;**94**:063424

[7] Tanuma S, Powell CJ, Penn DR. Calculations of electron inelastic mean free paths. IX. Data for 41 elemental solids over the 50 eV to 30 keV range. Surface and Interface Analysis. 2011;**43**: 689-713. DOI: 10.1002/sia.3522

[8] Roth F, Lupulescu C, Darlatt E, Gottwald A, Eberhardt W. Angle resolved photoemission from Cu single crystals: Known facts and a few surprises about the photoemission process. J. Electron Spectrosc. 2016;**208**:2-10. DOI: 10.1016/j.elspec.2015.09.006

[9] Hohenberg P, Kohn W. Inhomogeneous electron gas. Physical Review B. 1969;**136**:864. DOI: 10.1103/ PhysRev.136.B864

[10] Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Physical Review A. 1965;**140**:1133. DOI: 10.1103/PhysRev. 140.A1133

[11] Sahni V, Bohnen KP, Harbola MK. Analysis of the local-density approximation of density-functional theory. Physical Review A. 1988;**37**:1895

[12] Ziesche P, Kurth S, Perdew JP. Density functionals from LDA to GGA. Journal of Computational Materials Science. 1998;**11**(2):122-127

[13] Womack JC, Mardirossian N, Head-Gordon M. Self-consistent implementation of meta-GGA functionals for the ONETEP linear-scaling electronic structure package. The Journal of Chemical Physics. 2016;**145**(20):204114

[14] Koller D, Blaha P, Wien TU, Tran F. Hybrid functionals for solids with an optimized Hartree-Fock mixing parameter. Journal of Physics Condensed Matter. 2013;**25**(43):435503. DOI: 10.1088/0953-8984/25/43/435503

[15] Wigner E. On the interaction of electrons in metals. Physics Review. 1934;**46**:1001. DOI: 10.1103/ PhysRev.46.1002

[16] Perdew JP, Wang Y. Accurate and simple analytic representation of the

electron-gas correlation energy. Physical Review B. 1992;**45**:13244. DOI: 10.1103/ PhysRevB.45.13244

[17] Hedin L, Lundqvist B. Explicit local exchange-correlation potentials. Journal of Physics C. 1971;**4**:2064. DOI: 10.1088/ 0022-3719/4/14/022

[18] Ceperley DM, Alder BJ. Ground state of the electron gas by a stochastic method. Physical Review Letters. 1980; **45**:566. DOI: 10.1103/PhysRevLett. 45.566

[19] Loos PF, Gill PMW. Leading-order behavior of the correlation energy in the uniform electron gas. International Journal of Quantum Chemistry. 2012; **112**:1712

[20] Moller C, Plesset M. Note on an approximation treatment for manyelectron systems. Physics Review. 1934; **46**:618. DOI: 10.1103/PhysRev.46.618

[21] Handler GS. On the removal of the exchange singularity in extended systems. International Journal of Quantum Chemistry. 1988;**33**:173

[22] Zunger A, Ihm J, Cohen ML. Momentum-space formalism for the total energy of solids. Journal of Physics C. 1979;**12**(21):4409. DOI: 10.1088/ 0022-3719/12/21/009

[23] Runge E, Gross EKU. Densityfunctional theory for time-dependent systems. Physical Review Letters. 1984; **52**:997. DOI: 10.1103/PhysRevLett.52.997

[24] Schirmer J. Reexamination of the Runge-Gross action-integral functional. Physical Review A. 2012;**86**:012514. DOI: 10.1103/PhysRevA.86.012514

[25] Mermin ND. Contribution à la modélisation électro-thermique de la cellule de commutation MOSFET-Diode. Physique des solides. France: Universitaire, Institut des Sciences Appliquées de Lyon; 2003

[26] Payne MC, Teter MP, Joannopoulos JD. Iterative minimization techniques for *ab initio* total-energy calculations: Molecular dynamics and conjugate gradients. Reviews of Modern Physics. 1992;**64**(4):1045. DOI: 10.1103/ RevModPhys.64.1045

[27] Hendrik JM, James D. Special points for Brillouin-zone integrations. Pack Physical Review B. 1976;**13**:5188. DOI: 10.1103/PhysRevB.13.5188

[28] Chevreuil P-A, Brunner F, Thumm U, Keller U, Gallmann L. Breakdown of the single-collision condition for soft xray high harmonic generation in noble gases. Optica. 2022. DOI: 10.1364/ OPTICA.471084

[29] Saydanzad E, Li J, Thumm U. Strong-field ionization of plasmonic nanoparticles. Physical Review A. 2022. DOI: 10.1103/PhysRevA.106.033103

[30] Thumm U. Enhanced extreme ultraviolet high-harmonic generation from chromium- doped magnesium oxide. Applied Physics Letters. 2021. DOI: 10.1063/5.0047421

#### **Chapter 3**

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

*Masami Kusunoki*

#### **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 eternallyunknown 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 quasiparticle 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 watersplitting 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.

**Keywords:** spin density functional theory (SDFT), LCAO-natural molecular orbitals (NMO), principal exchange-correlation energy, Heisenberg spin Hamiltonian, secondary correlation energy, superconductivity

#### **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–6], must be stated as incomplete, during a number of active works motivated on it (e.g [7–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 -*KXC*. 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 quasiparticle's wave functions, the magnetic property such as the mean isotropic spinexchange coupling constants {*Ji,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 S1 state Mn4CaO5 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 Mn4 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 Mn4CaOx cluster being thermally distributed in some isomeric substates of any Kok-Si state.

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–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


*atomic spin densities [14].)*

#### *Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

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 (*Nα*, *Nβ*Þ of a pair of the expected values of their operators, (*N*b*α*, *N*b*α*Þ, 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 *θ*/kB (kB 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\Psi = E\Psi(\mathfrak{x}\_1, \mathfrak{x}\_2, \dots, \mathfrak{x}\_{\mathfrak{A}}); M \bowtie N,\tag{1}$$

where *H* is the Hamiltonian operator given by

$$H = T + V\_{\text{ref}} + V\_{\text{ee}},\tag{2}$$

$$=\sum\_{i=1}^{M}\left(-\frac{1}{2}\nabla\_i^2\right) + \sum\_{i=1}^{M}\nu(r\_i,\varepsilon) + \sum\_{i$$

Here, *x<sup>i</sup>* � (*ri*, *si*) represents the (orbital, spin) coordinates of the *i*th electron,*T* is the electron kinetic energy operator; *Vne* 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; *Vee* is the electron Coulomb interaction operator; and *ri,j* � |*ri*-*rj*|. 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 Ψ(*x*1,

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

*x*2, … , *xM*) 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), **<sup>G</sup>**A� [*g<sup>a</sup> <sup>l</sup>* (*r*)] (*a*, the atomic order number from 1 to *NA*; *l*, a set of plural AO quantum numbers). The maximum size of (*Mα*, *Mβ*)-dimensional Hilbert space of the NMO set **M**<sup>0</sup> 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}' = {}^a \mathbf{M}' + {}^\beta \mathbf{M}' \equiv \{ \Psi\_\tau'(\mathbf{x}), \epsilon\_\tau'; \tau = \mathbf{1}\_a, \dots, M\_a, \mathbf{1}\_\beta, \dots, M\_\beta \},\tag{4}$$

$$\sum\_{m\sigma}^{\prime \mathbf{M}^{\prime}} |\Psi^{\prime}\_{m\sigma}\rangle\langle\Psi^{\prime}\_{m\sigma}| \cong \mathbf{1} \text{ for } M\_{\sigma} > > N\_{\sigma}; \sigma = \mathbf{a}, \boldsymbol{\beta}. \tag{5}$$

$$\Psi'\_{m\sigma}(\mathfrak{x}) \equiv \Phi'\_{m\sigma}(\mathfrak{r}) \xi'\_{m\sigma}(\mathfrak{s}),\tag{6}$$

$$
\langle \Psi'\_{m\sigma} | \Psi'\_{m\sigma'} \rangle = \langle \Phi'\_{m\sigma} | \Phi'\_{m\sigma'} \rangle \langle \xi'\_{m\sigma} | \xi'\_{m\sigma'} \rangle = \delta\_{m\sigma, m\sigma'} \delta\_{\sigma, \sigma'}, \tag{7}
$$

$$\Phi\_{m\sigma}(\mathbf{r}) = \sum\_{a,l}^{\mathbf{G}\_{\mathbf{A}}} \mathbf{C}\_{a,l}^{m\sigma} \mathbf{g}\_{l}^{a}(\mathbf{r}) = \sum\_{a=1}^{N\_{A}} \sum\_{l} \mathbf{C}\_{a,l}^{m\sigma} \mathbf{g}\_{l}^{a}(\mathbf{r}). \tag{8}$$

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<sup>a</sup> <sup>l</sup>* <sup>j</sup>*g<sup>a</sup> <sup>l</sup>*0⟩ ¼ *δ<sup>l</sup>*,*l*<sup>0</sup> 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

$$
\hat{\Gamma}[p',\Psi'] \equiv \sum\_{N'}^{\approx} \sum\_{\tau}^{\mathbf{M}'} P\_{N',\tau}' |\Psi'\_{N',\tau}\rangle \langle \Psi'\_{N',\tau}|,\tag{9}
$$

$$=\sum\_{\tau}^{\mathbf{M}'}\sum\_{n\_{\tau}=0}^{1}p'\_{\tau,n\_{\tau}}\lfloor\Psi'\_{\tau,n\_{\tau}}\rfloor\langle\Psi'\_{\tau,n\_{\tau}}\rfloor,\tag{10}$$

with respect to a set of the distribution probabilities, designated *P*0 *N*0 ,*<sup>τ</sup>*; *<sup>N</sup>*<sup>0</sup> <sup>¼</sup> 0, … , <sup>∞</sup>; *<sup>τ</sup>* <sup>¼</sup> 1, … , *<sup>M</sup>* � � for fermions and bosons in Eq. (9) or *p*0 *<sup>τ</sup>*,*n<sup>τ</sup>* f g ; *τ* ¼ 1, … , M; *n<sup>τ</sup>* ¼ 0, 1 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):

$$\operatorname{Tr}\left[\widehat{\Gamma}\right] = \sum\_{\tau}^{\mathbf{M}'} \sum\_{\mathfrak{n}\_\tau=\mathbf{0}}^{1} p'\_{\tau,\mathfrak{n}\_\tau} = \mathbf{1},\tag{11}$$

*Density Functional Theory – New Perspectives and Applications*

$$\hat{N}\_{\sigma} = \sum\_{\mathfrak{r}(\sigma)}^{\prime \mathbf{M}^{\prime}} \hat{a}\_{\mathfrak{r}(\sigma)}^{\dagger} \hat{a}\_{\mathfrak{r}(\sigma)}; \sigma = \mathfrak{a}, \mathfrak{b} \tag{12}$$

$$\mathrm{Tr}\mathrm{Tr}\left[\hat{\mathcal{N}}\_{\sigma}\hat{\Gamma}\right] = \sum\_{\mathfrak{r}(\sigma)}^{\mathsf{''M'}} \sum\_{n\_{\mathfrak{r}(\sigma)}=0}^{1} n\_{\mathfrak{r}(\sigma)} p'\_{\mathfrak{r}(\sigma), n\_{\mathfrak{r}(\sigma)}} \equiv \overline{\mathcal{N}}\_{\sigma}[p', \Psi'], \tag{13}$$

where <sup>b</sup>*a<sup>τ</sup>* and <sup>b</sup>*a*† *<sup>τ</sup>* are the annihilation and creation operators of an electron in the *τ*th NMO-eigenstate, respectively. Subsequently, using the GC entropy

$$\Sigma[p',\Psi'] = -k\_B \text{Tr}\left(\widehat{\Gamma}\,\ln\widehat{\Gamma}\right) = -k\_B \sum\_{\mathfrak{r}}^{\mathbf{M}'} \sum\_{n\_\mathfrak{r}=0}^{1} p'\_{\mathfrak{r},n\_\mathfrak{r}} \ln p'\_{\mathfrak{r},n\_\mathfrak{r}},\tag{14}$$

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

$$
\hat{H} = \hat{H}\_0 + \hat{H}\_1 = \sum\_{\tau}^{\mathbf{M}'} \epsilon'\_{\tau} \hat{a}\_{\tau}^{\dagger} \hat{a}\_{\tau} + \hat{H}\_1,\tag{15}
$$

which in general consists of the principal part *H*b <sup>0</sup> and the secondary part *H*b1, 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:

$$\Omega[p',\Psi'] \equiv \text{Tr}\left[\hat{\Gamma}\left(\theta\ln\hat{\Gamma} + \hat{H}\_0 - \mu\_a\hat{N}\_a - \mu\_\beta\hat{N}\_\beta\right)\right],\tag{16}$$

$$=\sum\_{\tau}^{\mathbf{M}'} \sum\_{n\_{\tau}=0}^{1} p'\_{\tau,n\_{\tau}} \left[ \theta \ln p'\_{\tau,n\_{\tau}} + \langle \Psi'\_{\tau}, n\_{\tau} | \left( \boldsymbol{\varepsilon}'\_{\tau} - \boldsymbol{\mu}\_{\sigma(\tau)} \right) n\_{\tau} | \Psi'\_{\tau}, n\_{\tau} \rangle \right]. \tag{17}$$

In Eqs. (16) and (17) should be noted that *H is replaced by the principal part* b *<sup>H</sup>*<sup>b</sup> <sup>0</sup> <sup>¼</sup> <sup>P</sup> **M**0 *τ ϵ*0 *τ*b*a*† *<sup>τ</sup>*b*aτ:* The variational equations subject to the normalization condition Eq. (11) are given by

$$\frac{\partial}{\partial p'\_{\tau, n\_{\tau}}} \left( \boldsymbol{\Omega} [\{ p'\_{\tau, n\_{\tau}}, \Psi'\_{\tau} \}] + \lambda \left[ \sum\_{\tau}^{\mathbf{M}} \sum\_{n\_{\tau} = 0}^{1} p'\_{\tau, n\_{\tau}} \cdot \mathbf{1} \right] \right) = \mathbf{0},\tag{18}$$

with λ being a Lagrange's multiplier, satisfying

$$e^{-1-\lambda/\theta} = \sum\_{\tau}^{\mathbf{M}'} \sum\_{n\_{\tau}=0}^{1} \langle \Psi'\_{\tau}, n\_{\tau} | \epsilon'\_{\tau} - \mu\_a \overline{N}\_a - \mu\_\beta \overline{N}\_\beta | \Psi'\_{\tau}, n\_{\tau} \rangle. \tag{19}$$

Thus, we obtain an intermediate solution with fixed {Ψ<sup>0</sup> τ}:

$$\begin{aligned} p'\_{\mathsf{r},n\_{\mathsf{r}}} \stackrel{jelds}{\to} p^{0}\_{\mathsf{r},n\_{\mathsf{r}}}(\mathsf{Y}') = \frac{\exp\left[ \langle \Psi'\_{\mathsf{r}}, n\_{\mathsf{r}} | \left( \mu\_{\sigma(\mathsf{r})} - c'\_{\mathsf{r}} \right) n\_{\mathsf{r}} | \Psi'\_{\mathsf{r}}, n\_{\mathsf{r}} \rangle / \theta \right]}{\sum\_{n\_{\mathsf{r}}=0}^{1} \exp\left[ \langle \Psi'\_{\mathsf{r}}, n\_{\mathsf{r}} | \left( \mu\_{\sigma(\mathsf{r})} - c'\_{\mathsf{r}} \right) n\_{\mathsf{r}} | \Psi'\_{\mathsf{r}}, n\_{\mathsf{r}} \rangle / \theta \right]} \tag{20} \end{aligned} \tag{20}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

$$
\widehat{\Gamma}^{0}(\Psi') = \sum\_{\tau}^{\mathbf{M}'} \sum\_{n\_{\tau}=0}^{1} p\_{\tau,n\_{\tau}}^{0}(\Psi') \lfloor \Psi'\_{\tau,n\_{\tau}} \rangle \langle \Psi'\_{\tau,n\_{\tau}}| \,. \tag{21}
$$

Then, the GC potential decreases by the non-negative quantity (the equality appears when *p*<sup>0</sup> = *p*<sup>0</sup> ), namely

$$\begin{split} \Omega[\{p'\_{\tau,n\_{\tau}},\Psi'\_{\tau}\}] &- \Omega\left[\left\{p^{0}\_{\tau,n\_{\tau}},\Psi'\_{\tau}\right\}\right] \\ = \sum\_{\tau}^{\mathbf{M}} \sum\_{n\_{\tau}=0}^{1} p'\_{\tau,n\_{\tau}} \left(\theta \ln p'\_{\tau,n\_{\tau}} + \langle \Psi'\_{\tau}, n\_{\tau} \big| \left(\boldsymbol{\epsilon}'\_{\tau} - \boldsymbol{\mu}\_{\sigma(\boldsymbol{\tau})}\right) n\_{\tau} | \Psi'\_{\tau}, n\_{\tau} \rangle \right) \\ &+ \theta \ln \sum\_{\tau}^{\mathbf{M}} \sum\_{n\_{\tau}=0}^{1} \exp\left[ \langle \Psi'\_{\tau}, n\_{\tau} \big| \left(\boldsymbol{\mu}\_{\sigma(\boldsymbol{\tau})} - \boldsymbol{\epsilon}'\_{\tau} \right) n\_{\tau} | \Psi'\_{\tau}, n\_{\tau} \rangle / \theta \right] \\ = \theta \sum\_{\tau}^{\mathbf{M}} \sum\_{n\_{\tau}=0}^{1} p'\_{\tau,n\_{\tau}} \left[ \ln p'\_{\tau,n\_{\tau}} - \ln p\_{\tau,n\_{\tau}}^{0} (\Psi') \right] \ge 0, \end{split} \tag{22}$$

which is notably induced only by the increase of the GC entropy Σ [*p*<sup>0</sup> ,Ψ<sup>0</sup> ] -Σ [*p*<sup>0</sup> ,Ψ<sup>0</sup> ]. This initial-guess state of the GC ensemble **M**<sup>0</sup> = {Ψ<sup>0</sup> <sup>τ</sup>, *ϵ*<sup>0</sup> *<sup>τ</sup>*} 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} = {}^{a}\mathbf{M} + {}^{\beta}\mathbf{M} \equiv \left\{ \Psi\_{m\sigma}(\mathbf{x}), \varepsilon\_{m\sigma} f\_{m\sigma}(\mu\_{\sigma}); m\sigma \equiv \tau = \mathbf{1}\_{a}, \dots, \mathbf{M}\_{a}, \mathbf{1}\_{\emptyset}, \dots, \mathbf{M}\_{\emptyset} \right\},\tag{23}$$

$$\sum\_{m\sigma}^{\mathsf{TM}} |\Psi\_{m\sigma}\rangle\langle\Psi\_{m\sigma}| \cong \mathbf{1} \text{ for } \mathsf{M}\_{\sigma} > > \mathsf{N}\_{\sigma}; \sigma = \alpha, \beta. \tag{24}$$

$$\Psi\_{m\sigma}(\mathfrak{x}) \equiv \Phi\_{m\sigma}(\mathfrak{r}) \xi\_{m\sigma}(\mathfrak{s}),\tag{25}$$

$$
\langle \Psi\_{m\sigma} | \Psi\_{m\sigma'} \rangle = \langle \Phi\_{m\sigma} | \Phi\_{m\sigma'} \rangle \langle \xi\_{m\sigma} | \xi\_{m\sigma'} \rangle = \delta\_{m\sigma, m\sigma'} \delta\_{\sigma, \sigma'}, \tag{26}
$$

$$f\_{m\sigma}(\mu\_{\sigma}) = p\_{m\sigma, 1} = \left\{ \exp[(\epsilon\_{m\sigma} - \mu\_{\sigma})/\theta] + 1 \right\}^{-1} \equiv f(\epsilon\_{m\sigma} - \mu\_{\sigma}),\tag{27}$$

where it should be noted that only the populated distribution probability *p*<sup>0</sup> *<sup>m</sup>σ*,1 of Eq. (20) converges to the Fermi-Dirac distribution probability *f <sup>m</sup><sup>σ</sup> μσ* ð Þ 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*<sup>0</sup> ) and Γ2(*x*1*x*2, *x*<sup>0</sup> 1*x*0 2), respectively, for the GC ensemble **M** in thermal equilibrium state, as given by

$$\Gamma(\mathfrak{x}, \mathfrak{x}') = \sum\_{\mathfrak{r}}^{\mathbf{M}} f\_{\mathfrak{r}} |\Psi\_{\mathfrak{r}}(\mathfrak{x})\rangle\langle\Psi\_{\mathfrak{r}}(\mathfrak{x}')|,\tag{28}$$

$$
\Gamma\_2(\mathfrak{x}\_1 \mathfrak{x}\_2, \mathfrak{x}\_1' \mathfrak{x}\_2') = \frac{1}{2} [\Gamma(\mathfrak{x}\_1, \mathfrak{x}\_1')\Gamma(\mathfrak{x}\_2, \mathfrak{x}\_2') - \Gamma(\mathfrak{x}\_1, \mathfrak{x}\_2')\Gamma(\mathfrak{x}\_2, \mathfrak{x}\_1')],\tag{29}
$$

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*(*ϵmσ*�*μσ*) 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 NMOfermions for σ = (α, β) in Eq. (31):

$$\mathbf{G} = \mu\_a \overline{N}\_a + \mu\_\beta \overline{N}\_\beta,\tag{30}$$

$$\overline{N}\_{\sigma} = \sum\_{m\sigma}^{\text{"\mkern-1.2mu\sigma}} f(\epsilon\_{m\sigma} - \mu\_{\sigma}).\tag{31}$$

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

*ρσ <sup>r</sup>*, *<sup>r</sup>*<sup>0</sup> ð Þ¼ Tr*<sup>s</sup>* <sup>Γ</sup>*<sup>σ</sup> <sup>x</sup>*, *<sup>x</sup>*<sup>0</sup> ½ � ð Þ *<sup>s</sup>*¼*s*<sup>0</sup> <sup>¼</sup> <sup>X</sup>*<sup>σ</sup>***<sup>M</sup>** *mσ <sup>f</sup> <sup>m</sup><sup>σ</sup> μσ* ð ÞΦ*m<sup>σ</sup>*ð Þ*<sup>r</sup>* <sup>Φ</sup><sup>∗</sup> *<sup>m</sup><sup>σ</sup> r*<sup>0</sup> ð Þ; σ ¼ α, β, (32)

$$
\rho(\mathbf{r}, \mathbf{r}') = \rho\_a(\mathbf{r}, \mathbf{r}') + \rho\_\beta(\mathbf{r}, \mathbf{r}'), \tag{33}
$$

$$
\rho\_{\sigma}(\mathbf{r}) \equiv \rho\_{\sigma}(\mathbf{r}, \mathbf{r}),\\\rho(\mathbf{r}) = \rho\_{a}(\mathbf{r}) + \rho\_{\beta}(\mathbf{r}), \tag{34}
$$

where Tr*<sup>s</sup>* represents the trace on the spin coordinate *s*.

Next, we will prove that the internal energy *U* Γα, Γβ � � as a function of the reduced density matrix Γ = (Γα, Γ*β*) can be decoupled into two parts, as seen in Eq. (36): (1) the principal part *U*<sup>0</sup> *UHFD* Γ*α*, Γ*<sup>β</sup>* � � of Eq. (37) including the principal exchange-correlation energy functional -*KXC* [*ρα*, *ρ*β] defined by Eq. (40) and (2) a secondary part containing only a spin-dependent correlation energy functional Δ*Ecorr UHFD* Γ*α*, Γ*<sup>β</sup>* � � defined by Eq. (43). This ultimate decoupling scheme neglecting the secondary correlation term Δ*Ecorr UHFD* ½Γα, Γβ] of Eq. (43) will be tentatively called "the *Unrestricted Hartree-Fock-Dirac (UHFD) approximation*". because the Dirac's spin permutation operator *(σ<sup>i</sup>* is called spinor)

$$P\_{12}^{\sigma} = \frac{1}{2} [1 + (\sigma\_1, \sigma\_2)] = \frac{1}{2} [1 + 4(\mathfrak{s}\_1, \mathfrak{s}\_2)],\tag{35}$$

including an inner product of two of Pauli's spin operators ð Þ¼ σ1, σ<sup>2</sup> ð Þ 2*s*1, 2*s*<sup>2</sup> , 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[\Gamma\_a,\Gamma\_\beta\right] = U^0\_{UHFD}\left[\Gamma\_a,\Gamma\_\beta\right] + \Delta E^{corr}\_{UHFD}\left[\Gamma\_a,\Gamma\_\beta\right],\tag{36}$$

$$\mathcal{U}\_{\text{UHFD}}^{0}\left[\Gamma\_{a},\Gamma\_{\beta}\right] \equiv \mathbf{T}[\rho] + \mathbf{V}\_{n\epsilon}[\rho] + f[\rho] - \frac{\mathbf{1}}{2}\mathbf{K}\_{\text{XC}}[\rho\_{a},\rho\_{\beta}] + H\_{\text{ES}}[\Gamma\_{a},\Gamma\_{\beta}],\tag{37}$$

$$V\_{n\epsilon}[\rho] = \int d\mathbf{r} \upsilon(\mathbf{r}, \varepsilon)\rho(\mathbf{r}),\tag{38}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

$$J[\rho] = \frac{1}{2} \iiint \frac{dr\_1 dr\_2}{r\_{12}} \rho(r\_1)\rho(r\_2),\tag{39}$$

$$-K\_{\rm XC} \left[ \rho\_a, \rho\_\beta \right] = -\frac{1}{2} \left[ \right] \frac{dr\_1 dr\_2}{r\_{12}} \rho(r\_1, r\_2) \rho(r\_2, r\_1),\tag{40}$$

$$H\_{\rm ES} \left[ \Gamma\_a, \Gamma\_\beta \right] = -2 \iint \frac{dr\_1 dr\_2}{r\_{12}} Q(r\_1, r\_2) Q(r\_2, r\_1)(s\_1, s\_2),\tag{41}$$

$$Q(r\_1, r\_2) = \rho\_a(r\_1, r\_2) - \rho\_\beta(r\_1, r\_2),\tag{42}$$

$$\Delta E\_{\rm UHFD}^{\rm corr} \left[ \Gamma\_a, \Gamma\_\beta \right] = -2 \iint \frac{d\mathbf{r}\_1 d\mathbf{r}\_2}{r\_{12}} \left[ \rho\_a(\mathbf{r}\_1, \mathbf{r}\_2) \rho\_\beta(\mathbf{r}\_2, \mathbf{r}\_1) + \rho\_\beta(\mathbf{r}\_1, \mathbf{r}\_2) \rho\_a(\mathbf{r}\_2, \mathbf{r}\_1) \right] (\mathbf{s}\_1, \mathbf{s}\_2), \tag{43}$$

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 *HES* Γ*α*, Γ*<sup>β</sup>* � � 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{split} \Gamma(\mathbf{x}\_1, \mathbf{x}'\_2) \Gamma(\mathbf{x}\_2, \mathbf{x}'\_1) &= \sum\_{\sigma, m\sigma}^{\mathbf{M}} f\_{m\sigma}(\mu\_{\sigma}) \Phi\_{m\sigma}(r\_1) \Phi\_{m\sigma}(r'\_2) \,^\* \\ \times \sum\_{\sigma', m\sigma'}^{\mathbf{M}} f\_{m\sigma'}(\mu\_{\sigma'}) \Phi\_{m\sigma'}(r\_2) \Phi\_{m\sigma}(r'\_1) \,^\* \left| \xi\_{m\sigma}(s\_2) \xi\_{m\sigma'}(s\_1) \right\rangle \langle \xi\_{m\sigma}(s'\_1) \xi\_{m\sigma'}(s'\_2) \,| \,. \end{split} \tag{44}$$

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\_{m\sigma}(\mathfrak{s}\_2)\xi\_{m\sigma'}(\mathfrak{s}\_1) = P\_{12}^{\sigma}\xi\_{m\sigma}(\mathfrak{s}\_1)\xi\_{m\sigma'}(\mathfrak{s}\_2). \tag{45}
$$

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

$$\begin{aligned} \Gamma(\mathbf{x}\_1, \mathbf{x}\_2') \Gamma(\mathbf{x}\_2, \mathbf{x}\_1') \\ = \rho(\mathbf{r}\_1, \mathbf{r}\_2') \rho(\mathbf{r}\_2, \mathbf{r}\_1') \left[ \frac{1}{2} (1 + \sigma\_{1\varepsilon}\sigma\_{2\varepsilon}) + \frac{1}{4} (\sigma\_{1\varepsilon}\sigma\_{2\cdots} + \sigma\_{1\cdots 2\cdots}) \right], \end{aligned} \tag{46}$$

$$= \frac{1}{2} \left[ \rho\_a(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_a(\mathbf{r}\_2, \mathbf{r}\_1') + \rho\_\beta(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_\beta(\mathbf{r}\_2, \mathbf{r}\_1') \right]$$

$$\quad + \frac{1}{2} \left[ \rho\_a(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_\beta(\mathbf{r}\_2, \mathbf{r}\_1') + \rho\_\beta(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_a(\mathbf{r}\_2, \mathbf{r}\_1') \right]$$

$$\quad + \frac{1}{2} Q(\mathbf{r}\_1, \mathbf{r}\_2') Q(\mathbf{r}\_2, \mathbf{r}\_1') (\sigma\_1, \sigma\_2)$$

$$\quad + \frac{1}{2} \left[ \rho\_a(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_\beta(\mathbf{r}\_2, \mathbf{r}\_1') + \rho\_a(\mathbf{r}\_1, \mathbf{r}\_2') \rho\_\beta(\mathbf{r}\_2, \mathbf{r}\_1') \right] (\sigma\_1, \sigma\_2),$$

with the use of the spin density matrix *Q* of Eq. (42) and the off-diagonal spinors *σ<sup>j</sup>*� of Eq. (48):

$$
\sigma\_{\rm j\pm} = \sigma\_{\rm j\varepsilon} \pm i \sigma\_{\rm j\gamma}, j = \mathbf{1}, 2. \tag{48}
$$

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* Γ*α*, Γ*<sup>β</sup>* � � in 3D-spin space decomposed into the principal part *U*<sup>0</sup> *UHFD* Γ*α*, Γ*<sup>β</sup>* � � and the secondary part Δ*EUHFDcorr* in Eq. (36). (QED)

Since all the DFT calculations can be made in the binary-spin Hilbert space, we must take the trace of *HE*<sup>S</sup> [Γα, Γβ] 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\_{\rm ES}[\rho\_a, \rho\_\beta] \equiv \operatorname{Tr}\_{\rm S} \{ H\_{\rm ES} \left[ \Gamma\_a, \Gamma\_\beta \right] \} = -\frac{1}{2} K\_{\rm XC} [\rho\_a, \rho\_\beta] \,. \tag{49}$$

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

$$\mathcal{U}^{0}\_{\text{UHFD}}\left[\rho\_{a},\rho\_{\beta}\right] \equiv \text{Trs}\left(\mathcal{U}^{0}\_{\text{UHFD}}\left[\Gamma\_{a},\Gamma\_{\beta}\right]\right) = \text{T}[\rho] + V\_{\text{ref}}[\rho] + f[\rho] - K \text{xc}\left[\rho\_{a},\rho\_{\beta}\right],\tag{50}$$

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

$$
\hat{H}\_{\text{UHFD}}^{0} = \sum\_{\mathbf{r}}^{\mathbf{M}} \epsilon\_{\mathbf{r}} \hat{a}\_{\mathbf{r}}^{\dagger} \hat{a}\_{\mathbf{r}},\tag{51}
$$

leading to

$$\mathbb{E}\left[U\_{\text{UHFD}}^{0}\left[\rho\_{a},\rho\_{\beta}\right]\right] \equiv \text{Tr}\_{\mathbf{M}}\left[\Gamma\widehat{H}\_{0}\right] = \int d\mathbf{r} \sum\_{\tau}^{\mathbf{M}} \epsilon\_{\mathbf{f}} f\_{\tau} \left(\epsilon\_{\tau} - \mu\_{\sigma(\mathbf{t})}\right) \Phi\_{\tau}^{\*}\left(\mathbf{r}\right) \Phi\_{\tau}(\mathbf{r}).\tag{52}$$

Similarly, *U*<sup>0</sup> *UHFD* can be expanded as the GC ensemble average of a self-consistent effective Hamiltonian as given by

$$\begin{split} \mathcal{U}\_{\text{UHFD}}^{0}\left[\rho\_{\alpha},\rho\_{\beta}\right] &= \int d\mathbf{r} \left\{ \sum\_{\tau}^{\mathbf{M}} f\_{\tau} \left(\varepsilon\_{\tau} - \mu\_{\sigma(\tau)}\right) \Phi\_{\tau}^{\*}\left(\mathbf{r}\right) \right. \\ &\times \left[ -\frac{1}{2} \nabla^{2} \Phi\_{\tau}(\mathbf{r}) + \int \nu\_{\text{NMO}}(\mathbf{r},\mathbf{r}') \Phi\_{\tau}(\mathbf{r}') d\mathbf{r}' \right] \right\} \end{split} , \tag{53}$$

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

$$\upsilon\_{\rm NMO}(\mathbf{r}, \mathbf{r}') = \left[\upsilon(\mathbf{r}, \varepsilon) + \frac{1}{2} \int\_{\mathbf{M}} d\mathbf{r}'' \frac{\rho(\mathbf{r}'')}{|\mathbf{r} - \mathbf{r}''|} \right] \delta(\mathbf{r} - \mathbf{r}') - \frac{\rho(\mathbf{r}, \mathbf{r}')}{2|\mathbf{r} - \mathbf{r}'|},\tag{54}$$

$$\rho(\boldsymbol{r}) = \sum\_{\tau}^{\mathbf{M}} f\_{\tau} \left( \boldsymbol{e}\_{\tau} - \mu\_{\sigma(\tau)} \right) \left| \Phi\_{\tau}(\boldsymbol{r}) \right|^{2}, \tag{55}$$

$$\rho(\mathbf{r}, \mathbf{r}') = \sum\_{\mathbf{r}}^{\mathbf{M}} f\_{\mathbf{r}} \left( e\_{\mathbf{r}} - \mu\_{\sigma(\mathbf{r})} \right) \Phi\_{\mathbf{r}}^{\*} (\mathbf{r}') \Phi\_{\mathbf{r}}(\mathbf{r}).\tag{56}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

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

$$\begin{split} \boldsymbol{U}\_{\text{UHFD}}^{0} \left[ \rho\_{\alpha}, \rho\_{\beta} \right] - \text{Tr}\_{\mathbf{M}} \left[ \Gamma \hat{H}\_{0} \right] &= \int d\mathbf{r} \left\{ \sum\_{\tau}^{\mathbf{M}} f\_{\tau} \left( \boldsymbol{\varepsilon}\_{\tau} - \mu\_{\sigma(\boldsymbol{\varepsilon})} \right) \Phi\_{\tau}^{\*} (\boldsymbol{r}) \right. \\\\ \left[ -\frac{1}{2} \nabla^{2} \Phi\_{\tau} (\boldsymbol{r}) + \int \boldsymbol{\nu}\_{\text{NMO}} (\boldsymbol{r}, \boldsymbol{r}') \Phi\_{\tau} (\boldsymbol{r}') d\boldsymbol{r}' - \boldsymbol{\varepsilon}\_{\tau} \Phi\_{\tau} (\boldsymbol{r}) \right] &= \mathbf{0} \end{split} \tag{57}$$
 
$$\begin{split} -\frac{1}{2} \nabla^{2} \Phi\_{\tau} (\boldsymbol{r}) + \int d\mathbf{r}' \boldsymbol{\nu}\_{\text{NMO}} (\boldsymbol{r}, \boldsymbol{r}') \Phi\_{\tau} (\boldsymbol{r}') &= \boldsymbol{\varepsilon}\_{\tau} \Phi\_{\tau} (\boldsymbol{r}), \\ \text{and } \mathbf{h}. \mathbf{c}, \text{ for all } \,\boldsymbol{\varepsilon} \in \mathbf{M}. \end{split} \tag{58}$$

This central (not variational!) solution will be stocked as the presumed GC ensemble **M**, in which the eigenvalues f g *ϵτ* are usually assumed to be rearranged from the minimum *ϵ*1*<sup>σ</sup>* to the maximum *ϵ*<sup>M</sup><sup>σ</sup> in the order of increasing energies, first for α-spin NMOs and second for β-spin NMOs, as

$$\pi = \mathbf{1}\_a, \mathbf{2}\_a, \dots, \mathbf{M}\_a, \mathbf{1}\_\beta, \mathbf{2}\_\beta, \dots, \mathbf{M}\_\beta = \mathbf{1}, \mathbf{2}, \dots, \mathbf{M} \left( \mathbf{M} = \mathbf{M}\_a + \mathbf{M}\_\beta \right). \tag{59}$$

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, Δ*Ecorr UHFD* ½Γα, Γβ], 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

*H*b1 *UHFD* <sup>¼</sup> Tr*<sup>s</sup>* <sup>Δ</sup>*Ecorr UHFD* Γ*α*, Γ*<sup>β</sup>* � � � � ¼ � <sup>1</sup> 2 ðð *dr*1*dr*<sup>2</sup> *r*12 Tr*<sup>s</sup> ρα*ð Þ *r*1, *r*<sup>2</sup> *ρβ*ð Þþ *r*2, *r*<sup>1</sup> *ρβ*ð Þ *r*1, *r*<sup>2</sup> *ρα*ð Þ *r*2, *r*<sup>1</sup> �� �ð Þg *<sup>σ</sup>*1, *<sup>σ</sup>*<sup>2</sup> ¼ �X*<sup>α</sup>***<sup>M</sup>** *mα* X*β***M** *mβ f ϵm<sup>α</sup>* � *μα* ð Þ *f ϵm<sup>β</sup>* � *μβ* � � ðð *<sup>d</sup>r*1*dr*<sup>2</sup> *r*12 � <sup>Φ</sup><sup>∗</sup> *<sup>m</sup><sup>α</sup>*ð Þ *<sup>r</sup>*<sup>2</sup> <sup>Φ</sup>*m<sup>α</sup>*ð Þ *<sup>r</sup>*<sup>1</sup> <sup>Φ</sup><sup>∗</sup> *<sup>m</sup><sup>β</sup>*ð Þ *r*<sup>1</sup> Φ*m<sup>β</sup>*ð Þ *r*<sup>2</sup> h i *<sup>σ</sup>*þ,*mα*,*σ*�,*m<sup>β</sup>* <sup>þ</sup> *<sup>σ</sup>*�,*<sup>m</sup>ασ*þ,*m<sup>β</sup>* � � ¼ �X*<sup>α</sup>***<sup>M</sup>** *mα* X*β***M** *mβ σ*þ,*mα*,*σ*�,*m<sup>β</sup>* þ *σ*�,*<sup>m</sup>ασ*þ,*m<sup>β</sup>* � � ðð *<sup>d</sup>r*1*dr*<sup>2</sup> *r*12 � <sup>Φ</sup><sup>∗</sup> *<sup>m</sup><sup>α</sup>*ð Þ *<sup>r</sup>*<sup>2</sup> <sup>Φ</sup><sup>∗</sup> *<sup>m</sup><sup>β</sup>*ð Þ *r*<sup>1</sup> Φ*m<sup>β</sup>*ð Þ *r*<sup>2</sup> Φ*m<sup>α</sup>*ð Þ *r*<sup>1</sup> h ib*a*† *mβ*b*am<sup>β</sup>*b*a*† *<sup>m</sup><sup>σ</sup>*b*amα:* (60)

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 {*Si*, *i* = 1, … ,*n*}, which arise from transition metal (TM) cations, C/N/Oradicals, -C=C- bond radicals, and so on. These spins are quantum-mechanically interacting with each other via the exchange coupling constants (*Ji, j*) in the Heisenberg spin Hamiltonian defined by

$$H\_{\rm ex} = -2\sum\_{i=1}^{n-1} \sum\_{j=i+1}^{n} J\_{i,j}(\mathbf{S}\_i, \mathbf{S}\_j). \tag{61}$$

However, this *Hex* model takes account only the pure spin operators {*Si*} but does not contain any kind of polarized spins of ligand atoms, designated {*siL*}, "Why?" The most fundamentally important question is, "What is the origin of *Hex*?" 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 *HES* [Γα,Γβ] of Eqs. (41) and (42).

Since in the binary spin space appear 2*n*�<sup>1</sup> (*n* ≥ 2) mutually-independent up/down*i*ES arrangements (ESA), which represent one ferromagnetic and the other antiferromagnetic 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}^{(k)}, \mathbf{M}^{(k)} = \,^a \mathbf{M}^{(k)} + \,^\beta \mathbf{M}^{(k)};\tag{62}$$

$$\mathbf{^\sigma M^{(k)}} \equiv \left\{ \Psi^{(k)}\_{m\sigma}(\mathbf{x}), \epsilon^{(k)}\_{m\sigma}, \mathbf{f}^{(k)}\_{m\sigma} \Big(\mathbf{e}^{(k)}\_{m\sigma} - \mu^{(k)}\_{\sigma}\Big); m\sigma = \mathbf{1}, \dots, M\_{\sigma} \right\}; \sigma = \mathbf{a}, \emptyset. \tag{63}$$

Practically, we can calculate only a set of 2n�<sup>1</sup> principal internal energy functionals:

$$\mathcal{U}\_{\text{UHFD}}^{(k)}\left[\rho\_a^{(k)},\rho\_\beta^{(k)}\right] = T^{(k)}[\rho] + \mathcal{V}\_{\text{nc}}^{(k)}[\rho] + J^{(k)}[\rho] - K\_{\text{XC}}^{(k)}\left[\rho\_a^{(k)},\rho\_\beta^{(k)}\right]; k = 1,\ldots,2^{n-1}.\tag{64}$$

However, the origin of *Hex* must be traced to a set of 2*n*�<sup>1</sup> equality relationships between the principal exchange-correlation energy functional, �*K*ð Þ*<sup>k</sup> XC ρ* ð Þ*k <sup>α</sup>* , *ρ* ð Þ*k β* h i, and the projected value of the spin-dependent XC energy functional, i.e.

$$-K\_{\rm XC}^{(k)}\left[\boldsymbol{\rho}\_{a}^{(k)},\boldsymbol{\rho}\_{\boldsymbol{\beta}}^{(k)}\right] = \mathbf{2}H\_{\rm ES}^{(k)}\left[\boldsymbol{\rho}\_{a}^{(k)},\boldsymbol{\rho}\_{\boldsymbol{\beta}}^{(k)}\right] = \mathbf{2}\mathrm{Tr}\_{t}\left(H\_{\rm ES}^{(k)}\left[\boldsymbol{\Gamma}\_{a}^{(k)},\boldsymbol{\Gamma}\_{\boldsymbol{\beta}}^{(k)}\right]\right); k = \mathbf{1},\ldots,\mathbf{2}^{n-1}.\tag{65}$$

We note that the *k*th projected value of *Si*, *S<sup>j</sup>* � � onto the binary Hilbert space spanned by the *k*th GC ensemble **M**(k) must depend not only on the *k*th principal exchange-correlation energy between *i*ES and *j*ES, but also on the polarized spins of bridging and non-bridging ligand atoms, *iLaj* (*j* 6¼ *i*) and *iLnb*, respectively, via the conservation law of each projected spin number of *ni* (k) and *nj* (k), defined by

$$m\_i^{(k)} \equiv \mathfrak{Z} \mathfrak{S}\_i \sigma\_{x,i}^{(k)}; \sigma\_{x,i}^{(k)} = \pm \mathbf{1}, k = \mathbf{1}, \dots, \mathfrak{Z}^{n-1}; i = \mathbf{1}, \dots, n. \tag{66}$$

$$\sum\_{i=1}^{n} \mathbf{S}\_i \sigma\_{x,i}^{(k)} \ge \mathbf{0}, k = \mathbf{1}, \dots, \mathbf{2}^{n-1}. \tag{67}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

Although the arrangement order of *n* ESA-codes *σ*ð Þ*<sup>k</sup> <sup>z</sup>*,*<sup>i</sup>* ¼ 1 or � 1, *i* ¼ 1 � *n* n o leaves the choice of one's best depending on the arrangement order of *n* TM-cation spins {*Si*, �1 – *n*}, a non-negative sum rule of Eq. (67) must be satisfied owing to the timereversal symmetry.

For the *i*th transition metal (TM) cations, its non-bridging ligand atoms *iLnb* and its bridging ligand atoms *iLaj* between the *i*th and *j*th TM-cations, this wave-packet spin projection may be entrusted to the respective spin operator by itself, that is *Si*, *siLnb* and *siLaj*, 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:

$$\mathbf{S}\_{\mathsf{i}\mathbf{j}}\mathbf{g}^{a}\_{l}(\mathfrak{x}) = \delta\_{a,i}\delta\_{l,3d\_{i}}\mathbf{S}\_{\mathsf{i}\mathbf{j}}\mathbf{g}^{i}\_{\mathfrak{j}d}(\mathfrak{x}),\tag{68}$$

$$\mathfrak{s}\_{iLnb}\mathfrak{g}\_{l}^{a}(\mathfrak{x}) = \delta\_{a\in iLnb}\delta\_{l,2p}\mathfrak{s}\_{iLnb}\mathfrak{g}\_{2p}^{iLnb}(\mathfrak{x}),\tag{69}$$

$$s\_{iLaj}g\_l^{a}(\mathbf{x}) = \delta\_{a \in iLaj} \delta\_{l,2p} s\_{iLaj} g\_{2p}^{iLaj}(\mathbf{x}); f\_{iLajm\sigma} = \frac{f\_{im\sigma}}{\nu\_{iLaj}} = \frac{f\_{im\sigma}}{2},\tag{70}$$

where we have introduced the share frequency *νiLajm<sup>σ</sup>* among *n* different subsets (in this case, it is 2), *δ<sup>a</sup>* <sup>∈</sup> *<sup>A</sup>* = (1 for *a* ∈ A; 0 for otherwise) and *δx,*<sup>y</sup> is Kronecker's δ.

Concomitantly, the *k*th GC ensemble **M**(*k*) might be decomposed into *n* spindependent NMO-subsets associated with these elements, and the other spin-free subset as in Eq. (71), each **M**ð Þ*<sup>k</sup> <sup>i</sup>* further decomposed into three components as in Eqs. (72) and (73), finally to define the *i*th ES in terms of two components in Eq. (74).

$$\mathbf{M}^{(k)} = \sum\_{i=1}^{n} \mathbf{M}\_i^{(k)} + \mathbf{M}'^{(k)},\tag{71}$$

$$\mathbf{M}\_i^{(k)} = \mathbf{M}\_{id}^{(k)} + \mathbf{M}\_{iL}^{(k)} + \mathbf{M}\_{io}^{(k)},\tag{72}$$

$$\mathbf{M}\_{iL}^{(k)} = \sum\_{j \neq i}^{n} \mathbf{M}\_{iLaj}^{(k)} + \mathbf{M}\_{iLnb}^{(k)},\tag{73}$$

$$\mathbf{M}\_{i\text{ES}}^{(k)} = \mathbf{M}\_{id}^{(k)} + \sum\_{j \neq i}^{n} \mathbf{M}\_{iLaj}^{(k)},\tag{74}$$

where the *i*th subset in Eq. (72) consists of the subset **M**ð Þ*<sup>k</sup> id* associated with (3d, 4 s, 4p)-electron AO's in the *i*th TM cation, the subset **M**ð Þ*<sup>k</sup> iL* associated with 2p-valence electron AO's in the *iL*th assembly of ligand atoms, and the subset **M**ð Þ*<sup>k</sup> io* associated with other doubly-occupied core-shell AO's in the *i*th TM cation. The *iL*th subset in Eq. (73) can be further decomposed into two kinds of ligand assembly: (1) a *thermal equipartition* half of the *iLaj*th assembly of bridging ligand atoms between the *i*th and *j*th TM cations, which can mediate the antiferromagnetic super-exchange coupling, and (2) the *iLnb*th assembly of non-bridging ligand atoms around the *i*th TM cation, which can be either paramagnetically or diamagnetically polarized depending on the ligand C/N/O atomic structure and hence control the *i*th ES magnitude via *the spin number (ni) conservation law* governing the Mulliken atomic spin densities {*M*ð Þ*<sup>k</sup> <sup>a</sup>* } [14] of these magnetically-interacting atoms, finally to define the *i*th ES density *n*ð Þ*<sup>k</sup> iEF* so as to satisfy Eq. (76):

*Density Functional Theory – New Perspectives and Applications*

*iLnb*0

$$n\_{\rm iEF}^{(k)} \equiv \mathcal{M}\_{\rm id}^{(k)} + \sum\_{j \neq i}^{n} \sum\_{\rm iLaj} \frac{\mathcal{M}\_{\rm iLij}^{(k)} \mathcal{J} \left( \mathcal{M}\_{\rm iLij}^{(k)} \sigma\_{x,i}^{(k)} \right)}{\nu\_{\rm iLij}^{(k)}} + \delta\_{k > 1} \sum\_{j \neq i}^{n} \sum\_{\rm iLaj} \sigma\_{x,i}^{(k)} \frac{\mathcal{M}\_{\rm iLij}^{(1)} \mathcal{J} \left( -\sigma\_{x,j}^{(1)} \sigma\_{x,i}^{(1)} \right)}{\nu\_{\rm iLij}^{(1)}}; \tag{75}$$
 
$$n\_{\rm iLib}^{(k)} \equiv \sum \mathcal{M}\_{\rm iLib}^{(k)}; n\_{\rm iEF}^{(k)} + n\_{\rm iLib}^{(k)} = n\_{\rm i} k = \mathbf{1}(\mathbf{F}), \mathbf{2}(\mathbf{AF}), \dots, \mathbf{2}^{n-1}(\mathbf{AF}). \tag{76}$$

Here, ϑ(sign1\*sign2) is the Heaviside step function, and *ν* ð Þ*k iLaj* is the frequency of the *iLaj* spin density shared by plural ESs with the same sign, which can be calculated as

$$
\omega\_{iLaj}^{(k)} = \mathbf{1} + \mathfrak{G}\left(\sigma\_{x,i}^{(k)}\sigma\_j^{(k)}\right) \text{ for all } k'\mathbf{s}, \tag{77}
$$

(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*n*�<sup>1</sup> } from the spin-dependent XC energy functional in Eq. (65), as follows:

The *k*th spin density matrix for Eq. (42) can be decomposed into

$$Q^{(k)}(r,r') = \sum\_{\sigma=a}^{\beta} (-1)^{\sigma} \sum\_{m\sigma}^{\mathbf{\prime} \mathbf{M}^{(k)}} \Phi\_{m\sigma}(r) \Phi\_{m\sigma}(r)^{\*} ; (-1)^{a/\beta} = \pm 1,\tag{78}$$

$$\hat{\mathbf{y}} = \sum\_{i=1}^{n} \left[ \mathbf{Q}\_{i\text{ES}}^{(k)}(\mathbf{r}, \mathbf{r}') + \mathbf{Q}\_{i\text{Lub}}^{(k)}(\mathbf{r}, \mathbf{r}') + \mathbf{Q}\_{io}^{(k)}(\mathbf{r}, \mathbf{r}') \right] + \mathbf{Q}'^{(k)}(\mathbf{r}, \mathbf{r}'),\tag{79}$$

$$\begin{split} \mathcal{I} &= \sum\_{i=1}^{n} \Big[ \mathcal{Q}\_{iES}^{(k)}(\mathbf{r}, \mathbf{r}') + \mathcal{Q}\_{iLnb}^{(k)}(\mathbf{r}, \mathbf{r}') \Big] \\ &+ \sum\_{i=1}^{n} \Big[ \rho\_{i\alpha}^{(k)}(\mathbf{r}, \mathbf{r}') - \rho\_{i\alpha\beta}^{(k)}(\mathbf{r}, \mathbf{r}') \Big] + \rho\_{\alpha}^{\prime (k)}(\mathbf{r}, \mathbf{r}') - \rho\_{\beta}^{\prime (k)}(\mathbf{r}, \mathbf{r}'), \end{split} \tag{80}$$

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

$$\mathbf{2}H\_{\mathrm{ES1}}^{(k)}\left[\rho\_{\alpha}^{(k)},\rho\_{\beta}^{(k)}\right] \equiv \mathbf{2}E\_{\mathrm{ES1}}^{(k)} = -\sum\_{i=1}^{n} J\_{i,i}^{(k)} \mathbf{S}\_{i} (\mathbf{S}\_{i} + \mathbf{1}),\tag{81}$$

$$J\_{i,i}^{(k)} \equiv \frac{1}{2n\_i^2} \left[ \right] \frac{dr\_1 dr\_2}{r\_{12}} \left[ \rho\_{i \text{ES}}^{(k)}(r\_1, r\_2) + \rho\_{i \text{Lub}}^{(k)}(r\_1, r\_2) \right] \left[ \rho\_{i \text{ES}}^{(k)}(r\_2, r\_1) + \sigma\_{i \text{Lub}}^{(k)}(r\_2, r\_1) \right],\tag{82}$$

and

$$2H\_{\rm ES2}^{(k)}\left[\rho\_a^{(k)},\rho\_\beta^{(k)}\right] = 2E\_{\rm ES2}^{(k)} = \langle H\_{\rm ex}\rangle^{(k)} = -2\sum\_{i$$

$$\left< \left( \mathbf{S}\_i, \mathbf{S}\_j \right) \right>^{(k)} = \operatorname{Tr} \left[ \left( \mathbf{S}\_i, \mathbf{S}\_j \right) \Gamma\_2^{(k)} \left( \mathbf{x}\_1 \mathbf{x}\_2, \mathbf{x}\_1 \mathbf{x}\_2 \right) \right],\tag{84}$$

$$\mathbf{L} = \sum\_{M\_i=-\mathbf{S}\_i}^{\mathbf{S}\_i} \sum\_{M\_j=-\mathbf{S}\_j}^{\mathbf{S}\_j} \text{Tr}\left[ \langle \mathbf{M}\_i, \mathbf{M}\_j | \{ \mathbf{S}\_i, \mathbf{S}\_j \} | \mathbf{M}\_i, \mathbf{M}\_j \rangle \times \langle \mathbf{M}\_i, \mathbf{M}\_j | \Gamma\_{i,j}^{(k)} (\mathbf{x}\_1 \mathbf{x}\_2, \mathbf{x}\_1 \mathbf{x}\_2) | \mathbf{M}\_i, \mathbf{M}\_j \rangle \right], \tag{85}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

with

$$
\Gamma\_{i,j}^{(k)}(\mathfrak{x}\_1\mathfrak{x}\_2, \mathfrak{x}\_1\mathfrak{x}\_2) = \frac{2}{2} \left[ \Gamma\_{i\to 3}^{(k)}(\mathfrak{x}\_1, \mathfrak{x}\_1) \Gamma\_{j\to 3}^{(k)}(\mathfrak{x}\_2, \mathfrak{x}\_2) - \Gamma\_{i\to 3}^{(k)}(\mathfrak{x}\_1, \mathfrak{x}\_2) \Gamma\_{j\to 3}^{(k)}(\mathfrak{x}\_2, \mathfrak{x}\_1) \right],\tag{86}
$$

respectively, and the second and third terms to a major *Hex*-less. component in the principal XC energy, as given by

$$-K\_{\rm XC-H\_{\alpha}}^{(k)}\left[\rho\_{a}^{(k)},\rho\_{\beta}^{(k)}\right] = \mathbf{2}H\_{\rm ES0}^{(k)}\left[\rho\_{\rm ES0,a}^{(k)},\rho\_{\rm ES0,\beta}^{(k)}\right] = \mathbf{2}E\_{\rm ES0}^{(k)},\tag{87}$$

$$\begin{split} \mathcal{I} &= -\frac{1}{2} \Big[ \Big| \frac{dr\_1 dr\_2}{r\_{12}} \Big[ \rho^{(k)}(r\_1, r\_2) - \sum\_{i=1}^{n} \Big\{ \rho^{(k)}\_{i \to S}(r\_1, r\_2) + \rho^{(k)}\_{i \to nb}(r\_1, r\_2) \Big\} \Big] \\ &\times \Big[ \rho^{(k)}(r\_2, r\_1) - \sum\_{i=1}^{n} \Big\{ \rho^{(k)}\_{i \to S}(r\_2, r\_1) + \rho^{(k)}\_{i \to nb}(r\_2, r\_1) \Big\} \Big]; \end{split} \tag{88}$$

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

$$\mathbf{S}\_{\text{tot}} = \sum\_{i=1}^{n} \left( \mathbf{S}\_{i} + \sum\_{i \text{L}aj} s\_{i \text{L}aj} + \sum\_{i \text{L}ub} s\_{i \text{L}ub} \right), \tag{89}$$

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

$$\sum\_{M\_i=-S\_i}^{S\_i} |M\_i\rangle\langle M\_i| = \mathbf{1}\_i; i = \mathbf{1}, \dots, n,\tag{90}$$

and

$$
\langle a\_i | M\_i \rangle = \frac{\delta\_{\mathcal{M}\_i, \mathcal{S}\_i}}{\sqrt{|n\_i|}}, \\
\langle \beta | \mathcal{M}\_i \rangle = \frac{\delta\_{\mathcal{M}\_i, -\mathcal{S}\_i}}{\sqrt{|n\_i|}}; \tag{91}
$$

while the *iLaj*th and *iLnb*th ES spin operators, *siLaj* and *siLnb,* are assumed to automatically respond to the up or down state of *Si*. Then, the expected value of the zcomponent of **S**tot in the *k*th ESA state is given by

$$\langle \mathbf{S}\_{\text{tot},\mathbf{x}} \rangle^{(k)} = \sum\_{i=1}^{n} \text{Tr} \left\{ \mathbf{S}\_{i,\mathbf{x}} \left[ \Gamma\_{i\to \mathbf{S}}^{(k)}(\mathbf{x}, \mathbf{x}) + \Gamma\_{i\to \mathbf{b}}^{(k)}(\mathbf{x}, \mathbf{x}) \right] \right\}, \tag{92}$$

$$\mathcal{S} = \sum\_{i=1}^{n} \mathcal{S}\_{i} \sigma\_{x,i}^{(k)} \left( P\_{i \to \mathcal{S}}^{(k)} + P\_{i \to nb}^{(k)} \right) = \sum\_{i=1}^{n} \mathcal{S}\_{i} \sigma\_{x,i}^{(k)},\tag{93}$$

$$
\Gamma\_{i \to \mathbf{S}}^{(k)}(\mathbf{x}, \mathbf{x}') = \Gamma\_{id}^{(k)}(\mathbf{x}, \mathbf{x}') + \sum\_{j \neq i}^{n} \sum\_{i \text{Laj}'} \Gamma\_{i \text{Laj}'}^{(k)}(\mathbf{x}, \mathbf{x}'), \tag{94}
$$

$$P\_{i \to S}^{(k)} = \frac{n\_{i \to S}^{(k)}}{n\_i} = \frac{1}{n\_i} \int dr Q\_{i \to S}^{(k)}(r, r), \tag{95}$$

$$P\_{iLnb}^{(k)} = \frac{n\_{iLnb}^{(k)}}{n\_i} = \frac{1}{n\_i} \int dr Q\_{iLnb}^{(k)}(r, r),\tag{96}$$

$$P\_{i \to S}^{(k)} + P\_{i \to nb}^{(k)} = \left[ n\_{i \to S}^{(k)} + n\_{i \to nb}^{(k)} \right] / n\_i = \mathbf{1}. \tag{97}$$

It is important to remind that the last additional term in Eq. (75) in the antiferromagnet *kAF*th ESA state was absolutely required for a systematic better agreement with experimental *Ji,j* indicating that a large-positive *M*ð Þ <sup>1</sup>*<sup>F</sup> iLaj* density in the *1F*th-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*ð Þ *kAF iES* and *<sup>n</sup>*ð Þ *kAF jES* in the *kAF*th ESA state, as given by a programmatic equation [13],

$$\mathbf{M}\_{i\mathbf{I}\mathbf{a}\mathbf{j}}^{(k)} \Longleftrightarrow \frac{\mathbf{1}}{2} \sigma\_{x,i}^{(k)} \mathbf{M}\_{i\mathbf{I}\mathbf{a}\mathbf{j}}^{(1)} \boldsymbol{\theta} \left( -\sigma\_{x,i}^{(k)} \sigma\_{x,j}^{(k)} \right) + \mathbf{M}\_{i\mathbf{I}\mathbf{a}\mathbf{j}}^{(k)}.\tag{98}$$

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

$$
\langle \left( \mathbf{S}\_i, \mathbf{S}\_j \right) \rangle^{(\mathbf{k})} = n\_{i \to S}^{(k)} n\_{j \to S}^{(k)} \left( \mathbf{1} - \mathbf{S}\_{i \to S, j \to S}^{(k)} \right), \tag{99}
$$

$$S\_{iES,jES}^{(k)} \equiv \frac{\mathrm{Tr}\left[\Gamma\_{iES}^{(k)}(\mathbf{x}\_1, \mathbf{x}\_2)\Gamma\_{jES}^{(k)}(\mathbf{x}\_2, \mathbf{x}\_1)\right]}{\mathrm{Tr}\left[\Gamma\_{iES}^{(k)}(\mathbf{x}\_1, \mathbf{x}\_1)\right]\mathrm{Tr}\left[\Gamma\_{jES}^{(k)}(\mathbf{x}\_1, \mathbf{x}\_1)\right]}, k = 1, \ldots, 2^{n-1}. \tag{100}$$

Here, {*S* ð Þ*k iES*,*jES*} 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*�<sup>2</sup> *n*(*n*-1) XC/CC-DOI ratios, we could find out a reasonable solution of Eq. (103) by imposing 2*n*�<sup>1</sup> equations to eliminate all the residues {Δ<sup>2</sup> *S*<sup>2</sup> *tot* � �ð Þ*<sup>k</sup>* } from a set of the expected values of { *<sup>S</sup>*<sup>2</sup> *tot* � �ð Þ*<sup>k</sup>* }, which are given by

$$
\left\langle \left\langle \mathbf{S}\_{\text{tot}}^{2} \right\rangle^{(k)} = \left\langle \mathbf{S}\_{\text{tot},x} \right\rangle^{(k)} \left[ \left\langle \mathbf{S}\_{\text{tot},x} \right\rangle^{(k)} + \mathbf{1} \right] + \sum\_{i=1}^{n} \mathbf{S}\_{i} \left( \mathbf{1} - \sigma\_{x,i}^{(k)} \right) + \Delta^{2} \left\langle \mathbf{S}\_{\text{tot}}^{2} \right\rangle^{(k)},\tag{101}
$$

$$\Delta^2 \left< \mathbf{S}\_{\rm tot}^2 \right>^{(k)} = \sum\_{i=1}^n \mathbf{S}\_i^2 \left[ \mathbf{1} - \left( P\_{i \to \mathbf{S}}^{(k)} \right)^2 \right] - \frac{\mathbf{1}}{2} \sum\_{i$$

$$\mathbf{S}\_{i \to \mathbf{S}, \text{iES}}^{(k)} = \frac{4}{n(n-1)} \left\{ \sum\_{i'=1}^{n} \mathbf{S}\_{i'}^{2} \left[ \mathbf{1} - \left( n\_{i \to \mathbf{S}}^{(k)} / 2 \mathbf{S}\_{i'}^{(k)} \right)^{2} \right] \right\} \left( n\_{i \to \mathbf{S}}^{(k)} n\_{j \to \mathbf{S}}^{(k)} \right)^{-1}; k = 1, \ldots, 2^{n-1}, \tag{103}$$

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:

$$\langle H\_{\rm ex} \rangle^{(k)} = -\frac{1}{2} \sum\_{i$$

$$U\_{\rmUHFD-H\_{\rm cr}}^{(k)}[\rho] = T^{(k)}[\rho] + V\_{\rm en}^{(k)}[\rho] + J^{(k)}[\rho] - K\_{\rm XC-Hex}^{(k)}\left[\rho - \sum\_{i=1}^{n} \{\rho\_{i \rm ES} + \rho\_{i \rm Liab}\}\right], \tag{105}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

$$U\_{\rm UHFD}^{(k)}\left[\rho\_a,\rho\_\beta\right] = U\_{\rm UHFD-H\_{ex}}^{(k)}\left[\rho\right] + \left^{(k)}.\tag{106}$$

Notably, one may expect that the *Hex*-less internal energy function defined by Eq. (106) will become almost constant over 2*n*�<sup>1</sup> 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 ESexchange coupling constants, {*Ji,j*; *i*(*j* > *i*) = 1, … ,*n*}, since *n*2*n*�<sup>1</sup> effective spin densities {*niES* (*k*) }, 2*n*�<sup>2</sup> *n*(*n*-1) XC/CC-DOI ratios defined by Eq. (103) and (2*n*�<sup>1</sup> -1) energydifference equations (107) could be quantitatively calculated using UB3LYP/PBS(ε6)/ lacvp\*\* method:

$$
\Delta U\_{\rm UHFD}^{(k)}[\rho\_a, \rho\_\beta] \equiv U\_{\rm UHFD}^{(k)}[\rho\_a, \rho\_\beta] - U\_{\rm UHFD}^{(1)}[\rho\_a, \rho\_\beta] \cong \langle H\_{\rm ex} \rangle^{(k)} - \langle H\_{\rm ex} \rangle^{(1)}; k = 2, \dots, 2n - 1,\tag{107}
$$

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

$$\mathbf{A}^{T}\mathbf{A}\mathbf{X}=\mathbf{A}^{T}\mathbf{B},\tag{108}$$

$$\mathbf{A}=-\frac{1}{2}\begin{bmatrix}\Delta\left[\left(\mathbf{1}-\mathbf{S}\_{\mathrm{iES,ES}}^{(2)}\right)n\_{\mathrm{iES}}^{(2)}n\_{\mathrm{iES}}^{(2)}\right] & \cdots & \Delta\left[\left(\mathbf{1}-\mathbf{S}\_{\mathrm{(n-1)ES,nES}}^{(2)}\right)n\_{\mathrm{(n-1)ES}}^{(2)}n\_{\mathrm{iES}}^{(2)}\right] \\ \vdots & \ddots & \vdots \\ \Delta\left[\left(\mathbf{1}-\mathbf{S}\_{\mathrm{(n-1)ES,nES}}^{\left(\mathbf{2}^{n-1}\right)}\right)n\_{\mathrm{(n-1)ES}}^{\left(\mathbf{2}^{n-1}\right)}\right] & \cdots & \Delta\left[\left(\mathbf{1}-\mathbf{S}\_{\mathrm{(n-1)ES,nES}}^{\left(\mathbf{2}^{n-1}\right)}\right)n\_{\mathrm{(n-1)ES}}^{\left(\mathbf{2}^{n-1}\right)}\right] \end{bmatrix},\tag{109}$$

$$\mathbf{B} = (B\_2, \dots, B\_{2^{n-1}}), B\_k = \Delta U\_{\text{UHFD}}^{(k)} - \Delta U\_{\text{UHFD}-H\_{\text{ox}}}^{(k)}; k = 2, \dots, 2^{n-1}, \tag{110}$$

$$\mathbf{X} = \left(\mathbf{X}\_1, \dots, \mathbf{X}\_{n(n-1)/2}\right) \equiv \left(\mathbf{X}\_{\vec{\eta}}\right); \mathbf{X}\_{\vec{\eta}} = f\_{i,j}, \tag{111}$$

where *A<sup>T</sup>* is the transpose of *A*. Thus, we get a unique solution:

$$J\_{i,j} = \left[ \left( \mathbf{A}^T \mathbf{A} \right)^{-1} \mathbf{A}^T \mathbf{B} \right]\_{ij}. \tag{112}$$

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

$$J\_{1,2} \cong \frac{2\Delta\mathbf{U}\_{\mathrm{UHFD}}^{(2)}}{\left(\mathbf{1} - \mathbf{S}\_{\mathrm{1ES,2ES}}^{(1)}\right)n\_{\mathrm{1ES}}^{(1)}n\_{\mathrm{2ES}}^{(1)} - \left(\mathbf{1} - \mathbf{S}\_{\mathrm{1ES,2ES}}^{(2)}\right)n\_{\mathrm{1ES}}^{(2)}n\_{\mathrm{2ES}}^{(2)}}.\tag{113}$$

In **Table 1**, we show again the results of benchmark-test calculations of the ESexchange coupling constants *J*1,2, designated (*J* ð Þ *m* 1,2*=***<sup>a</sup>**, *J* ð Þ *m* 1,2*=***<sup>c</sup>**, … , *J* ð Þ *m* 1,2*=***<sup>k</sup>**), for 10 biomimetic binuclear Cu, Mn and Fe complexes, named (**a**, **c**, … , **k**), were made using 13 conventional *m*XC/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 *J exp* 1,2*=***<sup>a</sup>**, *J exp* 1,2*=***<sup>c</sup>**, … *:*, *J exp* 1,2*=***k** � �, to show all the excellent quantitative agreements between the theoretical values *J* ð Þ 4 1,2*=***<sup>a</sup>**, *J* ð Þ 4 1,2*=***<sup>c</sup>**, … , *J* ð Þ 4 1,2*=***k** � � 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 geometryoptimization 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 S1 (0) state Mn4CaO5 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 Mn4CaO5 cluster in photosystem II (PSII). The electronabstracting and proton-releasing reactions from the so-called oxygen-evolving complex (OEC) are considered to occur serially via five redox states, called Kok's Si-states (i = 0, 1, … , 4), where S1 is the dark-stable state, and S4 spontaneously decays to the initial S0-state after releasing two protons and evolving dioxygen: the generalized reaction schemes are symbolically given by

$$\mathcal{S}\_1^{(0)} \stackrel{-\epsilon}{\rightarrow} \mathcal{S}\_2^{(+1)} \stackrel{-\epsilon \leftarrow H^+}{\rightarrow} \mathcal{S}\_3^{(+1)} \stackrel{-\epsilon \leftarrow 2H^+}{\rightarrow} \mathcal{S}\_4^{(0)} \stackrel{+2H\_2O \rightarrow O\_2}{\rightarrow} \mathcal{S}\_0^{(0)} \stackrel{-\epsilon \leftarrow H^+}{\rightarrow} \left\{ \mathcal{S}\_1^{(0)} \mathcal{O} \right\},\tag{114}$$

where the figure *k* in the superfix parentheses of S*<sup>i</sup>* (*k*) represents a formal charge of the *i*th OEC, �e above an arrow (!) indicate one electron transfer from OEC to P680(+), an oxidized PSII reaction center intermittingly generated by every �<sup>10</sup> <sup>μ</sup><sup>s</sup> light-pulse, �*H*<sup>+</sup> 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 S1 (0)-state Mn4CaO5 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 S1 (0)-state. Although we can'nt exclude the possibility that the XFEL model [21] may reflect a photo-reduced S0 (�1) state of the S1 (0)-state Mn4CaO5 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 S1 (0)-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 S1A and S1B in [23], which coexist with quasi-degenerate lowest energies in the respective S1 (0)-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

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

model [26]. The highest *Tc* 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. *θ*/kB 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-*Tc* than those of normal metals: (1) higherfrequency-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,

$$\lim\_{\kappa \to 0^{+}} \frac{e^{2}}{r} e^{-\kappa r} = \int V\_{\kappa}(\mathbf{q}) e^{iq \cdot \mathbf{r}} \mathbf{d} \mathbf{q}; V\_{\kappa^{\prime}}(\mathbf{q}) = \lim\_{\kappa \to 0^{+}} \frac{4\pi e^{2}}{q^{2} + \kappa^{2}},\tag{115}$$

where *<sup>κ</sup>* is called "Thomas-Fermi wave number" and the limit of *<sup>κ</sup>* ! <sup>0</sup><sup>+</sup> implies a bare-Coulomb interaction. In order to treat such a many-body effect for Blochelectrons 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

$$
\Omega = -pV = U - \mu\_F \overline{N} - \theta \Sigma. \tag{116}
$$

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* � *μFN* 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-Tc superconductors [27].

In contrast to the conventional idea of repulsive Coulomb force, the principal exchange-correlation energy Hamiltonian *<sup>H</sup>*<sup>b</sup> <sup>0</sup> *UHFD* 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 *<sup>f</sup>*(*ϵ<sup>k</sup>* � *<sup>μ</sup>F*), designated *<sup>σ</sup> <sup>B</sup>***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:

$$\mathbf{^B}\_B \mathbf{\mathbf{\color{red}{U}}} \mathbf{\color{red}{U}} = \mathbf{\_B\mathbf{\color{red}{U}}} \mathbf{\color{red}{U}} \mathbf{\color{red}{U}} \mathbf{\color{red}{F}} + \mathbf{\_B\mathbf{\color{red}{U}}} \mathbf{\color{red}{U}} \mathbf{\color{red}{U}} + \mathbf{\_B\mathbf{\color{red}{U}}} \mathbf{\color{red}{U}} \mathbf{\color{red}{U}} \mathbf{\color{red}{D}} + \cdots,\tag{117}$$

$$\hat{\boldsymbol{a}}\_{B}\hat{\boldsymbol{\Omega}}\_{\text{UHFD}}^{0}\boldsymbol{\cong}\_{\text{B}}\hat{\boldsymbol{H}}\_{\text{UHFD}}^{0} - \mu\_{\text{F}}\sum\_{\mathbf{k}\sigma}^{\text{''M}}\boldsymbol{a}\_{\mathbf{k}\sigma}^{\dagger}\boldsymbol{a}\_{\mathbf{k}\sigma} = \sum\_{\sigma=\mathbf{a}\_{2}\beta}\sum\_{\mathbf{k}\sigma}^{\text{''M}}[\boldsymbol{c}\_{\mathbf{k}\sigma} - \mu\_{\text{F}}]\boldsymbol{a}\_{\mathbf{k}\sigma}^{\dagger}\boldsymbol{a}\_{\mathbf{k}\sigma},\tag{118}$$

$$\hat{\mathbf{a}}\_{B}^{\text{T}}\hat{\mathbf{a}}\_{\text{UHFD}}^{\text{T}} = \sum\_{\sigma=\alpha,\beta} \sum\_{\substack{\mathbf{k}\sigma,\mathbf{k}\sigma}}^{\mathbf{\tilde{c}}\mathbf{M}} V\_{\text{UHFD}}^{\mathbf{k}\sigma,\mathbf{k}\overline{\sigma}} V\_{\text{k}\sigma,\mathbf{k}\overline{\sigma}} a\_{\mathbf{k}\sigma}^{\dagger} a\_{\mathbf{k}\overline{\sigma}}^{\dagger} a\_{\mathbf{k}'\overline{\sigma}} a\_{\mathbf{k}'\sigma},\tag{119}$$

$$\begin{split} V\_{\text{UHFD}}^{\mathsf{h}\sigma,\mathsf{k}\overline{\sigma}} \equiv & -e^{2} \iint \frac{d\boldsymbol{r}\_{1} d\boldsymbol{r}\_{2}}{2r\_{12}} \Phi\_{\mathsf{k}\cdot\sigma}^{\*}(\boldsymbol{r}\_{2}) \Phi\_{\mathsf{k}\overline{\sigma}}^{\*}(\boldsymbol{r}\_{1}) \Phi\_{\mathsf{k}\overline{\sigma}}(\boldsymbol{r}\_{2}) \Phi\_{\mathsf{k}\cdot\sigma}(\boldsymbol{r}\_{1}) \\ & \times (\boldsymbol{\sigma}\_{+,\mathsf{k}\ \sigma,\mathsf{}} \boldsymbol{\sigma}\_{-,\mathsf{k}\prime\overline{\sigma}} + \boldsymbol{\sigma}\_{-,\mathsf{k}\prime\sigma} \boldsymbol{\sigma}\_{+,\mathsf{k}\prime\overline{\sigma}}) \end{split} \tag{120}$$

$$\approx -\lim\_{\kappa \to 0^{+}} \frac{4\pi e^{2}V}{\left|\mathbf{k} - \mathbf{k}'\right|^{2} + \kappa^{2}} (\sigma\_{+,\mathbf{k}'\sigma,\mathbf{}}\sigma\_{-,\mathbf{k}'\overline{\sigma}} + \sigma\_{-,\mathbf{k}'\sigma}\sigma\_{+,\mathbf{k}'\overline{\sigma}}),\tag{121}$$

where we put *α* ¼ *β* and *β* ¼ *α*, and used the plane-wave approximation for the conductive Bloch-wave functions: <sup>Φ</sup>*k<sup>σ</sup>*ð Þ *<sup>r</sup>*<sup>2</sup> <sup>≈</sup> exp <sup>½</sup>*ik*•*r*2], <sup>Φ</sup><sup>∗</sup> *<sup>k</sup> <sup>σ</sup>*ð Þ *r***<sup>2</sup>** ≈exp ½�*ik*•*r*2] et c. to make the double integrations on the coordinate vectors *r*<sup>1</sup> and *r*<sup>2</sup> after transformed into *R*<sup>12</sup> = (*r*<sup>1</sup> + *r*2)/2 and *r*<sup>12</sup> = *r*2-*r*1, together with introducing the screened-Coulomb damping factor exp.(�*κr*12). Note that the volume of the system *V* appears from the integral on the center-of-mass coordinate *R*12.

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

$${}\_{B}\Omega\_{UHFD}^{2} = \sum\_{\sigma=\alpha,\beta} \sum\_{k\_{\rm l\sigma}}^{\tilde{\bf s}\mathbf{M}} \sum\_{\mathbf{k}\_{\rm l\sigma}}^{\tilde{\bf s}\mathbf{M}} \sum\_{k\_{\rm l\sigma}}^{\tilde{\bf s}\mathbf{M}} \sum\_{k\_{\rm l\sigma}}^{\tilde{\bf s}\mathbf{M}} \frac{\mathbf{V}^{\mathbf{k}\_{1}\mathbf{k}\_{\rm l\sigma}} \left[ V^{\mathbf{k}\_{1}\mathbf{k}\_{\rm l\sigma}}\_{UHFD} \right]^{\*}}{\mathbf{c}\_{\mathbf{k}\_{1s}} + \mathbf{c}\_{\mathbf{k}\_{\rm l\sigma}} - \mathbf{c}\_{\mathbf{k}\_{2\sigma}} - \mathbf{c}\_{\mathbf{k}\_{2\sigma}}} \tag{122}$$

$${}\_{B}\times f(\boldsymbol{c}\_{\mathbf{k}\_{\rm l\sigma}} - \boldsymbol{\mu}) f\left(\boldsymbol{c}\_{\mathbf{k}'\mathbf{j}\_{\rm l\sigma}} - \boldsymbol{\mu}\right) [1 - f(\boldsymbol{c}\_{\mathbf{k}\_{2\sigma}} - \boldsymbol{\mu})] \left[1 - f\left(\boldsymbol{c}\_{\mathbf{k}'\mathbf{j}\_{\rm l\sigma}} - \boldsymbol{\mu}\right)\right],$$

$$V^{\mathbf{k},\mathbf{k}\nu}\_{UHFD} = -\lim\_{\kappa \to 0^{+}} \frac{4\pi e^{2}V}{\left|\mathbf{k} - \mathbf{k}'\right|^{2} + \kappa^{2}}.\tag{123}$$

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

$$\begin{split} \hat{\mathbf{a}}\_{B}^{2} \hat{\mathbf{a}}\_{\text{UHFD}}^{2} &= \sum\_{\sigma=\alpha,\beta} \sum\_{k\_{1\sigma}}^{\hat{\mathbf{a}}} \sum\_{k\_{2\sigma}}^{\hat{\mathbf{a}}} \sum\_{k\_{3\sigma}}^{\hat{\mathbf{a}}} \sum\_{k\_{4\sigma}}^{\hat{\mathbf{a}}} \sum\_{k\_{5\sigma}}^{\hat{\mathbf{a}}} V\_{\text{UHFD}}^{k\_{1}k\_{1};k\_{2},k\_{2}} \\ &\times \hat{a}\_{k'\varpi} \hat{a}\_{k'\varpi}^{\dagger} \hat{a}\_{k\nu} \hat{a}\_{k\omega}^{\dagger} \hat{a}\_{k'\varpi}^{\dagger} \hat{a}\_{k'\varpi} \hat{a}\_{k\nu}^{\dagger} \hat{a}\_{k\omega}^{\dagger} \hat{a}\_{k\nu} \end{split} \tag{124}$$
 
$$V\_{UHFD}^{k\_{1}k\_{2};k\_{2},k\_{2}} = \lim\_{\kappa \to 0^{+}} \frac{\left(4\pi\kappa^{2}V\right)^{2}}{\left(|k\_{1} - k'\_{1}|^{2} + \kappa^{2}\right)\left(|k\_{2} - k'\_{2}|^{2} + \kappa^{2}\right)(c\_{k\_{1}} + c\_{k\_{1}} - c\_{k\_{2}} - c\_{k\_{2}})} < 0. \tag{125}$$

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*0 *<sup>i</sup><sup>β</sup>* = �*kiβ*; *i* = 1, 2 and multiplying twice the state number in each spherical-shell volume, 4π*kF* 3 *ωD/μF*, as given by *N*(*kF*) ≈ 4π*kF* 3 *ωD*/*μ<sup>F</sup>* (2π) 3 *V* = *kF* 3 *ωD*/2π<sup>2</sup> *VμF*):

$$\left\|\widehat{\boldsymbol{a}}\_{B}\right\|\_{\text{UHFD}}^{2} \approx \sum\_{\sigma=\mathbf{a}\_{\tau}\boldsymbol{\beta}} \sum\_{\mathbf{k}\_{1\sigma}}^{\hat{\mathbf{a}}\mathbf{M}} \sum\_{\mathbf{k}\_{2\sigma}}^{\hat{\mathbf{a}}\mathbf{M}} V\_{\mathrm{UHFD}}^{\mathbf{k}\_{1}-\mathbf{k}\_{1};\mathbf{k}\_{2\tau}-\mathbf{k}\_{2}} \hat{\boldsymbol{a}}\_{-\mathbf{k}\_{2\sigma}} \hat{\boldsymbol{a}}\_{-\mathbf{k}\_{2\sigma}}^{\dagger} \hat{\boldsymbol{a}}\_{\mathbf{k}\_{2\sigma}} \hat{\boldsymbol{a}}\_{\mathbf{k}\_{2\sigma}}^{\dagger} \hat{\boldsymbol{a}}\_{-\mathbf{k}\_{1\sigma}}^{\dagger} \hat{\boldsymbol{a}}\_{-\mathbf{k}\_{2\sigma}} \hat{\boldsymbol{a}}\_{\mathbf{k}\_{1\sigma}}^{\dagger} \hat{\boldsymbol{a}}\_{\mathbf{k}\_{1\sigma}} \right. \tag{126}$$

$$V\_{UHFD}^{k\_1, -k\_1; k\_2, -k\_2} \approx \lim\_{\kappa \to 0^+} \frac{(4\pi\epsilon^2)^2 \left(k\_F^3 \alpha\_D / 2\pi^2 \mu\_F\right)^2}{2\left(4k\_1^2 + \kappa^2\right)\left(4k\_2^2 + \kappa^2\right)\left(\epsilon\_{k\_1} - \epsilon\_{k\_2}\right)} < 0;\tag{127}$$

which is an attractive potential under the BCS restrictions:

$$-\alpha\_{\rm D} < \varepsilon\_{\rm k\_1} - \mu\_F < 0 < \varepsilon\_{\rm k\_2} - \mu\_F < \alpha\_{\rm D},\tag{128}$$

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

where ω<sup>D</sup> is the Debye frequency, and note that the matrix element Eq. (126) does not contain *V* and it is proportional to *ω<sup>D</sup>* 2 /(*ϵ<sup>k</sup>*<sup>2</sup> � *ϵ<sup>k</sup>*<sup>1</sup> ), which diverges as *ϵ<sup>k</sup>*<sup>2</sup> � *ϵ<sup>k</sup>*<sup>1</sup> ð Þ! 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 Mn4CaO5-cluster in photosystem II. This one-body quasiparticle world picture has been completed by our two revolutionary discoveries of the *principal exchange-correlation energy functional,* that is, a non-local exchangecorrelation 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 Mn4CaOx 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.

#### **Acknowledgements**

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.

### **Author details**

Masami Kusunoki Meiji University, Kawasaki, Japan

\*Address all correspondence to: kusunoki@meiji.ac.jp

© 2023 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

*Distinct Roles of the Principal Exchange-Correlation Energy and the Secondary Correlation… DOI: http://dx.doi.org/10.5772/intechopen.111746*

#### **References**

[1] Hohenberg P, Kohn W. Inhomogeneous electron gas. Physical Review. 1964;**136**(B):864-871. DOI: 10.1103/PhysRev.136.B864

[2] Kohn W, Sham LJ. Self-consistent equations including exchange and correlation effects. Physical Review. 1965;**140**(A):1133-1137. DOI: 10.1103/ PhysRev.140.A1133

[3] Sham LJ. Exchange and correlation in density-functional theory. Physical Review B. 1985;**32**:3876-3882. DOI: 10.1103/PhysRevB.32.3876

[4] Parr RG, Yang W. Density-functional theory of atoms and molecules. In: Breslow R, Goodenough JB, Halpern J, Rowlinson JS, editors. International Series of Monographs on Chemistry. New York: Oxford Science; 1989. vol. 16, pp. 1-333. ISBN: 0-19-509276-7

[5] Dreizler RM, Gross EKU. Density Functional Theory: An Approach to the Quantum Many-Body Problem. Berlin: Springer; 1990. pp. 1-302. DOI: 10.1007/ 978-3-642-86105-5

[6] Kock W, Holthausen MC. A Chemist's Guide to Density Functional Theory. Weinheim: WILEY-VCH; 2000. pp. 1-294. ISBN: 3-527-29918-1

[7] Becke AD. Perspective: Fifty years of density-functional theory in chemical physics. Journal of Chemical Physics. 2014;**140**(18):18A301. DOI: 10.1063/ 1.4869598

[8] Sundararaman R, Goddard WA, Alias TA. Grand canonical electronic density functional theory: Algorithms and applications to electrochemistry. Journal of Chemical Physics. 2017; **146**:114104-114114. DOI: 10.1063/ 1.4978411

[9] Evangelista FA. Perspective: Multireference coupled cluster theories of dynamical electron correlation. Journal of Chemical Physics. 2018;**149**: 030901-030914. DOI: 10.1063/ 1.5039496

[10] Mardirossian N, Head-Gordon M. Thirty years of density functional theory in computational chemistry: An overview and extensive assessment of 200 density functionals. Molecular Physics. 2017;**115**:2315-2372. DOI: 10.1080/00268976.2017.1333644

[11] Kitakawa CK, Maruyama T, Oonari J, Mitsuta Y, Kawakami T, Okumura M, et al. Linear response functions of densities and spin densities for systematic Modeling of the QM/MM approach for mono- and poly-nuclear transition metal systems. Molecules. 2019;**24**:821-849. DOI: 10.3390/molecules24040821

[12] Zhou C, Hermes MR, Wu D, Bao JJ, Pandharkar R, King DS, et al. Electronic structure of strongly correlated systems: Recent developments in multiconfiguration pair density functional theory and multiconfiguration nonclassical-energy functional theory. Chemical Science. 2022;**13**:685-7706. DOI: 10.1039/ d2sc01022d

[13] Kusunoki M. Heisenberg spin Hamiltonian derived from a multiple grand canonical spin density functional theory with a principal nonlocal exchangecorrelation energy functional. Journal of Physical Society of Japan. 2022;**91**:014702- (1-28). DOI: 10.7566/JPSJ.91.014702

[14] Mulliken RS. Electronic population analysis on LCAO–MO molecular wave functions. I. Journal of Chemical Physics. 1955;**23**:1833-1840. DOI: 10.1063/ 1.1740588

[15] Lee C, Yang W, Parr RG. Development of the colic-Salvetti correlation-energy formula into a functional of the electron density. Physical Review B. 1988;**98**:785-789. DOI: 10.1103/ PhysRevB.37.785

[16] Becke AD. A new mixing of Hartree-Fock and local density-functional theories. Journal of Chemical Physics. 1993;**98**: 1372-1378. DOI: 10.1063/1.464304

[17] Najafpour MM, Renger G, Holynska M, Moghaddam AN, Aro EM, Carpentier R, et al. Manganese compounds as water-oxidizing catalysts: From the natural water-oxidizing complex to Nanosized manganese oxide structures. Chemicl Reviews. 2016;**116**:2886-2936. DOI: 10.1021/acs.chemrev.5b00340

[18] Dexheimer SL, Klein MP. Detection of a paramagnetic intermediate in the S1 state of the photosynthetic oxygenevolving complex. Journal of American Chemical Society. 1992;**114**:2821-2826. DOI: 10.1021/ja00034a010

[19] Yamauchi T, Mino H, Matsukawa T, Kawamori A, Ono T. Parallel polarization electron paramagnetic resonance studies of the S1-state manganese cluster in the photosynthetic oxygen-evolving system. Biochemistry. 1997;**36**:7520-7526. DOI: 10.1021/bi962791g

[20] Campbell KA, Gregor W, Pham DP, Peloquin JM, Debus RJ, Britt RD. The 23 and 17 kD extrinsic proteins of photosystem II modulate the magnetic properties of the S1-state manganese cluster. Biochemistry. 1998;**37**: 5039-5045. DOI: 10.1021/bi9800552

[21] Suga M, Akita F, Hirata K, Ueno G, Murakami H, Nakajima Y, et al. Native structure of photosystem II at 1.95A° resolution viewed by femtosecond X-ray pulses. Nature. 2015;**517**:99-103. DOI: 10.1038/nature13991

[22] Kusunoki M. S1-state Mn4Ca complex of photosystem II exists in equilibrium between the two most-stable isomeric substates: XRD and EXAFS evidences. Journal of Photochemistry and Photobiology B. 2011;**104**:100-1010. DOI: 10.1016/j.jphotobiol.2011.03.002

[23] Bardeen J, Cooper LN, Schrieffer JR. Theory of superconductivity. Physical Review. 1957;**108**:1175-1204. DOI: 10.1103/PhysRev.108.1175

[24] Schrieffer JR. Theory of Superconductivity. New York: WA Benjamin, Inc; 1964. pp. 1-282 ISBN: 9780429975332

[25] Fröhlich H. Theory of the superconducting state. I. the ground zero of temperature. Physical Review. 1950; **79**:845-846. DOI: 10.1103/PhysRev. 79.845

[26] Bogoliubov NN. A new method in the theory of superconductivity. Soviet Physics JETP. 1958;**34**:41-46

[27] Drozdov AP, Eremets MI, Troyan IA, Ksenofontov V, Shylin SI. Conventional superconductivity at 203 kelvin at high pressures in the sulfur hydride system. Nature. 2015;**525**:73. DOI: 10.1038/nature14964
