Computational Physics Course 2020 summer

Brief introduction of lecturers
Schedule
In the first half of the course, I will introduce the density functional theory (DFT) and thesoftware package “Quantum Espresso”. The contents are as follows
1. Density functional theory (DFT)
2. Relation with many-body theory
3. Quantum espresso
In the second half Prof. Arita will show you
4. Introduction to strong correlation
5. Dynamical mean field theory (DMFT)
6. DFT+DMFT

1

General introduction
I will introduce the density functional theory (DFT) and its extensions, which have beendeveloped as a tool to simulate the electronic structure of materials and to predict propertiesof materials. The importance of DFT is in its practicality; the theory provides an approximate
but practical way to understand the electronic structure, which has had, and will continue tohave, impact in the research of materials. (In my definition, materials mean those matterstargeted by condensed matter physics, biophysics, or chemical physics.) In this context, I
consider it important for beginners to practice a simulation, so that students are requested topractice a few simulations.
Why is practicality so important?
To understand the electronic structures, one can use the many-body wave function, or
alternatively the Green’s function. The many-body wave function contains full information
but has impractically large degrees of freedom. This may be understood by counting thenumber of ways to assign ten electrons to 𝑀 orbitals, which amounts to 𝑀10 ! The number
can be reduced by considering symmetry but is still too large. Of course, all the degrees of
freedom are not essentially important, and it is possible to considerably reduce the number
with information theory for example: Research in that direction is on-going. The need forpractical theory has led to the invention (development) of the density functional theory. This
point was emphasized by W. Kohn in his Nobel prize lecture:https://www.nobelprize.org/uploads/2018/06/kohn-lecture.pdf
The invent of DFT turned out to be revolutionary in that the electronic structure
calculation has made an easy access of nonexperts, including experimentalists as well as
theoreticians working on abstract models. From a practical point of view, experimentalists areoften required to do a DFT simulation to support their experiment by reviews of journals andsimilar things happen for theoreticians. Experiencing a simulation at the start of his/herresearch life should be advantageous.
Simulation in this course
I will spend a few hours to show students how to install a simulation software and how to
use it. I will focus on the program package called Quantum Expresso (QE). This is currently
the most popular one in the world and there are many how-to’s available via the internet.
Students of this course are requested to use their own computers and will be guided to do a
few simulations.

2Introduction to Density functional theory (DFT)
I will follow the article written by J. Toulouse
http://www.lct.jussieu.fr/pagesperso/toulouse/enseignement/introduction_dft.pdf
which was prepared for the summer school in France.
Some recommended articles
The most popular introduction to DFT was written by Parr and Yang.
R. G. Parr and W. Yang, Density-Functional Theory of Atoms and Molecules (Oxford
University Press, New York, 1989).
which covers the development until mid-1980. Practical aspect of DFT simulation is summarizedinW. Koch and M. C. Holthausen, A Chemist’s Guide to Density Functional Theory (WileyVCH, New York, 2001),
which covers the technology on the exchange-correlation that had been developed until 2000
although further developments were made thereafter. Most recent textbook can be found, forexample, in
T. Helgaker, P. Jørgensen and J. Olsen, Density–Functional Theory: A Convex Treatment(WileyBlackwell, 2016)

The starting Hamiltonian
Here we on the Born-Oppenheimer approximation and non-relativistic approximation for
𝑁-electron systems. Then, the target of study is the electronic Hamiltonian, consisting of the
kinetic energy of electrons, electron-electron interaction, and the Coulomb field from thenuclear charge
̂el = 𝑇̂el + 𝑊
̂ee + 𝑉̂ne ,
𝐻

(1.4)

for the details see Eqs. (1.1)-(1.4) of the article by J. Toulouse (hereafter, I will call it just“article”). Instead of the wave function |𝛹⟩, the density 𝑛(𝒓) can be used to characterize theground state: This is what DFT is telling about.
What does it mean?
The density is defined by𝑛(𝒓) = 𝑁 ∭|𝛹(𝒙, 𝒙2 , ⋯ , 𝒙𝑁 )|2 𝑑𝜎𝑑𝒙2 ⋯ 𝑑𝒙𝑁 ,(1.6)

where 𝒙 is the space-spin coordinate and satisfies ∫ 𝑛(𝒓)𝑑𝒓 = 𝑁. This three-dimensional
quantity is the “descriptor” of the 𝑁 -electron system although |𝛹⟩ is, of course, thedescriptor as well. This surprising statement was shown by Hohenberg and Kohn about fifty3years ago. (This is surprising when considering large difference in the degrees of freedom
between the density and the wave function.)
Hohenberg-Kohn theorem
This theorem relates or maps the electron density 𝑛(𝒓) to the external potential 𝑣(𝒓). In
considering the mapping, it is trivial that mapping from 𝑣 to 𝑛 exits simply because, given
the external potential 𝑣, one can solve, in principle, the Schrödinger equation and the density
𝑛 is obtained from Eq. (1.6). What is nontrivial is the reverse mapping, i.e., mapping from 𝑛
to 𝑣 , which can be derived as follows. (Before going in detail on the Hohenberg-Kohn
theorem, I note that this may be counterintuitive in that the precise measurement of theelectron density yields uniquely the position of nuclei.)
The proof can be given by different ways; (1) two-step proof by contradiction and (2)proof using the Levy’s constraint search. Here I will show the first one only.
Two-step proof by contradiction (a modified version)
We will derive a contradiction for the assumption that there exist two local potentials(differing by more than an additive constant) which have the same ground-state density.
From the two local potentials, 𝑣1 (𝒓) and 𝑣2 (𝒓) which differ more than an additivê1 = 𝑇̂ + 𝑊
̂𝑒𝑒 + V
̂1 and 𝐻
̂2 = 𝑇̂ + 𝑊
̂𝑒𝑒 + V
̂2 .
constant, one can construct the Hamiltonians, 𝐻
Let us assume that they share the same ground state |𝛹⟩. Then, one gets(𝑉̂1 − 𝑉̂2 )|𝛹⟩ = (𝐸1 − 𝐸2 )|𝛹⟩,(1.9)

namely𝑁∑[𝑣1 (𝒓𝑖 ) − 𝑣2 (𝒓𝑖 )]𝛹(𝒙1 , 𝒙2 , ⋯ , 𝒙𝑁 ) = (𝐸1 − 𝐸2 )𝛹(𝒙1 , 𝒙2 , ⋯ , 𝒙𝑁 ).

(1.10)

𝑖=1

holds globally in the phase space, leading to 𝑣1 (𝒓) − 𝑣2 (𝒓) = const when the potentials
reasonably well behave almost everywhere so that 𝛹(𝒙1 , 𝒙2 , ⋯ , 𝒙𝑁 ) ≠ 0 except for somepoints. This contradicts with the initial hypothesis.
̂1 and 𝐻
̂2,
We will then assume that different states |𝛹1 ⟩ and |𝛹2 ⟩ corresponding to 𝐻
respectively, have the same density 𝑛(𝒓). These states are different and therefore |𝛹2 ⟩̂1,cannot be a ground-state wave function of 𝐻
̂1|𝛹1 ⟩ < ⟨𝛹2 |𝐻
̂1|𝛹2⟩ = ⟨𝛹2 |𝐻
̂2 + 𝑉̂1 − 𝑉̂2 |𝛹2⟩ = 𝐸2 + ∫ [𝑣1 (𝒓) − 𝑣2 (𝒓)]𝑛(𝒓)𝑑𝒓 (1.11)
𝐸1 = ⟨𝛹1|𝐻
holds. In parallel, we can show the following:
𝐸2 < 𝐸1 + ∫ [𝑣2 (𝒓) − 𝑣1 (𝒓)]𝑛(𝒓)𝑑𝒓(1.12)

Then, we arrive at an inconsistent inequality 𝐸1 + 𝐸2 < 𝐸1 + 𝐸2 .
Therefore, we are lead to the conclusion that 𝑣1 (𝒓) and 𝑣2 (𝒓) are the same except for4the additive constant and the mapping from 𝑛 to 𝑣 holds.
From the discussion made in the proof it is also clear that the ground-state energy 𝐸 canbe obtained as the minimum of the expectation value of the Hamiltonian, or there exist a totalenergy functional 𝐸[𝑛] which is given aŝee + ∫ 𝑣[𝑛](𝒓)𝑛(𝒓)𝑑𝒓|𝛹[𝑛]⟩}.
𝐸[𝑛] = min{⟨𝛹[𝑛]|𝑇̂ + 𝑊(1.15)

We can also define a functional called Hohenberg-Kohn functional aŝee|𝛹[𝑛]⟩}.
𝐹[𝑛] = min{⟨𝛹[𝑛]|𝑇̂ + 𝑊(1.14)

𝑛𝑛

Once the functional has been obtained, the ground state density and the ground state energy
can be obtained without solving the Schrödinger equation. The advantage is very large.

𝒗-representability
It is reminded that 𝛹[𝑛] is not uniquely determined when the ground-state is degenerate.
Nevertheless, the functionals 𝐹 and 𝐸 are uniquely determined. These values can be
determined by finding a minimum once the functionals 𝐹[𝑛] and 𝐸[𝑛] are established. This,
however, does not mean that the statement is of practical importance, especially because thesearch needs to be restricted to those represented from the ground-state wave function via
Eq. (1.6), namely those that satisfy the 𝑣-representable condition. All three-dimensional
functions are not necessarily 𝑣-representable!
Levy-Lieb constrained minimization
The Hohenberg-Kohn functional 𝐹[𝑛] can be made more practical when we introduce analternative definition (Levy and Lieb)̂ee|𝛹⟩
𝐹[𝑛] = min⟨𝛹|𝑇̂ + 𝑊(1.17)

𝛹→𝑛in terms of the constrained minimization of Ψ. This should provide the same functional asthat provided by Eq. (1.14) but does not require the 𝑣-representability, so that the minimum
can be found by simply searching all functions that satisfy ∫ 𝑛(𝒓)𝑑𝒓 = 𝑁: namely assuming
𝑁 -representability. Of course, the obtained density does not necessarily satisfy the 𝑣 representability, which needs to be checked afterwards. When Eq. (1.17) is introduced, it isstraightforward to show that the ground-state energy can be obtained by minimizing 𝑛(𝒓):̂ee + 𝑉̂en|𝛹⟩ = min min⟨𝛹|𝑇̂ + 𝑊
̂ee + 𝑉̂en|𝛹⟩𝐸0 = min⟨𝛹|𝑇̂ + 𝑊
𝛹

𝑛

𝛹→𝑛

̂ee|𝛹⟩ + ∫ 𝑣(𝒓)𝑛(𝒓)𝑑𝒓}= min {min⟨𝛹|𝑇̂ + 𝑊𝑛

𝛹→𝑛= min{𝐹[𝑛] + ∫ 𝑣(𝒓)𝑛(𝒓)𝑑𝒓}𝑛Additional comment of Hohenberg-Kohn theorem
5(1.18)

It should be emphasized that, so far, more than fifty years were spent in trying to find areasonable functional form appropriate for materials; the accuracy is still not enough for many
materials. In the early stage, people searched a form for the kinetic energy functional 𝑇[𝑛],but the efforts have not been so successful. Instead, method of Kohn-Sham was introduced aswill be explained below. People have also searched a form for the electron-electron interaction,𝑉ee [𝑛]; equivalently, a form for the exchange-correlation functional that is define as 𝑉ee [𝑛] =𝑉H [𝑛] + 𝑉xc [𝑛], where the first term is the Hartree-energy defined as
1 𝑛(𝒓)𝑛(𝒓′ )𝑉H [𝑛] = ∫𝑑𝒓𝑑𝒓′,|𝒓 − 𝒓′ |2
and is sometimes called as the classical Coulomb.
Exercise
To understand the difficulty of constructing the kinetic energy functional, let us test the one
constructed using the Thomas Fermi (TF) model. Note that TF model is appropriate forsystems where the electron density has a small spatial dependence, i.e., for a nearly uniform
system. According this model, the kinetic functional is given as𝑇=

35/3 𝜋 4/3∫ 𝑛(𝒓)5/3 𝑑𝒓
10

for a three-dimensional system and as𝑇=

𝜋2 𝐿∫ 𝑑𝑥 𝑛(𝑥)3𝑑𝑥
6 0

for a one-dimensional system. From this, one can naively define the kinetic energy density
functional as𝑡[𝑛(𝑥)] =𝜋2𝑛(𝑥)3 .
6

This is a local density approximation to the kinetic energy.
The integration can be done analytically for a 1D potential well of infinite wall height,
where the wave function has a form𝑛𝜋𝑥
𝜑𝑛 (𝑥) = √2/𝐿 sin (),
𝐿and the kinetic energy is given as
𝑁

𝐿

𝑛

0

2 1 𝑛𝜋 2𝑛𝜋𝑥𝜋21
𝑇 = ∑ ∫ ( ) sin2 () 𝑑𝑥 = 2 𝑁(𝑁 + 1) (𝑁 + ).
𝐿 2 𝐿
𝐿
6𝐿
2
When using the Thomas Fermi model, the kinetic energy is approximated as

6

𝑇

TF

𝐿

𝑁

0

𝑛

3

𝜋2
2
𝑛𝜋𝑥
𝜋2
9
3
=
∫ { ∑ sin2 ()} 𝑑𝑥 = 2 𝑁 (𝑁 2 + 𝑁 + ).
6
𝐿
𝐿
6𝐿
8
8

Therefore, the ratio 𝑇 TF /𝑇 behaves as a function of the density 𝑁/𝐿 likewhich is reasonably well for higher density, or large value for 𝑁/𝐿, but deteriorates withdecreasing density. The reason can be found by comparing the TF kinetic density and theexact one.
0.025
0.020
0.015
0.010
0.005
10

20

30

40

50

The values for 𝑁 = 10, 𝑁/𝐿 = 6 are comparable in the region between the potential walls but
are quite different near the end points. This reflects the general difficulty of constructing the
kinetic energy density functional in the region where the density rather abruptly changes.
Note that the potential changes very abruptly near the nuclear position as well and thedensity also changes in the core region. This is the reason why Thomas Fermi model fails; themodel cannot reproduce the shell structure of an atom!
There have been much efforts to develop the kinetic energy functional as can be foundfrom the textbook written by Parr and Yang. None of the effects have been very successful as
far as I know. The reason may be due to the complexity of the functional. There is a trial toconstruct the functional 𝑇[𝑛] by machine learning (Brockherde et atl. “Bypassing the KohnSham equations with machine learning”, Nat. Comm. 8, 872 (2017)).
Kohn-Sham method (modern version)
To overcome the problem of constructing the kinetic energy functional, Kohn and Sham
proposed a method to replace the kinetic energy functional with that of a non-interacting
7system that shares the same density distribution 𝑛(𝒓). More precisely, one assumes that there
always exists a non-interacting system under the influence of an effective potential 𝑣eff (𝒓)whose ground-state electron density is the same as that of the system of interest. Since thewave function of the non-interacting system can be relatively easily computed numerically,
the kinetic energy functional can also be easily constructed, which will be denoted as 𝑇𝑠 [𝑛].
The proposition of Kohn and Sham is to use 𝑇𝑠 [𝑛] instead of 𝑇[𝑛] and approximate 𝐹[𝑛] as
𝐹[𝑛] = 𝑇𝑠 [𝑛] + 𝑉H [𝑛] + 𝑉xc [𝑛],where
𝑇𝑠 [𝑛] = min⟨Φ|𝑇̂|Φ⟩ = ⟨Φ[𝑛]|𝑇̂|Φ[𝑛]⟩.

(1.21)

Φ→𝑛

The last equation reflects the fact that the non-interacting system can be represented by the
Slater determinant Φ. (Here we use the Hartree-Fock state as the reference system but it canbe anything in principle. Importantly, the reference system should be such that the groundstate wave function can be prepared.) (You may wonder if the kinetic energy can beapproximated by that of Hartree-Fock; this is not the case always especially for stronglyinteracting systems.)
By tracing Eq. (1.18) reversely, we find that the total-energy can be obtained by searchingfor the density representable by the Slater determinant that minimizes the HK functional as
𝐸0 = min{𝐹[𝑛] + ∫ 𝑣(𝒓)𝑛(𝒓)𝑑𝒓} ⟹ min {min⟨Φ|𝑇̂|Φ⟩ + 𝐸Hxc [𝑛] + ∫ 𝑣(𝒓)𝑛(𝒓)𝑑𝒓}
𝑛

𝑛

Φ→𝑛

= min {min⟨Φ|𝑇̂ + 𝑉̂ne |Φ⟩ + 𝐸Hxc [𝑛Φ ]}𝑛

Φ→𝑛= min{⟨Φ|𝑇̂ + 𝑉̂ne|Φ⟩ + 𝐸Hxc [𝑛Φ ]}.
Φ

(1.22)

By minimizing the functional with respect to the Slater determinant Φ, the density can beobtained as 𝑛Φ . Note that the Hartree-exchange-correlation functional can be decomposed
into
𝐸Hxc [𝑛] = 𝐸H [𝑛] + 𝐸x [𝑛] + 𝐸c [𝑛]wherêee|Φ[𝑛]⟩ − 𝐸H [𝑛],𝐸x [𝑛] ≡ ⟨Φ[𝑛]|𝑊(1.26)

̂ee|Ψ[𝑛]⟩ − ⟨Φ[𝑛]|𝑇̂ + 𝑊
̂ee|Φ[𝑛]⟩.
𝐸c [𝑛] ≡ ⟨Ψ[𝑛]|𝑇̂ + 𝑊(1.27)

andFrom these definitions, we find that the difference in the true kinetic energy functional andthat of the non-interacting system, 𝑇[𝑛] − 𝑇𝑠 [𝑛], is absorbed by the correlation functional andthus the correlation functional is comprised of the kinetic energy contribution and potential
energy contribution.
𝐸c [𝑛] = 𝑇c [𝑛] + 𝑈c [𝑛],where8(1.27)

𝑇c [𝑛] ≡ ⟨Ψ[𝑛]|𝑇̂|Ψ[𝑛]⟩ − ⟨Φ[𝑛]|𝑇̂|Φ[𝑛]⟩̂ee|Ψ[𝑛]⟩ − ⟨Φ[𝑛]|𝑊̂ee|Φ[𝑛]⟩.
𝑈c [𝑛] ≡ ⟨Ψ[𝑛]|𝑊Note that this is a rather old definition of the exchange-correlation functional but recently
this “dirty” definition, dirty in that the kinetic energy and the potential energy having different
scaling properties are mixed, is avoided by introducing the adiabatic connection formula, as
explained below.
Kohn-Sham equation
It is well described in many textbooks on DFT how to determine the Slater determinant.
It is determined by solving a Schrödinger-like equation called Kohn-Sham equation. The
equation can be derived by rewriting the minimization with respect to Φ by the minimization
with respect to the single-electron orbitals 𝜓𝑖 that constitute the Slater determinant. Within
the scheme of the unrestricted closed-shell description of the Slater determinant, 𝜓𝑖 (𝒙) =𝜙𝑖 (𝒓)𝜒𝜎𝑖 (𝜎), the functional to be minimized is𝑁

𝐸[{𝜙𝑖 }] = ∑ ∫ 𝜙𝑖∗ (𝒓) (−𝑖

𝛁2+ 𝑣ne (𝒓)) 𝜙𝑖 (𝒓)𝑑𝒓 + 𝐸Hxc [𝑛]2(1.28)

which is very similar to the Hartree-Fock theory. The variation is usually done by keeping theorthonormal condition and the relation𝑁𝑛(𝒓) = ∑|𝜙𝑖 (𝒓)|2 .

(1.29)𝑖

When the variation is done using the method of Lagrange multiplier,𝑁ℒ[{𝜙𝑖 }] = 𝐸[{𝜙𝑖 }] − ∑ 𝜀𝑖 (∫ 𝜙𝑖∗ (𝒓)𝜙𝑖 (𝒓)𝑑𝒓 − 1) ,(1.30)𝑖as𝛿ℒ[{𝜙𝑖 }]= 0,𝛿𝜙𝑖∗ (𝒓)(1.31)

we can obtain the differential equation1
𝛿𝐸Hxc [𝑛](− ∇2 + 𝑣ne (𝒓)) 𝜙𝑖 (𝒓) += 𝜀𝑖 𝜙𝑖 (𝒓).
2𝛿𝜙𝑖∗ (𝒓)(1.32)The left-hand side of the equation can be rewritten as1
𝛿𝐸Hxc [𝑛] 𝛿𝑛(𝒓′ ) ′1 2(− ∇2 + 𝑣ne (𝒓)) 𝜙𝑖 (𝒓) + ∫∗ (𝒓) 𝑑𝒓 = (− ∇ + 𝑣ne (𝒓) + 𝑣Hxc (𝒓)) 𝜙𝑖 (𝒓).
′
)
2
𝛿𝑛(𝒓
𝛿𝜙𝑖
2
Therefore, the differential equation becomes1(− ∇2 + 𝑣ne (𝒓) + 𝑣Hxc (𝒓)) 𝜙𝑖 (𝒓) = 𝜀𝑖 𝜙𝑖 (𝒓).
2This is called as the Kohn-Sham equation.
9(1.35)

Self-consistent field
Note that the effective potential 𝑣eff (𝒓) ≡ 𝑣ne (𝒓) + 𝑣Hxc (𝒓) is a functional of the density𝑛(𝒓), so that the Kohn-Sham equation needs to be solved by keeping the self-consistency
between the density used to construct the effective potential and that resulting from (1.29).
This makes the calculation much more difficult than the one-electron Schrödinger equation.
At this point, it is very instructive to compare the Kohn-Sham equation and the Hartree-Fockequation: The latter is1(− ∇2 + 𝑣ne (𝒓) + 𝑣H (𝒓) + 𝑣x (𝒓)) 𝜙𝑖 (𝒓) = 𝜀𝑖 𝜙𝑖 (𝒓),
2
where
𝑁

𝑣x (𝒓)𝜙𝑖 (𝒓) = − ∑ ∫𝑗∑𝑗 𝜙𝑗∗ (𝒓′ )𝜙𝑗 (𝒓)𝜙𝑗∗ (𝒓′ )𝜙𝑖 (𝒓′ ) ′
(𝒓)
𝑑𝒓
𝜙
=
−
𝜙𝑖 (𝒓′ )𝑑𝒓′∫𝑗|𝒓 − 𝒓′ ||𝒓 − 𝒓′ |≡ ∫ 𝑣𝑥 (𝒓, 𝒓′ )𝜙𝑖 (𝒓′ )𝑑𝒓′
is non-local. The non-local nature of the Hartree-Fock has a significant meaning in that
contributions of the 𝑖-th orbital to the Hartree potential 𝑣H (𝐫) and the exchange potential𝑣x (𝒓, 𝒓′ ) are the same except for the different sign. That is, the contributions of the 𝑖-th
orbital are canceled out when determining the 𝑖-th orbital. This property of the effective
potential is called self-interaction free.
Local density approximation for the exchange-correlation
The exchange-correlation functional, given by Eqs. (1.26-1.27), are generally non-local;
non-locality of the exchange was shown explicitly in the above and that of the correlation ismore complicated as shown below. Non-locality makes the self-consistent calculation very
complicated. In the early stage of DFT, therefore, people searched for a local form of theexchange-correlation.
The simplest form for the exchange can be found in the literature as the Thomas Fermi
exchange or Dirac exchange as
1

4
3 3 3
𝐸x [𝑛(𝑟)] = − ( ) ∫ 𝑛(𝒓)3 𝑑𝒓.
(1.53)
4 𝜋
Later, elaborate calculation was done for homogeneous electron gas, where the electron

density is constant, using various values for the electron density 𝑛 and the calculated energy
density was related to 𝑛 to find a better local form for the exchange-correlation. The resulting
functional was called “local density approximation (LDA)”. LDA was found to reasonably well10

reproduce the energy of many weakly interacting electrons and was used to predict stability,
phonon, dynamics of materials.
How to deal with spin
You may wonder how to deal with spin. Or, you may alternatively want to use anunrestricted Hartree-Fock form for the non-interacting system (see discussion around Eq.
(1.28)). The corresponding Hohenberg Kohn functional can be defined aŝee|Ψ⟩.
𝐹[𝑛↑ , 𝑛↓ ] = min ⟨Ψ|𝑇̂ + 𝑊
Ψ→𝑛↑ ,𝑛↓(1.57)The search is done for normalized antisymmetric functions Ψ that yields𝑛↑ (𝒓) = 𝑁 ∭ |Ψ(𝒓 ↑, 𝒙2 , ⋯ , 𝒙𝑁 )𝑑𝒙2 ⋯ 𝑑𝒙𝑁(1.58𝑎)𝑛↓ (𝒓) = 𝑁 ∭ |Ψ(𝒓 ↓, 𝒙2 , ⋯ , 𝒙𝑁 )𝑑𝒙2 ⋯ 𝑑𝒙𝑁 .

(1.58𝑏)and

A KS scheme can be derived by decomposing 𝐹[𝑛↑ , 𝑛↓ ] into the kinetic, Hartree, andexchange-correlation parts, where the kinetic part is defined as
𝑇𝑠 [𝑛↑ , 𝑛↓ ] = min ⟨Φ|𝑇̂|Φ⟩.
Φ→𝑛↑ ,𝑛↓(1.60)The Slater determinant is composed from the one-electron orbital of the form
𝜓𝑖 (𝒙) = 𝜙𝑖𝜎𝑖 (𝒓)𝜒𝜎𝑖 (𝜎)and then the KS equation can be derived separately for each spin state1(− ∇2 + 𝑣ne (𝒓) + 𝑣H (𝒓) + 𝑣xc,↑ (𝒓)) 𝜙𝑖↑ (𝒓) = 𝜀𝑖↑ (𝒓)𝜙𝑖↑ (𝒓)2(1.61𝑎)1(− ∇2 + 𝑣ne (𝒓) + 𝑣H (𝒓) + 𝑣xc,↓ (𝒓)) 𝜙𝑖↓ (𝒓) = 𝜀𝑖↓ (𝒓)𝜙𝑖↓ (𝒓)2(1.61𝑏)The exchange-correlation potentials are given in terms of the derivative with respect to theelectron densities as shown in Eqs. (1.62)-(1.63) of the article. It was shown by Oliver and
Perdew that the spin-dependent exchange functional 𝐸x [𝑛↑ , 𝑛↓ ] can be expressed rigorouslyas1
𝐸x [𝑛↑ , 𝑛↓ ] = (𝐸x [2𝑛↑ ] + 𝐸x [2𝑛↓ ]),2(1.64)

because, in the non-relativistic case,𝐸x [𝑛↑ , 𝑛↓ ] = 𝐸x [𝑛↑ , 0] + 𝐸x [0, 𝑛↓ ]holds and, therefore, when using the spin-unpolarized case, 𝑛↑ = 𝑛↓ = 𝑛/2, we get

11(S. 8)

1
𝐸x [𝑛, 0] = 𝐸x [0, 𝑛] = 𝐸x [2𝑛].
2(S. 9)

using the spin-scaling relation which states that the up-spin density and down-spin density
are uncoupled to form the exchange functional.
This formalism for the spin density functional theory suggests a possibility to constructextended formalisms as well. Indeed, by describing the HK functional in terms of the density
and the current density, the current density functional theory was developed and was applied
to electrons under excited states. One can also describe the HK functional using the electricfield and magnetic flux to simulate the electrons under inhomogeneous electromagnetic field.
In those cases, it is not obvious if one can prepare appropriate exchange-correlation functional.
If one wants to follow the procedure adopted in LDA, one needs to do an accurate simulation
under the inhomogeneous conditions and establish the functional relation. This is not always
possible especially when one needs to construct a non-local functional.
Density matrix formalism
As a step to extend LDA, it is convenient to adopt the density matrix formalism or theGreen’s function formalism. Let us, for the moment, focus on the former formalism. In thatcase, we need to investigate up to the second order density matrix because we are studying
electrons interacting through two-body interaction, that is the Coulombic interaction.
The two-body density matrix is defined as𝑛2 (𝒓1 , 𝒓𝟐 ) = 𝑁(𝑁 − 1) ∭|Ψ[𝑛](𝒙1 , ⋯ , 𝒙𝑁 )|2 𝑑𝜎1 𝑑𝜎2 𝑑𝒙3 ⋯ 𝑑𝒙𝑁 ,(2.1)

although factor 1/2 is multiplied in many literatures. This density matrix corresponds to theprobability of finding two electrons simultaneously at 𝒓1 and 𝒓2 . This density matrix isparticularly important because the electron-electron interaction functional can be written
using 𝑛2 as
1
𝑛 (𝒓 , 𝒓 )̂ee|Ψ[𝑛]⟩ = ∬ 2 1 2 𝑑𝒓1 𝑑𝒓2 .
⟨Ψ[𝑛]|𝑊|𝒓1 − 𝒓2 |
2

(2.2)

The two-body density matrix is conventionally decomposed into the product of the density
and the rest as
𝑛2 (𝒓1 , 𝒓2 ) = 𝑛(𝒓1 )𝑛(𝒓2 ) + 𝑛2,xc (𝒓1, 𝒓2 ).

(2.3)

The second term in the right-hand side is called the exchange-correlation pair density, which
can be further decomposed as𝑛2,xc (𝒓1 , 𝒓2 ) = 𝑛(𝒓1 )𝑛xc (𝒓1, 𝒓2 ),(2.4)

where 𝑛xc is called the exchange-correlation hole. We call this as a hole because of theexistence of the following sum rule12∫ 𝑛xc (𝒓1 , 𝒓2 )𝑑𝒓2 = −1.

(2.6)

That is, the value is negative when integrated over the whole space. The electron-electroninteraction energy can be rewritten as1𝑛(𝒓1 )𝑛(𝒓2 )1𝑛(𝒓1 )𝑛xc (𝒓1 , 𝒓2 )
̂ee|Ψ[𝑛]⟩ = ∬𝑑𝒓1 𝑑𝒓2 + ∬𝑑𝒓1 𝑑𝒓2 .
⟨Ψ[𝑛]|𝑊|𝒓1 − 𝒓2 ||𝒓1 − 𝒓2 |22This indicates that the (potential energy contribution) of exchange-correlation energy is givenby the interaction with its exchange-correlation hole.
Exercise
Let us solve a Hooke’s atom problem and compute the exchange-correlation hole. The system
consists of two electrons confined in a Harmonic potential and the Hamiltonian is written as1111̂ = − ∇12 − ∇22 + 𝑘(𝑟12 + 𝑟22 ) +𝐻|𝒓1 − 𝒓2 |222
which is a model He atom. Using the center of mass coordinate 𝑹 = (𝒓1 + 𝒓2 )/2 and therelative coordinate 𝒖 = 𝒓2 − 𝒓1 , the Hamiltonian becomes111
̂ = − ∇2𝑹 − ∇2𝒖 + 𝑘𝑅 2 + 𝑘𝑢2 + ,𝐻
4
4𝑢so that the wave function can be decomposed into a product Ψ(𝐫1 , 𝒓2 ) = 𝜒(𝑹)Φ(𝒖) and the
Schrödinger equation becomes1(− ∇2𝑹 + 𝑘𝑅 2 ) 𝜒(𝑹) = 𝐸𝑹 𝜒(𝑹)41
1(−∇2𝒖 + 𝑘𝑢2 + ) Φ(𝒖) = 𝐸𝒖 Φ(𝒖).
4𝑢
The first differential equation provides the well-known solution32√𝑘 4 −√𝑘𝑅2
𝜒𝒏 (𝑹) =() 𝑒𝐻𝑛𝑥 (21/2 𝑘1/4 𝑅𝑥 )𝐻𝑛𝑦 (21/2 𝑘1/4 𝑅𝑦 )𝐻𝑛𝑧 (21/2 𝑘1/4 𝑅𝑧 )√2𝑛𝑥 𝑛𝑦𝑛𝑧 𝑛𝑥 ! 𝑛𝑦 ! 𝑛𝑧 ! 𝜋
1

3
𝐸𝑹,𝒏 = (𝑛𝑥 + 𝑛𝑦 + 𝑛𝑧 + ) √𝑘
2
The second one is centrosymmetric and the solution should have a form𝑅𝑙 (𝑢)𝑌𝑙𝑚 (𝑢̅)and the equation for the radial part is
(−

1 𝑑
𝜕
𝑙(𝑙 + 1) 1 2 1(𝑢2 ) ++ 𝑘𝑢 + ) 𝑅𝑙 (𝑢) = 𝐸𝒖,𝑙 𝑅𝑙 (𝑢)
2
𝑢 𝑑𝑢
𝑑𝑢
𝑢2
4
𝑢

By rewriting 𝑅𝑙 (𝑢)

13

√𝑘 2

𝑒 − 4 𝑢 𝑇𝑙 (𝑢)𝑅𝑙 (𝑢) =,𝑢
we obtain(−𝑑2
𝑑 𝑙(𝑙 + 1) √𝑘 1
+ √𝑘𝑢
+
+
+ ) 𝑇𝑙 (𝑢) = 𝐸𝒖,𝑙 𝑇𝑙 (𝑢).
2
𝑑𝑢
𝑑𝑢
𝑢2
2
𝑢

Following the Frobenius method, let us express the regular solution for 𝑇𝑙 (𝑢) as∞𝑇𝑙 (𝑢) = 𝑢

𝑙+1∑ 𝑎𝑖 𝑢𝑖
𝑖=0

and require that the recursion of the coefficient 𝑎𝑖−(𝑖 + 𝑙 + 1)(𝑖 + 𝑙)𝑎𝑖 𝑢𝑖+𝑙−1 + √𝑘(𝑖 + 𝑙 + 1)𝑎𝑖 𝑢𝑖+𝑙+1 + 𝑙(𝑙 + 1)𝑎𝑖 𝑢𝑖+𝑙−1 + (√𝑘− 𝐸𝒖,𝑙 ) 𝑎𝑖 𝑢𝑖+𝑙+12+ 𝑎𝑖 𝑢𝑖+𝑙 = 0
to terminate at some point of 𝑖. The recursion is rewritten as[−(𝑖 + 𝑙 + 2)(𝑖 + 𝑙 + 1)𝑎𝑖+1 + √𝑘(𝑖 + 𝑙)𝑎𝑖−1 + 𝑙(𝑙 + 1)𝑎𝑖+1 + (√𝑘− 𝐸𝒖,𝑙 ) 𝑎𝑖−1 + 𝑎𝑖 ] 𝑢𝑖−1 = 0
2

so that𝑎𝑖+1 =𝑎1 =1
𝑎𝑖 + (√𝑘 (𝑖 + 𝑙 + 2) − 𝐸𝒖,𝑙 ) 𝑎𝑖−1
(𝑖 + 1)(𝑖 + 2𝑙 + 2)

𝑎0
𝑎0
1
3
, 𝑎2 =(+ √𝑘 (𝑙 + ) − 𝐸𝒖,𝑙 ).
2(𝑙 + 1)2(2𝑙 + 3) 2(𝑙 + 1)
2
𝑎3 =

5
𝑎2 + (√𝑘 (𝑙 + 2) − 𝐸𝒖,𝑙 ) 𝑎1

3(2𝑙 + 4)
The termination condition depends on the value of √𝑘, and the simplest solution can befound for those parameters which vanish 𝑎2 and 𝑎3 : The result is√𝑘 =12(𝑙 + 1)𝐸𝒖,𝑙 =

2𝑙 + 5.
4(𝑙 + 1)andNamely, for the harmonic potential of strength 𝑘 = 4(𝑙 + 1)2, analytical solution can be found.
(Other solutions can be found by vanishing 𝑎 with higher indices.)
The radial wave function is2

𝑅𝑙 (𝑢) = 𝑢𝑙 (1 +𝑢1−𝑢) 𝑒 8(𝑙+1)2(𝑙 + 1)and thus, the wave function for the relative coordinate is

14

2

𝑢
1
−
̅ ).
𝑢) 𝑒 8(𝑙+1) 𝑌𝑙𝑚 (𝒖2(𝑙 + 1)

Φ(𝒖) = 𝑢𝑙 (1 +When multiplied by the state with the lowest energy for the center of mass 𝜒𝟎 (𝑹) is
3

2

2

𝑢
4 − 𝑅
1
1
−
̅)
𝜒𝟎 (𝑹)Φ(𝒖) = () 𝑒 2(𝑙+1) 𝑢𝑙 (1 +𝑢) 𝑒 8(𝑙+1) 𝑌𝑙𝑚 (𝒖𝜋(𝑙 + 1)2(𝑙 + 1)32

2

4 −(𝑟1 +𝑟2 )1
1̅̅̅̅̅̅̅̅̅̅).
=() 𝑒 4(𝑙+1) |𝒓2 − 𝒓1 |𝑙 (1 +|𝒓 − 𝒓1 |) 𝑌𝑙𝑚 (𝒓2 − 𝒓1𝜋(𝑙 + 1)2(𝑙 + 1) 2

After normalization,
Ψ(𝒓1 , 𝒓2 ) = 𝑁𝑙 𝑒(𝑟 2 +𝑟 2 )− 1 2
4(𝑙+1) |𝒓 − 𝒓 |𝑙 (1 +211̅̅̅̅̅̅̅̅̅̅).
|𝒓 − 𝒓1 |) 𝑌𝑙𝑚 (𝒓2 − 𝒓12(𝑙 + 1) 2

Pair density is then obtained as𝑛2 (𝒓1 , 𝒓2 ) = 𝑁𝑙2 𝑒−(𝑟12 +𝑟22 )2(𝑙+1) |𝒓 − 𝒓 |2𝑙 (1 +
2
1

𝑙

2
1
2
|𝒓2 − 𝒓1 |) ∑ |𝑌𝑙𝑚 (𝒓̅̅̅̅̅̅̅̅̅̅)|2 − 𝒓12(𝑙 + 1)𝑚=−𝑙2
2
1 +𝑟2 )
2𝑙 + 1 −(𝑟2(𝑙+1)2= 𝑁𝑙𝑒|𝒓2 − 𝒓1 |2𝑙 (1 +4𝜋21|𝒓 − 𝒓1 |)2(𝑙 + 1) 2

after taking average over the degenerate states. The electron density is obtained as2

2𝑁𝑙2 ∭

2

2
1 +𝑟 )
2𝑙 + 1 −(𝑟2(𝑙+1)1|𝒓 − 𝒓1 |) 𝑑𝒓1𝑒|𝒓 − 𝒓1 |2𝑙 (1 +4𝜋
2(𝑙 + 1)2
−2𝑟𝑢 cos 𝜃
2𝑙 + 1 −2𝑟 +𝑢2(𝑙+1)1𝑒𝑢2𝑙 (1 +𝑢) 2𝜋𝑢2 sin 𝜃 𝑑𝑢𝑑𝜃.
4𝜋2(𝑙 + 1)2

= 2𝑁𝑙2 ∭2

Therefore,𝑟2 𝜋 7𝑟2
1
𝑟2𝑛(𝒓) = 2𝑁02 (𝑒 −𝑟 + 𝑒 − 2 √ ( + + (𝑟 + ) erf ( )))
2 4 4
𝑟
√2

for 𝑙 = 0 and2𝑟𝑟3√𝜋𝑒 − 4 (156𝑟 + 36𝑟 3 + 𝑟 5 + 8(12 + 12𝑟 2 + 𝑟 4 ) erf (2)−22𝑛(𝒓) = 2𝑁1 (6𝑒 2 (10 + 𝑟 ) +)8𝑟
𝑟2

for 𝑙 = 1. The exchange-correlation hole is then obtained from𝑛xc (𝒓2 ; 𝒓1 ) =𝑛2 (𝒓1 , 𝒓2 ) − 𝑛(𝑟1 )𝑛(𝑟2 ).
𝑛(𝑟1 )For 𝑙 = 0, it looks like

15

Here the first electron is located at 𝒓1 = (1,0,0) and the coordinate of the second electron istaken as 𝒓2 = (𝑥, 𝑦, 0). The aim of KS-DFT is to establish the relationship between theelectron density (right; unnormalized) and the (exchange-)correlation hole (left). For 𝑙 = 1,
the exchange-correlation hole looks like𝑙 = 1 is an excited state, so that this is outside the scope of DFT strictly.
Returning to the main steam, let us continue the discussion on the exchange-correlation hole.
We can separate the exchange-correlation functional into exchange and correlation. This canbe done by using the pair density of the Kohn-Sham single-determinant Φ[𝑛] as
𝑛2,KS (𝒓1 , 𝒓2 ) = 𝑛(𝒓1 )𝑛(𝒓2 ) + 𝑛(𝒓1 )𝑛x (𝒓1, 𝒓2 )The exchange hole can be written using the Kohn-Sham orbitals as

16

𝑁

2

𝑛x (𝑟1 , 𝑟2 )𝑛(𝑟1 ) = − ∑ |∑ 𝜙𝑗∗ (𝑟1 𝜎)𝜙𝑗 (𝑟2 𝜎)| ,𝜎𝑗so that it is negative everywhere.
In weakly correlation cases, the exchange hole is much larger in magnitude than thecorrelation hole in general. Therefore, the major task of DFT in those systems is to find an
appropriate functional form for the exchange plus relatively small contribution from thecorrelation. In practice, finding the form for the correlation is much more difficult. Note that,when using the above form for the exchange, which is implicitly the functional of the electron
density via the Kohn-Sham orbitals, one needs to try to find the correlation functional only;
when, on the other hand, when using a local (and approximate) value, one usually tries to findthe form for the exchange-correlation without separating it into pieces.
Adiabatic connection
Here I introduce the important technique called adiabatic connection, with which toredefine the exchange-correlation. In my view, this is extremely important in connecting DFT
to formally rigorous many-body theories.
In this scheme, we introduce a path connecting between the non-interacting KS systemand the target system of interest while keeping the density constant. The density is kept equal
to the exact one. The Hamiltonian changes via the path aŝ 𝜆 = 𝑇̂ + 𝜆𝑊
̂ee + V̂𝜆,𝐻

(2.18)
𝜆̂
where V is the external local potential operator playing a role to keep the density constant.
Note that the Hamiltonian is the Kohn-Sham one when 𝜆 is zero and is the physical one when
it is 1. Along the line we have discussed, we can introduce the HK functional for each 𝜆.
̂ee|Ψ⟩ = ⟨Ψ λ [𝑛]|𝑇̂ + 𝜆𝑊
̂ee|Ψ λ [𝑛]⟩𝐹 𝜆 [𝑛] = min⟨Ψ|𝑇̂ + 𝜆𝑊
(2.19)
Ψ→𝑛This functional can be decomposed into𝜆𝐹 𝜆 [𝑛] = 𝑇𝑠 [𝑛] + 𝐸H𝜆 [𝑛] + 𝐸xc[𝑛](2.20)

1𝑛(𝒓1 )𝑛(𝒓2 )𝐸H𝜆 [𝑛] = ∬𝜆𝑑𝒓1 𝑑𝒓2 = 𝜆𝐸H [𝑛],2|𝒓1 − 𝒓2 |(2.21)

̂ee|Φ[𝑛]⟩ − 𝐸H𝜆 [𝑛] = 𝜆𝐸x [𝑛]𝐸x𝜆 [𝑛] = ⟨Φ[𝑛]|𝜆𝑊(2.22)

̂ee|Ψ λ [𝑛]⟩ − ⟨Φ[𝑛]|𝑇̂ + 𝜆𝑊
̂ee|Φ[𝑛]⟩.
𝐸c𝜆 [𝑛] = ⟨Ψ λ [𝑛]|𝑇̂ + 𝜆𝑊

(2.23)

withThe last equation can be used to derive𝜕𝐸c𝜆 [𝑛]̂ee|Ψ λ [𝑛]⟩ − ⟨Φ[𝑛]|𝑊̂ee|Φ[𝑛]⟩,= ⟨Ψ λ [𝑛]|𝑊𝜕𝜆
where the Hellmann Feynman theorem has been used. Then, we can derive17(2.24)

1

̂ee|Ψ λ [𝑛]⟩𝑑𝜆 − ⟨Φ[𝑛]|𝑊̂ee|Φ[𝑛]⟩.
𝐸𝑐 [𝑛] = ∫⟨Ψ λ [𝑛]|𝑊(2.25)

0

Therefore, only the interaction energy (not the kinetic energy) contribute to the formulation;
we have thus successfully avoided the treatment of the kinetic energy part. We can further
rewrite the correlation energy as
1

1𝑛(𝒓1 )𝑛𝑐𝜆 (𝒓1 , 𝒓2 )𝐸𝑐 [𝑛] = ∫ 𝑑𝜆 ∬𝑑𝒓1 𝑑𝒓2|𝒓𝟏 − 𝒓2 |2

(2.26)

1𝑛(𝒓1 )𝑛̅𝑐 (𝒓1 , 𝒓2 )𝐸𝑐 [𝑛] = ∬𝑑𝒓1 𝑑𝒓2|𝒓𝟏 − 𝒓2 |2

(2.27)

0

or

where
1

𝑛̅𝑐 (𝒓1 , 𝒓2 ) = ∫ 𝑑𝜆𝑛𝑐𝜆 (𝒓1 , 𝒓2 ).
0This indicates that it is the 𝜆-averaged density that needs to be handled. Later, I will showthat one can handle it using the fluctuation dissipation theorem.
Fractional number of electrons
With DFT, it is plausible if one can consider the system under finite temperatures, where
the system contains thermally excited state and thus the theory is no more a ground-state one.
It is also plausible if one can study a subsystem interacting with the host material, and in
addition, exchanging electrons with the host. Extension of DFT in those directions werediscussed in the past. I will briefly explain the density functional theory for fractional numberof electrons as a step towards those problems.
Let us denote the number as 𝒩 which is also written as 𝒩 = 𝑁 − 1 + 𝑓 with 𝑁 ∈ 𝕀 and
0 ≤ 𝑓 ≤ 1. The total-energy is described aŝee + 𝑉̂ne )],𝐸0𝒩 = min Tr[Γ̂(𝑇̂ + 𝑊̂Γ

(2.30)

using the ensemble density matrix defined asΓ̂ = (1 − 𝑓)|Ψ𝑁−1 ⟩⟨Ψ 𝑁−1 | + 𝑓|Ψ 𝑁 ⟩⟨Ψ 𝑁 |.

(2.31)(One may use many states with different number of electrons as well.) The variation is donefor Ψ 𝑁−1 and Ψ 𝑁 , by which the ensemble density matrix is obtained asΓ̂0 = (1 − 𝑓)|Ψ0𝑁−1 ⟩⟨Ψ0𝑁−1 | + 𝑓|Ψ0𝑁 ⟩⟨Ψ0𝑁 |.
Then the total-energy can be described as18(2.32)

𝐸0𝒩 = (1 − 𝑓)𝐸0𝑁−1 + 𝑓𝐸0𝑁 .
From this definition, it is clear that the total-energy has a derivative discontinuity as
𝑑𝐸0𝒩
𝐸0𝑁 − 𝐸0𝑁−1 ≡ −𝐼𝑁 when 𝑁 − 1 < 𝒩 < 𝑁= { 𝑁+1𝑑𝒩
𝐸0 − 𝐸0𝑁 ≡ −𝐴𝑁 when 𝑁 < 𝒩 < 𝑁 + 1
Usually, the ionization 𝐼𝑁 and the affinity 𝐴𝑁 are different, so that the electronic chemical
potential 𝑑𝐸0𝒩 ⁄𝑑𝒩 ≡ 𝜇 is discontinuous at the integer number: This is called integer
discontinuity. Note that the fundamental gaps, or the HOMO-LUMO gap, is defined as
𝑁𝐼𝑁 − 𝐴𝑁 ≡ 𝐸gap.

DFT for fractional number of electrons assumes existence of universal HK functional̂ee )].
𝐹[𝑛] = min Tr[Γ̂(𝑇̂ + 𝑊
̂→𝑛Γ(2.38)

This can be decomposed, in the KS theory, into𝐹[𝑛] = 𝑇𝑠 [𝑛] + 𝐸Hxc [𝑛](2.39)

where the Kohn-Sham non-interacting kinetic energy is given as
𝑇𝑠 [𝑛] = min Tr[Γ̂𝑠 𝑇̂],̂s →𝑛
Γand the ensemble non-interacting density matrix isΓ̂𝑠 = (1 − 𝑓) |Φ𝑁−1,𝑓 ⟩⟨Φ𝑁−1,𝑓 | + 𝑓|Φ𝑁,𝑓 ⟩⟨Φ𝑁,𝑓 |.

(2.41)The total-energy can also be described using the Kohn-Sham orbitals in a way very similar tothe one shown above, but here I show the resulting equations. By introducing partial
occupation number of each KS orbitals 𝑛𝑖 , the total-energy can be described as𝑁

1
𝐸 = ∑ 𝑛𝑖 ∫ 𝜙𝑖∗ (𝒓) (− ∇2 + 𝑣ne (𝒓)) 𝜙𝑖 (𝒓)𝑑𝒓 + 𝐸Hxc [𝑛],2(2.42)𝑖with the density𝑁𝑛(𝒓) = ∑ 𝑛𝑖 |𝜙𝑖 (𝒓)|2 .

(2.43)

𝑖

Note that the occupation number satisfies the condition that 𝑛𝑖 = 1 for 𝑖 ≤ 𝑁 − 1 and 𝑛𝑁 =𝑓. The KS equation is thereby unchanged1(− ∇2 + 𝑣s (𝒓)) 𝜙𝑖 (𝒓) = 𝜀𝑖 𝜙𝑖 (𝒓)2(2.44)

𝛿𝐸Hxc [𝑛].
𝛿𝑛(𝒓)(2.45)𝑣s (𝒓) = 𝑣ne (𝒓) +It should be reminded, however, the derivative with respect to the electron density becomesambiguous when change in the number of electrons is allowed in the formalism. Therefore
𝛿𝐸Hxc [𝑛] should be understood in terms of19𝛿𝐸Hxc [𝑛] = ∫ (𝛿𝐸Hxc [𝑛]+ const) 𝛿𝑛(𝒓)𝑑𝒓𝛿𝑛(𝒓)(2.46)and the change in the number of electrons should be carefully taken into account to determine
the constant, as will be shown in the next subsection. From the above equations, we can derive𝜕𝐸= 𝜀𝑖 .
∂𝑛𝑖(2.47)

This is called Janak theorem. Recall the Koopmans theorem of the Hartree-Fock, which issimilar to Eq. (2.47) but is different in that the differentiation is replaced by finite differencein the Koopmans.
LUMO energy and the derivative discontinuity
From the Janak theorem, the derivative at a number slightly larger than an integer is,
(

𝜕𝐸0𝒩
𝑁+𝛿
)
= 𝜀HOMO,𝜕𝒩 𝑁+𝛿

(2.54)

where the right-hand side is equal to the HOMO energy at 𝑁 + 1-th orbital. Since Eq. (2.54)
is identical to the definition of the affinity, we have
𝑁+𝛿
𝜀HOMO
= −𝐴𝑁 .

(2.55)

This means that the affinity can be calculated using the HOMO energy. But, actually, this isunrelated to the LUMO energy. To see it, let us write the HOMO orbital in terms of the KS
orbitals as
1
𝑁+𝛿
𝑁+𝛿 (𝒓)∗𝑁+𝛿 (𝒓)𝑑𝒓.
𝜀HOMO= ∫ 𝜙HOMO(− ∇2 + 𝑣𝑠𝑁+𝛿 (𝒓)) 𝜙HOMO
2(2.56)The naive counterpart will be the LUMO energy of a system with a number slightly smaller
than 𝑁. The corresponding LUMO energy is
1𝑁−𝛿𝑁−𝛿 (𝒓)∗𝑁−𝛿 (𝒓)𝑑𝒓.
𝜀LUMO= ∫ 𝜙LUMO(− ∇2 + 𝑣𝑠𝑁−𝛿 (𝒓)) 𝜙LUMO
2(2.57)

With decreasing 𝛿 → 0+ , 𝑣𝑠𝑁+𝛿 (𝒓) and 𝑣𝑠𝑁−𝛿 (𝒓) approach the potentials that are equal toeach other except for the additive constant, 𝑣𝑠𝑁+𝛿 (𝒓) − 𝑣𝑠𝑁−𝛿 (𝒓) ≡ Δ𝑁
xc . Therefore,
𝑁+𝛿
𝑁−𝛿
𝜀HOMO
= 𝜀LUMO+ Δ𝑁
xc(2.58)

holds. Namely, after taking the limit (𝛿 → 0),𝑁
𝜀LUMO= −𝐴𝑁 − Δ𝑁
xc .

(2.60)

This allows to rewrite the fundamental gap in terms of the derivative discontinuity as
𝑁
𝑁
𝑁
𝐸gap
= 𝜀LUMO− 𝜀HOMO+ Δ𝑁
xc .

(2.63)

The constant appeared in the derivation of HK theorem thus turns out to be the derivative
discontinuity.

20

Existing approximate functionals
There are many different functionals of different level of approximations. Usually, there isa trend that the more complex the functional is the more accurate the result is, but this does
not mean superiority in all respects. A functional can be more accurate in one respect but isless accurate in the other respects. Let us briefly survey popular functionals.
Local density approximation (LDA): The exchange-correlation of LDA can be written as
uniform gas

LDA [𝑛]
𝐸xc= ∫ 𝑛(𝑟)𝜀xc(𝑛(𝒓))𝑑𝒓.

𝜀 is obtained by doing a Monte Carlo simulation of uniform electron gas of a fixed electron
density. The result of the exchange-correlation energy is fitted to a form1

3 3 3 1
− ( ) 𝑛3 (𝒓) + 𝜀𝑐 (𝑟𝑠 (𝑛(𝒓)))4 𝜋
𝐴 ln 𝑟𝑠 + 𝐵 + 𝐶 𝑟𝑠 ln 𝑟𝑠𝑏𝜀𝑐 (𝑟𝑠 ) = { 𝑎
+ 3/2
𝑟𝑠 𝑟

when 𝑟𝑠 ≤ 𝑟𝑠,0
when 𝑟𝑠 > 𝑟𝑠,0

𝑠

3

1
3

where 𝑟𝑠 = (4𝜋𝑛) .
Generalized gradient approximation (GGA): LDA is based on the uniform electron gas andtherefore is more accurate for a system having more uniform electron density. This means
LDA favors too much metallic systems over insulating systems and favors also condensed
phases over isolated phases like an atom. This is a serious drawback when comparing phasesof different electronic structure. It is also known that the adsorption energy on the surface oratomization energy is too large and is outside of the required chemical accuracy. As a method
to overcome it, gradient of the electron density is incorporated into the functional. A naturalform for the gradient corrected exchange-correlation functional is obtained by expanding by
the dimensionless parameter ∇𝑛(𝒓)⁄𝑛(𝒓)4/3 as2
4 ∇𝑛(𝒓)
GEA [𝑛]LDA [𝑛]𝐸xc= 𝐸xc+ ∫ 𝐶xc (𝑛(𝒓))𝑛(𝒓)3 (4 ) 𝑑𝒓,𝑛(𝒓)3(3.10)The expansion parameter, however, can be large in some part of the target system making the
expansion breaks down. In this context, people have tried to construct a functional so that
GGA [𝑛]𝐸xc= ∫ 𝑓(𝑛(𝒓), ∇𝑛(𝒓))𝑑𝒓.

(3.11)

Note that LDA uses a local approximation to the exchange-correlation, in this sense GGA issaid to be a semi-local approximation scheme for the exchange-correlation. Some works nicely
at crystalline phases and others at surfaces, and so on. For details, please read the textbooks
21

shown at the beginning.
GGA[#1] B88 exchange functional: Becke 88 was constructed to reproduce the nonlocal
exchange functional (Eq. (1.26)), such as the asymptotic behavior of the exchange energy
per particle and the value of the nonlocal exchange energy of rare-gas atoms.
GGA[#2] LYP correlation functional: Lee-Yang-Parr (LYP) correlation functional wasconstructed based on an approximation of the Hartree-Fock pair density.
GGA[#3] PW91 exchange-correlation functional: Perdew-Wang 91 exchangecorrelation functional is based on a model of the exchange hole (Eq. (2.14)) and thecoupling-constant-averaged correlation hole (Eq. (2.29)).
GGA[#4] PBE exchange-correlation functional: Perdew-Burke-Ernzerhof exchangecorrelation functional is a simplified PW91.
Meta GGA: Meta GGA is based on the description of the functional using the second
derivative of the density and the newly introduced quantity 𝜏(𝒓) defined by Eq. (3.14).

mGGA
𝐸xc
= ∫ 𝑓(𝑛(𝒓), ∇𝑛(𝒓), ∇2 𝑛(𝒓), 𝜏(𝒓))𝑑𝒓(3.13)

𝑁1
𝜏(𝒓) = ∑|∇𝜙𝑖 (𝒓)|22(3.14)𝑖𝜏(𝒓) is constructed by the square of the KS gradient. This quantity is distinct from the local
functional of electron density and its derivatives. The ground-state energy can then be given
formally as
𝐸0 = min⟨Φ|𝑇̂ + 𝑉̂ne |Φ⟩ + 𝐸H [𝑛Φ ] + 𝐸xc [𝑛Φ , 𝜏Φ ].
Φ

Therefore, it is a slight extension of the original KS-DFT formalism. Meta GGA is generally
superior to GGA but requires dedicated basis functions or finer real-space mesh.
Hybrid DFT: Becke proposed to mix the nonlocal exchange energy, or the Hartree-Fockexchange energy,
𝑁𝜎

𝐸xHF = −∗ (𝒓 )𝜙 (𝒓 )𝜙 ∗ (𝒓 )𝜙 (𝒓 )
𝜙𝑖𝜎
1
1
𝑗𝜎 𝟏
𝑖𝜎 2
𝑗𝜎 2
∑∑∬
𝑑𝒓1 𝑑𝒓2 ,2|𝒓1 − 𝒓2 |𝜎𝑖𝑗with the existing forms for the exchange and correlation, such as22(3.18)

3H
𝐸xc= 𝑎𝐸xHF + 𝑏𝐸xGGA + (1 − 𝑎 − 𝑏)𝐸xLDA + 𝑐𝐸cGGA + (1 − 𝑐)𝐸cLDA .

(3.17)The most famous parameterization is known as B3LYP, which was constructed by using theB88 exchange for the GGA exchange and the LYP correlation for the GGA correlation.
Another form is
1H
𝐸xc= 𝑎𝐸x𝐻𝐹 + (1 − 𝑎)𝐸xDFA + 𝐸cDFA ,(3.19)

where DFA stands for any semilocal density functional approximation (DFA). The most
famous one is called PBE0 where the fitting parameter was set to be 0.25. Hybrid DFT israther free from the self-interaction error. The mixing parameters, 𝑎 to 𝑐, are determinedby comparing with the experimental data for selected matters: Because of its highly empiricalnature, the range of applicability should be limited. It is also noted that hybrid DFT generally
works for insulators but sometimes fails for metals.
Double hybrid DFT: Grimme constructed
2DH
𝐸xc= 𝑎x 𝐸x𝐻𝐹 + (1 − 𝑎x )𝐸xDFA + (1 − 𝑎c )𝐸cDFA + 𝑎c 𝐸cMP2,(3.22)

where MP2 stands for the second order approximation to the correlation energy𝑁1
𝐸cMP2 = − ∑
4

2𝑀

∑

𝑖𝑗=1 𝑎𝑏=𝑁+12|⟨𝜓𝑖 𝜓𝑗 |𝜓𝑎 𝜓𝑏 ⟩ − ⟨𝜓𝑖 𝜓𝑗 |𝜓𝑏 𝜓𝑎 ⟩|,𝜀𝑎 + 𝜀𝑏 − 𝜀𝑖 − 𝜀𝑗(3.23)

where the bracket appearing in the numerator indicates the Coulomb integral⟨𝜓𝑖 𝜓𝑗 |𝜓𝑎 𝜓𝑏 ⟩ = ∬𝜓𝑖∗ (𝒓1 )𝜓𝑗∗ (𝒓𝟐 )𝜓𝑎 (𝒓1 )𝜓𝑏 (𝒓𝟐 )𝑑𝒙1 𝑑𝒙2 .
|𝒓1 − 𝒓𝟐 |(3.24)

Note that 𝑎 and 𝑏 is large than 𝑁 indicating that unoccupied KS orbitals are included in
the formalism contrary to others introduced so far, where only occupied orbitals are included.
This way, the computation of the functional are becoming more and more demanding
although the accuracy is generally improving by that.
Range-separated hybrid DFT: Iikura, Tsuneda, Yanai, and Hirao proposed a long-range
correction (LC) scheme to assign different forms for the exchange depending on the distance
between the electrons aslong range,HF

LC
𝐸xc
= 𝐸x

short range,DFA+ 𝐸x+ 𝐸cDFA .

The range-separation is done using the error function as𝑁𝜎∗ (𝒓 )𝜙 (𝒓 )erf
∗ (𝒓 )𝜙 (𝒓 )𝜙𝑖𝜎(𝜇|𝒓1 − 𝒓2 |)𝜙𝑗𝜎
1
1
𝑗𝜎 𝟏
2
𝑖𝜎 2
long range,HF
𝐸x= − ∑∑∬
𝑑𝒓1 𝑑𝒓2 ,2|𝒓1 − 𝒓2 |𝜎 𝑖𝑗

(3.27)

This is based on the idea that the error of the semilocal approximation to the exchange isserious when two electrons are apart and the error is effectively corrected by the HF exchange.
23

The cutoff distance for the range-separation needs to be introduced empirically.
Exact exchange (HF exchange)
In this scheme, the one-electron orbitals 𝜙𝑖𝜎 appearing in Eq. (3.18) is not treated as anindependent variable but is regarded as a functional of the electron density 𝑛(𝒓). That is, theorbitals are treated as the solution of a KS equation where a (semi-)local form is used for theexchange-correlation functional. The exchange potential thereby generated is the derivativeof Eq. (3.18),
𝛿𝐸x
𝛿𝐸x 𝛿𝑛(𝒓′ )=∫𝑑𝒓′.
𝛿𝑣s (𝒓)𝛿𝑛(𝒓′ ) 𝛿𝑣𝑠 (𝒓)(4.2)

The derivative of the electron density can be obtained by using the definition of the response
function
χ0 (𝒓′ , 𝒓) =𝛿𝑛(𝒓′ )
𝛿𝑣𝑠 (𝒓)as∫ 𝑣x (𝒓′ )𝜒0 (𝒓′, 𝒓)𝑑𝒓′ =𝛿𝐸x.
𝛿𝑣s (𝐫)(4.3)

That is, using the response function, change in the electron density with respect to theeffective potential is expressed using the KS orbital as
𝑁𝜎

χ0

(𝒓′

𝑀

, 𝒓) = − ∑ ∑ ∑ ∬𝜎 𝑖=1 𝑎=𝑁𝜎 +1∗ (𝒓′ )𝜙 (𝒓)𝜙 ∗ (𝒓)𝜙 (𝒓′ )
𝜙𝑖𝜎
𝑎𝜎
𝑎𝜎
𝑖𝜎
𝑑𝒓1 𝑑𝒓2 + c. c.
𝜀𝑎𝜎 − 𝜀𝑖𝜎

(4.5)

When this is used, we can derive the exchange potential using the two-electron integral as
𝑁𝜎

𝛿𝐸x
=∑∑
𝛿𝑣s (𝒓)𝑀∑ (𝜙𝑎𝜎 𝜙𝑗𝜎 |𝜙𝑗𝜎 𝜙𝑖𝜎 )𝜎 𝑖𝑗=1 𝑎=𝑁𝜎 +1∗ (𝒓)𝜙𝑎𝜎 (𝒓)𝜙𝑖𝜎+ c. c.
𝜀𝑎𝜎 − 𝜀𝑖𝜎

(4.6)

This method for expressing the HF exchange is called optimized effective potential (OEP)
method.
Second-order Görling-Levy perturbation theorŷ 𝜆 = 𝑇̂ + 𝜆𝑊
̂ee + V
̂ 𝜆 , isThe Hamiltonian appearing in the adiabatic connection, 𝐻̂ 𝜆 = 𝑉̂s − 𝜆𝑉̂Hx − 𝑉̂c𝜆 asdecomposed using V̂ 𝜆 = (𝑇̂ + 𝑉̂s ) + 𝜆(𝑊̂ee − 𝑉̂Hx ) − 𝑉̂c𝜆 ,𝐻following the second order expansion of 𝑉̂c𝜆 with respect to 𝜆, which starts from the second(2)order as 𝑉̂c𝜆 = 𝜆2 𝑉̂c . Then the perturbation theory is applied to the KS state Φ𝑛=0 as24

|Ψ (1) ⟩ = − ∑ |Φ𝑛 ⟩𝑛≠0̂ee − 𝑉̂Hx|Φ𝑛 ⟩⟨Φ𝑛 |𝑊,ℇ𝑛 − ℇ0and this expression is inserted into Eq. (2.23)
̂ee|Ψ λ ⟩ − ⟨Φ|𝑇̂ + 𝜆𝑊
̂ee|Φ⟩
𝐸c𝜆 = ⟨Ψ λ |𝑇̂ + 𝜆𝑊
to obtain the correlation energy by expanding it in powers of 𝜆. The result is(2)𝐸ĉee|Ψ (1) ⟩ = ⟨Φ|𝑊̂ee − 𝑉̂Hx |Ψ (1) ⟩ = − ∑ |= ⟨Φ|𝑊𝑛≠0̂ee − 𝑉̂Hx|Φ𝑛 ⟩|2|⟨Φ|𝑊.
ℇ𝑛 − ℇ0

It was found this does not improve the correlation satisfactorily, indicating that it isnecessary to go beyond the second order approximation.
There are so many types of approximation. As we have seen above, some of them are
introduced to approximate the HF exchange, which takes large computational time, with semi
local one, for which the computational time is generally much smaller. Others are introduced
to better describe the correlation functionals semilocally. In many cases, the accuracy has been
improved empirically in the sense that fitting parameters are introduced to better reproduce
experiments or accurate simulations. This idea is recently combined with machine learning
methods which enable to use much flexible fitting functions like neural network. Contrary tothis approach, there have been trials to derive plausible functionals theoretically.
Fluctuation dissipation formula
To derive the functional form theoretically, I will introduce the general idea of many-bodytheories called fluctuation dissipation formula that relates the correlation described using thepair density 𝑛2 (𝒓1 , 𝒓2 ) with the fluctuation ⟨Ψ|𝑛̂2 (𝒓1 , 𝒓2 )|Ψ⟩ described by the expectation
value of the pair density operator. The fluctuation is known to be more easily evaluated thanthe correlation when the fluctuation is described by using the response functions. Note that
this is a step toward combining the density functional theory with the many-body theory.
Density operator
To compactly describe the fluctuation dissipation theorem, let us introduce the density
operator without using the terminology of the field theory. The density operator isdefined as𝑁𝑛̂1 (𝒓) = ∑ 𝛿(𝒓 − 𝒓𝑖 )𝑖where 𝒓𝑖 is the coordinate of the constituent 𝑁 electrons and 𝒓 is a point in real-space.
The density can be obtained as𝑛(𝒓) = ⟨Ψ|𝑛̂1 (𝒓)|Ψ⟩.
25

The pair density operator can also be defined as𝑛̂2 (𝒓1 , 𝒓2 ) = 𝑛̂1 (𝒓1 )𝑛̂1 (𝒓2) − 𝑛̂1 (𝒓1 )𝛿(𝒓1 − 𝒓2 ).
Those are used to obtain1
𝑇 = − ∫[∇2𝒓 𝑛̂1 (𝒓, 𝒓′ )]𝒓′ =𝒓 𝑑𝒓
2
1
𝑉ee = ∬ 𝑤ee (𝒓1 , 𝒓2 )𝑛̂2 (𝒓1 , 𝒓2 )𝑑𝒓1 𝑑𝒓2
2
𝑉ne = ∫ 𝑣ne (𝒓)𝑛̂(𝒓)𝑑𝒓

To obtain the correlation energy using the adiabatic connection,1
1 1
𝑛 𝜆 (𝒓 , 𝒓 )̂ee|Ψ λ ⟩ − ⟨Φ|𝑊
̂ee|Φ⟩ = ∫ 𝑑𝜆 ∬ 2,c 1 2 𝑑𝒓1 𝑑𝒓2 ,𝐸c = ∫ 𝑑𝜆 ⟨Ψ λ |𝑊|𝒓1 − 𝒓2 |2 0
0

(4.15)

we need the quantity
𝜆 (𝒓
𝜆
𝑛2,c
1 , 𝒓2 ) = 𝑛2 (𝒓1 , 𝒓2 ) − 𝑛2,KS (𝒓1 , 𝒓2 ).

Using the relations𝑛2𝜆 (𝒓1 , 𝒓2 ) = ⟨Ψ λ|𝑛̂2 (𝒓1 , 𝒓2 )|Ψ λ ⟩ = ⟨Ψ λ |𝑛̂1 (𝒓1 )𝑛̂1 (𝒓2 )|Ψ λ ⟩ − 𝛿(𝒓1 − 𝒓2 )⟨Ψ λ |𝑛̂1 (𝒓1 )|Ψ λ ⟩𝑛2,KS (𝒓1 , 𝒓2 ) = ⟨Φ|𝑛̂1 (𝒓1 )𝑛̂1 (𝒓2 )|Φ⟩ − 𝛿(𝒓1 − 𝒓2 )⟨Φ|𝑛̂1 (𝒓1 )|Φ⟩

(4.16 − 4.17)

and the fact that the density is independent of 𝜆 (by definition), we have
𝜆 (𝒓
λ
𝑛2,c
̂1 (𝒓1 )𝑛̂1 (𝒓2 )|Ψ λ ⟩ − ⟨Φ|𝑛̂1 (𝒓1 )𝑛̂1 (𝒓2 )|Φ⟩.
1 , 𝒓2 ) = ⟨Ψ |𝑛(4.18)

Now we relate the correlation function to the response function, i.e. response of the density
to the external perturbation. Equation (4.5) is a static response function
𝛿𝑛(𝒓′ )
𝛿𝑣𝑠 (𝒓)

χ0 (𝒓′ , 𝒓) =but it is more convenient to use a time-dependent one like
χ0 (𝒓′ 𝑡′, 𝒓𝑡) =𝛿𝑛(𝒓′ 𝑡′)𝛿𝑣𝑠 (𝒓𝑡)and to consider the static limit. We need to remind that eq. (4.5) describes a response of the
non-interacting system. What we need is the response of the interacting systems
𝜒(𝒓′ 𝑡′, 𝒓𝑡) =𝛿𝑛(𝒓′ 𝑡′).
𝛿𝑣ne (𝒓𝑡)We must relate 𝜒 with χ0. For this purpose, we use a many-body Green’s function formalism
for the time-dependent Schrödinger equation: Afterwards we will return to the timedependent DFT formalism. So, let us introduce the many-body Green’s function formalism
for a while and then return to the discussion on the fluctuation dissipation theorem.
26

Many-body Green’s function (article 2 of the references shown above)
We briefly introduce the many-body Green’s function following another article written by
Toulouse. For this purpose, we first change the language to that of the field theory. The
Hamiltonian can be rewritten as11
̂ = ∫ 𝑑𝒙 𝜓̂ † (𝒙) (− ∇2 + 𝑣ne (𝒓)) 𝜓̂(𝑥) + ∬ 𝑑𝒙1 𝑑𝒙2 𝜓̂ † (𝒙1 )𝜓̂ † (𝒙2 )𝑣(𝒓1 , 𝒓2 )𝜓̂(𝒙2 )𝜓̂(𝒙1),𝐻22where 𝑣(𝒓1 , 𝒓2 ) is the Coulombic interaction 1/|𝒓1 − 𝒓2 | and the term parenthesized in the1

first term − 2 ∇2 + 𝑣ne (𝒓) will be denoted as ℎ(𝒓). Here we have introduced the field operator
𝜓̂(𝒙) ≡ ∑ 𝜓𝑖 (𝒙)𝑐̂𝑖 ,iwhich will be time-evolved using the Heisenberg representation as𝜓̂(1) ≡ 𝜓̂(𝒙1 , 𝑡1 ) = 𝑒 𝑖𝐻̂𝑡1 𝜓(𝒙1 )𝑒 −𝑖𝐻̂𝑡1 .
This is just a rewrite of equations but will help to make equations much simpler. I intuitively
note that the “physical” Green’s function is the retarded one but “convenient” Green’s
function is the time-ordered one; there is a way to make the former from the latter. The timeordered Green’s function is given for an 𝑁 electron system byi𝐺(1,2) = ⟨𝑁|𝑇[𝜓̂(1)𝜓̂ † (2)]|𝑁⟩⟨𝑁|𝜓̂(1)𝜓̂ † (2)|𝑁⟩ when t1 > 𝑡2={−⟨𝑁|𝜓̂ † (2)𝜓̂(1)|𝑁⟩ when t1 < 𝑡2(2)

Here operators are ordered from the more past ones to the more recent ones. The Green’s
function follows the equation of motion(i𝜕 ∇2+ ) 𝐺(1,2) = 𝛿(1 − 2) − 𝑖 ∑ 𝑣(3,1) ⟨𝑁|𝑇[𝜓̂ † (3)𝜓̂(3)𝜓̂(1)𝜓̂† (2)]|𝑁⟩.
𝜕𝑡 2
𝜎3This can be derived as follows:First, consider the equation of motion of the field operator, which is𝜕𝜓̂(1)
̂ ]𝑒 −i𝐻̂𝑡1 .
= −i𝑒 i𝐻̂𝑡1 [𝜓̂(𝒙1 ), 𝐻𝜕𝑡1
By calculating the commutation relation, we get̂ ] = ℎ(𝒓1 )𝜓̂(𝒙1 ) + ∫ 𝑑𝒙2 𝜓̂ † (𝒙2 )𝑣(𝒓1 , 𝒓2 )𝜓̂(𝒙2)𝜓(𝒙1 )[𝜓̂(𝒙1 ), 𝐻
and then𝜕𝜓̂(1)= −iℎ(1)𝜓̂(1) − 𝑖 ∫ 𝑑3𝜓̂(3)𝑣(1,3)𝜓̂(3)𝜓̂(1).
𝜕𝑡1This can be used to take the time derivative of the one-body Green’s function as27𝜕𝜕𝜓̂(1) ̂ †𝐺(1,2) = 𝛿(𝑡1 − 𝑡2 )⟨𝑁|𝜓̂(1)𝜓̂† (2)|𝑁⟩ + 𝜃(𝑡1 − 𝑡2 ) ⟨𝑁|𝜓 (2)|𝑁⟩
𝜕𝑡1
𝜕𝑡1

i

𝜕𝜓̂(1)− 𝛿(𝑡1 − 𝑡2 )⟨𝑁|𝜓̂ † (2)𝜓̂(1)|𝑁⟩ − 𝜃(𝑡1 − 𝑡2 ) ⟨𝑁|𝜓̂ † (2)|𝑁⟩
𝜕𝑡1This can be simplified asi𝜕𝐺(1,2) = 𝛿(1,2) + ℎ(1)𝐺(1,2)
𝜕𝑡1− 𝑖 ∫ 𝑑3 𝑣(1,3)[𝜃(𝑡1 − 𝑡2 )⟨𝑁|𝜓̂ † (3)𝜓̂(3)𝜓̂(1)𝜓̂ † (2)|𝑁⟩ − 𝜃(𝑡2
− 𝑡1 ) ⟨𝑁|𝜓̂ † (2)𝜓̂† (3)𝜓̂(3)𝜓̂(1)|𝑁⟩

Using the definition of the two-body Green’s function as will appear below
i2 𝐺2 (1,2; 1′ , 2′) = ⟨𝑁|𝑇[𝜓̂(1)𝜓̂(2)𝜓̂† (1′)𝜓̂† (2′)]|𝑁⟩,further simplification is possible and the result isi𝜕
𝐺(1,2) = 𝛿(1,2) + ℎ(1)𝐺(1,2) + i ∫ 𝑑3𝑣(1,3)𝐺2 (1,3+ ; 2,3++ ),𝜕𝑡1

where time is ordered as 𝑡3++ > 𝑡3+ > 𝑡3 and the difference in the time is infinitesimal.
The one-body density matrix can be written as
𝑛1 (𝒙1, 𝒙2 ) = ⟨𝑁|𝜓̂ † (𝒙2)𝜓̂(𝒙1 )|𝑁⟩ = −i𝐺(𝒙1 𝑡1 , 𝒙2 𝑡1+ ).

(3)

and the one-body electron density as
𝑛1 (𝒓1 , 𝒓2 ) = −i ∑ 𝐺(𝒓1 𝜎𝑡1 , 𝒓2 𝜎𝑡1+ ) .

(4)σThe Green’s function can be rewritten using the completeness relation of 𝑁 ± 1 body statesasi𝐺(𝒙1 , 𝒙2 , 𝜏) = 𝜃(𝜏) ∑⟨𝑁|𝜓̂(𝒙1)|𝑁 + 1, 𝑎⟩⟨𝑁 + 1, 𝑎|𝜓̂ † (𝒙2)|𝑁⟩𝑒 −i(𝐸𝑁+1,𝑎−𝐸𝑁 )𝜏
𝑎− 𝜃(−𝜏) ∑⟨𝑁|𝜓̂ † (𝒙2)|𝑁 − 1, 𝑎⟩⟨𝑁 − 1, 𝑎|𝜓̂(𝒙1 )|𝑁⟩𝑒 −i(𝐸𝑁 −𝐸𝑁−1,𝑎 )𝜏 , (5)𝑎so that the Fourier component is
𝐺(𝒙1 , 𝒙2 , 𝜔) = ∑𝑎⟨𝑁|𝜓̂(𝒙1 )|𝑁 + 1, 𝑎⟩⟨𝑁 + 1, 𝑎|𝜓̂ † (𝒙2 )|𝑁⟩
𝜔 − 𝐸𝑁+1,𝑎 + 𝐸𝑁 + 𝑖0++∑𝑎⟨𝑁|𝜓̂ † (𝒙2 )|𝑁 − 1, 𝑎⟩⟨𝑁 − 1, 𝑎|𝜓̂(𝒙1 )|𝑁⟩.
𝜔 − 𝐸𝑁−1,𝑎 + 𝐸𝑁 − 𝑖0+

This is called the Lehman representation of the one-particle Green’s function. It is simply
rewritten as

28

=∑
𝑎

𝑓𝑎 (𝒙1 )𝑓𝑎∗ (𝒙2 )𝑓𝑖 (𝒙1 )𝑓𝑖∗ (𝒙2 )+∑,𝜔 + 𝐴𝑎 + 𝑖0+𝜔 − 𝐼𝑖 − 𝑖0+(6)𝑖where the affinity and ionization potential has been used. The function 𝑓 appearing in the
denominator is called as the (quasi particle) wave function.
Linear response and response function
At this moment, let us definite the linear-response function. For this, we will considera system applied with a perturbation likê ′ (t) = ∫ 𝑛̂(𝒓)𝑉ext (𝒓, 𝑡)𝑑𝒓𝐻and evaluate the expectation value of the density. We can use the equation of motion for
the time-evolution operator in the interaction representationi∂𝑈𝐼 (𝑡, 𝑡0 )
̂ ′𝐼 (𝑡)𝑈𝐼 (𝑡, 𝑡0 )=𝐻𝜕𝑡which is solved as
t

t

t0

t0

′

̂ ′𝐼 (𝑡 ′ ) = 1 − i ∫ 𝑑𝑡 ′ 𝑒 i𝐻̂0 𝑡 𝐻
̂ (𝑡 ′ )𝑒 −i𝐻̂0𝑡′ .
𝑈𝐼 (𝑡, 𝑡0 ) ≃ 1 − i ∫ 𝑑𝑡 ′ 𝐻
Therefore, the expectation value of the density is given as𝛿𝑛(𝒓𝑡) = 〈𝑈𝐼† (𝑡, 𝑡0 )𝑛̂𝐼 (𝒓𝑡)𝑈𝐼 (𝑡, 𝑡0 )〉.
The terms linear with respect to the perturbation is𝑡′𝛿𝑛(𝒓𝑡) = 𝑖 ∫ 𝑑𝑡 ′ ∫ 𝑑𝒓 〈[𝑛̂𝐼 (𝒓𝑡), 𝑛̂𝐼 (𝒓′ 𝑡 ′ )]〉𝑉ext (𝒓′, 𝑡′)𝑡0

The response function is thus given by the density-density correlation function as𝜒(𝒓𝑡, 𝒓′ 𝑡 ′ ) = i〈[𝑛̂𝐼 (𝒓𝑡), 𝑛̂𝐼 (𝒓′ 𝑡 ′ )]〉.
Here the subscript 𝐼 may be regarded as the Heisenberg representation. Depending on
the boundary condition, we can make the causal or time-ordered response function. Notethat this is a plausible extension of the static definition 𝑛̂2 (𝒓1 , 𝒓2 ) = 𝑛̂1 (𝒓1 )𝑛̂1 (𝒓2) −𝑛̂1 (𝒓1 )𝛿(𝒓1 − 𝒓2 )One can likewise introduce the two-body Green’s function asi2 𝐺2 (1,2; 1′ , 2′) = ⟨𝑁|𝑇[𝜓̂(1)𝜓̂(2)𝜓̂ † (1′)𝜓̂ † (2′)]|𝑁⟩,(7)

which describes the propagation of a pair of particles; electrons, holes, or electron-hole. A
related quantity called four-point linear response function is defined asi𝜒(1,2; 1′ , 2′) = i2 𝐺2 (1,2; 1′ , 2′ ) − 𝑖𝐺(1,1′ )𝑖𝐺(2,2′ ).
For a special time-ordering relation, it can be rewritten as

29

(8)

i𝜒(1,2; 𝒙′1 𝑡1+ , 𝒙′2 𝑡2+ )= ⟨𝑁|𝑇[𝜓̂ † (1′)𝜓̂(1)𝜓̂ † (2′)𝜓̂(2)]|𝑁⟩− ⟨𝑁|𝜓̂ † (1′ )𝜓̂(1)|𝑁⟩⟨𝑁|𝜓̂ † (2′ )𝜓̂(2)|𝑁⟩,(9)and the Fourie transformed equation is then given by𝜒(𝒙1, 𝒙2 ; 𝒙1′ , 𝒙′2 ; 𝜔)= ∑[⟨𝑁|𝜓̂ † (1′)𝜓̂(1)|𝑁, 𝑛⟩⟨𝑁, 𝑛|𝜓̂ † (2′ )𝜓̂(2)|𝑁⟩𝜔 − (𝐸𝑁,𝑛 − 𝐸𝑁 ) + 𝑖0+𝑛≠0−⟨𝑁|𝜓̂ † (2′ )𝜓̂(2)|𝑁, 𝑛⟩⟨𝑁, 𝑛|𝜓̂ † (1′ )𝜓̂(1)|𝑁⟩𝜔 + (𝐸𝑁,𝑛 − 𝐸𝑁 ) − 𝑖0+

]

(11)

This indicates that the four-point linear response function in this special time-orderingcorresponds to the linear-response of particle-hole pairs.
One can further introduce the independent-particle four-point linear-response function𝜒IP (1,2; 1′ , 2′) = −𝑖𝐺(1,1′)𝐺(2,2′ )(12)and the corresponding response function𝜒IP (𝒙1 , 𝒙2 ; 𝒙1′ , 𝒙′2 ; 𝜏) = −𝑖𝐺(𝒙1 , 𝒙′2 ; 𝜏)𝐺(𝒙2, 𝒙1′ ; −𝜏)𝜒IP (𝒙1 , 𝒙2 ; 𝒙1′ , 𝒙′2 ; 𝜔) = ∑ [
𝑖≤𝑁<𝑎

(13)

𝑓𝑖∗ (𝒙1′)𝑓𝑎 (𝒙1 )𝑓𝑎∗ (𝒙′2 )𝑓𝑖 (𝒙2 ) 𝑓𝑖∗ (𝒙′2 )𝑓𝑎 (𝒙2 )𝑓𝑎∗ (𝒙1′ )𝑓𝑖 (𝒙1)−
] . (14)𝜔 − (ℇ𝑎 − ℇ𝑖 ) + 𝑖0+𝜔 + (ℇ𝑎 − ℇ𝑖 ) − 𝑖0+

Having introduced the many-body Green’s function, let us return to the discussion on thefluctuation dissipation theorem. The response function (Eqs. (9-10)) can be written asiχλ (𝒓1 𝑡1 , 𝒓2 𝑡2 ) = ⟨Ψ 𝜆 |𝑇[𝑛̂ 𝜆 (𝒓1 𝑡1 )𝑛̂𝜆 (𝒓2𝑡2)]|Ψ 𝜆 ⟩ − ⟨Ψ 𝜆 |𝑇[𝑛̂ 𝜆 (𝒓1𝑡1 )]|Ψ 𝜆 ⟩⟨Ψ 𝜆 |𝑇[𝑛̂ 𝜆 (𝒓2 𝑡2 )]|Ψ 𝜆 ⟩= ⟨Ψ 𝜆 |𝑇[𝑛̂ 𝜆 (𝒓1 𝑡1 )𝑛̂ 𝜆 (𝒓2 𝑡2 )]|Ψ 𝜆 ⟩ − ⟨Ψ 𝜆 |𝑛̂(𝒓1 )|Ψ 𝜆 ⟩⟨Ψ 𝜆 |𝑛̂(𝒓2 )|Ψ 𝜆 ⟩(4.19)+

If we set 𝑡2 = 𝑡1 + 0 , the Heisenberg phase factor is canceled out andiχλ (𝒓1 , 𝒓2 ; 𝜏 = 0− ) = ⟨Ψ 𝜆 |𝑛̂(𝒓1 )𝑛̂(𝒓2 )|Ψ 𝜆 ⟩ − ⟨Ψ 𝜆 |𝑛̂(𝒓1 )|Ψ 𝜆 ⟩⟨Ψ 𝜆 |𝑛̂(𝒓2 )|Ψ 𝜆 ⟩.

(4.20)

It is possible to derive similar equation for λ = 0, givingiχ0 (𝒓1, 𝒓2 ; 𝜏 = 0− ) = ⟨Φ|𝑛̂(𝒓1 )𝑛̂(𝒓2 )|Φ⟩ − ⟨Φ|𝑛̂(𝒓1 )|Φ⟩⟨Φ|𝑛̂(𝒓2 )|Φ⟩.

(4.21)

This makes it possible to write the two-body density correlation in terms of the response
function as
𝜆 (𝒓
−
−
𝑛2,c
1 , 𝒓2 ) = i[𝜒𝜆 (𝒓1 , 𝒓2 ; 𝜏 = 0 ) − 𝜒0 (𝒓1 , 𝒓2 ; 𝜏 = 0 )],(4.22)

and using the Fourier transformed form as
∞
𝜆 (𝒓
𝑛2,c
1 , 𝒓2 ) = −∫−∞𝑑𝜔 𝑖𝜔0+[𝜒𝜆 (𝒓1 , 𝒓2 ; 𝜔) − 𝜒0 (𝒓1 , 𝒓2 ; 𝜔)],𝑒
2𝜋i

(4.23)

which is known as the fluctuation dissipation theorem. When this is combined with the
adiabatic connection formula, we get
30

1

∞

0

−∞

1
𝑑𝜔 𝑖𝜔0+𝜒𝜆 (𝒓1 , 𝒓2 ; 𝜔) − 𝜒0 (𝒓1 , 𝒓2 ; 𝜔)𝐸c = − ∫ 𝑑𝜆 ∫𝑒
𝑑𝒓1 𝑑𝒓2∬|𝒓1 − 𝒓2 |22𝜋iThis is a particularly important formula which provides a practical way to compute thecorrelation functional without the many-body wave function.
The problem is now reduced to the evaluation of the response function 𝜒𝜆 . This can bedone either by using the time-dependent DFT or the many-body perturbation within the GW
approximation. For the moment, let us concentrate the former.
Time-dependent density functional theory (TD-DFT)To explain the response function within DFT, I will briefly introduce TDDFT. According
to Runge and Gross, there exists a one-to-one correspondence between the time-dependent
density 𝑛(𝒓, 𝑡) and the time-dependent effective potential 𝑣s (𝒓, 𝑡), so that 𝑛(𝒓, 𝑡) can beobtained by solving the time-dependent Kohn-Sham equation:𝑖𝜕1
𝜑𝑖 (𝒓, 𝑡) = [− ∇2 + 𝑣s (𝒓, 𝑡)] 𝜑𝑖 (𝒓, 𝑡),
𝜕𝑡
2

and
𝑁

𝑛(𝒓, 𝑡) = ∑|𝜑𝑖 (𝒓, 𝑡)|2𝑖This is a natural extension of the Kohn-Sham DFT.
According to the many-body theory, in the linear response regime, the change in thedensity and change in the external potential can be related via the response function as∞𝜆𝛿𝑛𝜆 (𝒓, 𝑡) = ∫ 𝑑𝑡 ′ ∫ 𝑑𝒓′ 𝜒𝜆 (𝒓𝑡, 𝒓′ 𝑡 ′ )𝛿𝑣ext(𝒓′ 𝑡 ′ )0where𝜒𝜆 (𝒓𝑡, 𝒓′ 𝑡 ′ ) =𝛿𝑛(𝒓𝑡)𝜆 (𝒓′ ′ )𝛿𝑣ext𝑡

the response function can be given (from the definition) as𝜒𝜆 (𝒓, 𝒓′ , 𝜔) = ∑ [
𝑚≠0

𝜆
𝜆
𝜆
𝜆
|𝑛̂(𝒓′ )|Ψ 𝜆 ⟩ ⟨Ψ 𝜆 |𝑛̂(𝒓′ )|Ψ𝑚|𝑛̂(𝒓)|Ψ 𝜆 ⟩⟨Ψ 𝜆 |𝑛̂(𝒓)|Ψ𝑚
⟩⟨Ψ𝑚
⟩⟨Ψ𝑚
−
].
𝜔 − 𝜔𝑚 + 𝑖0+𝜔 + 𝜔𝑚 + 𝑖0+

In TDDFT, the equations are replaced by those of a noninteracting system∞δ𝑛(𝒓, 𝑡) = ∫ 𝑑𝑡 ′ ∫ 𝑑𝑟 ′ 𝜒𝜆,KS (𝒓𝑡, 𝒓′ 𝑡 ′ )𝛿𝑣s𝜆 (𝒓′ 𝑡 ′ ),031

𝜒𝜆,KS (𝒓𝑡, 𝒓′ 𝑡 ′ ) =δ𝑛(𝒓𝑡),
δ𝑣s𝜆 (𝒓′ 𝑡 ′ )and𝜒𝜆,KS (𝒓, 𝒓′ , 𝜔) = ∑The effective potential,

𝜑𝑘𝜆∗ (𝒓)𝜑𝑗𝜆 (𝒓)𝜑𝑗𝜆∗ (𝒓′ )𝜑𝑘𝜆 (𝒓′)𝜔 − 𝜔𝑗𝑘 + 𝑖𝜂𝑘𝑗
𝜆 (𝒓′ ′ )𝜆 (𝒓′ ′ )𝜆 (𝒓′ ′ ),𝑣s𝑡 = 𝑣ext 𝑡 + 𝑣H (𝒓′ 𝑡) + 𝑣xc
𝑡 is related to the density

via the exchange-correlation kernel 𝐾 as
𝜆 (𝒓𝑡, ′ ′ )
𝐾xc
𝒓𝑡 =

𝜆 (𝒓𝑡)
𝛿𝑣xc,
𝛿𝑛(𝒓′ 𝑡 ′ )so that𝜆 (𝒓𝑡)
δ𝑣s𝜆 (𝒓𝑡) = 𝛿𝑣ext+ ∫ d𝐫′𝛿𝑛(𝒓’𝑡)𝜆+ ∫ 𝑑𝒓′ 𝑑𝑡 ′ K xc
(𝒓𝑡, 𝐫 ′ t ′ )𝛿𝑛(𝒓′ 𝑡 ′ )|𝐫 − 𝐫 ′ |Since the same density response results in both cases, we get𝜒𝜆 (𝒓𝑡, 𝒓′ 𝑡 ′ ) = 𝜒λ,KS (𝒓𝑡, 𝒓′ 𝑡 ′ ) + ∭ 𝜒𝜆,KS (𝒓𝑡, 𝒓1 𝑡1 ) [𝛿(𝑡1 − 𝑡2 )
𝜆 (𝒓
′ ′
+ 𝐾xc
1 𝑡1 , 𝒓2 𝑡2 )] 𝜒𝜆 (𝒓2 𝑡2 , 𝒓 𝑡 )|𝒓1 − 𝒓2 |Namely, once functional form for the exchange-correlation kernel is determined, we get the
response function, with which one can calculate the correlation functional. The problem has
𝜆not been solved by this, instead, the problem has been transferred to the construction of 𝐾xc,

which is practically easier. Also, there is a way to construct it from a rigorous many-bodytheory.
Let us then introduce the many-body theory for constructing the response function. Again,
we will follow the article written by Toulouse.
Self-energy
Returning to the equation of motion, which was derived above asi𝜕𝐺(1,2) = 𝛿(1,2) + ℎ(1)𝐺(1,2) + i ∫ 𝑑3𝑣(1,3)𝐺2 (1,3+ ; 2,3++ ),
𝜕𝑡1

(15)

let us introduce the self-energy operator. For this purpose, let us assume that we can find asolution for the following equation∫ 𝑑3ΣHxc (1,3)𝐺(3,2) = −𝑖 ∫ 𝑑3𝑣(1,3)𝐺2 (1,3+ ; 2,3++ )(17)

or equivalently let us introduce the following quantity
ΣHxc (1,4) = −𝑖 ∫ 𝑑3𝑣(1,3)𝐺2 (1,3+ ; 2,3++ )𝐺 −1 (2,4).
Then, we obtain a formally closed form for the Green’s function
32

(18)

[𝑖

𝜕
− ℎ(1)] 𝐺(1,2) = 𝛿(1,2) + ∫ 𝑑3ΣHxc (1,3)𝐺(3,2).
𝜕𝑡1

(16)

This means that with the self-energy (Eq.(18)) we do not have to introduce the complicatedmulti-body Green’s function.
Note that Green’s function is given for independent electron approximation as[𝑖𝜕− ℎ(1)] 𝐺indep (1,2) = 𝛿(1,2)
𝜕𝑡1

(19)

−1 (4,3)𝐺(3,2),Using 𝐺(1,2) = ∫ 𝐺indep (1,4)𝐺indep
we get−1 (1,3)− ΣHxc (1,3)]G(3,2) = 𝛿(1,2)∫ 𝑑3[𝐺indep(20)

𝐺(1,2) = 𝐺indep (1,2) + ∬ 𝐺indep (1,3) ΣHxc (3,4)𝐺(4,2),(21)−1 (1,2)
𝐺 −1(1,2) = 𝐺indep− ΣHxc (3,4).

(22)

This is rewritten as

or as
Equation (21) is known as the Dyson equation for the one-body Green’s function.
Fourie transform is (assuming homogeneity of the system with respect to time) given as[𝜔 − ℎ(𝒓1 )]𝐺(𝒙1 , 𝒙𝟐 ; 𝜔) = 𝛿(𝒙1 , 𝒙2 ) + ∫ 𝑑𝒙3 ΣHxc (𝒙1 , 𝒙3 ; 𝜔)𝐺(𝒙3 , 𝒙2 ; 𝜔)(23)

Using the pole ℇ𝑘 and the wave function 𝑓𝑘 (𝒙) of the Green’s function, we can rewrite Eq.
(23) at the pole as[ℇ𝑘 − ℎ(𝒓1 )]𝑓𝑘 (𝒙1 )𝑓𝑘∗ (𝒙2) = ∫ 𝑑𝒙3 ΣHxc (𝒙1 , 𝒙3 ; ℇ𝑘 )𝑓𝑘 (𝒙3 )𝑓𝑘∗ (𝒙2) = 0,(24)when there is no degeneracy. This can be further simplified asℎ(𝒓1 )𝑓𝑘 (𝒙1 ) + ∫ 𝑑𝒙3 ΣHxc (𝒙1 , 𝒙3 ; ℇ𝑘 )𝑓𝑘 (𝒙3 ) = ℇ𝑘 𝑓𝑘 (𝒙1 ).

(25)

If the imaginary part of the self-energy is zero, Eq. (25) is equivalent to a one-particle
Schrödinger equation. Because of the imaginary part, however, it becomes similar to theequation with absorption potential, indicating decreasing amplitude of the wave function withtime. It is known in some cases that the imaginary part is small for the energy states close to
the Fermi level. So, Eq.(25) is called quasiparticle equation.
Note that the self-energy is apparently equivalent to the exchange-correlation potential of
TD-DFT. It is therefore possible, in principle, to derive the exchange-correlation from the
many-body theory although it is determined from empirical argument.
Exchange-correlation kernel from first principles
33

1. Green’s function using the interaction representation
Suppose the system is perturbed by an external field 𝑢(𝒓), or
𝐻 ′ (𝑡1 ) = ∫ 𝑑1 𝑢(1)𝑛(1).
The Green’s function can be described using the interaction representation
𝐺(1,1′) = −i

⟨𝑁|𝑇[𝑆𝜓̂(1)𝜓̂† (1′)]|𝑁⟩,⟨𝑁|𝑇[𝑆]|𝑁⟩

where 𝑆 is the time-evolution operator in the interaction representation𝑡𝑎

𝑆(𝑡𝑎 , 𝑡𝑏 ) = exp [−i ∫ 𝑑𝑡1 𝐻 ′ 𝐼 (1)] ;S = S(−∞, ∞).

𝑡𝑏

The two-body Green’s function is also
𝐺2 (1,1′ ; 2,2′) = (−i)2⟨𝑁|𝑇[𝑆𝜓̂(1)𝜓̂(2)𝜓̂† (2′)𝜓̂ † (1′)]|𝑁⟩,⟨𝑁|𝑇[𝑆]|𝑁⟩similarly to the one in the limit of 𝑢 → 0.
2. Equation of motion
The Green’s function follows an equation of motion[i𝜕− ℎ(1) − 𝑢(1)] 𝐺(1,1′) + 𝑖 ∫ 𝑑2𝑣(1,2)𝐺2 (1,1′ , 2+ , 2++ ) = 𝛿(1,1′ ).
𝜕𝑡1When the Green’s function is differentiated with respect to 𝑢, we get the equation known as
Schwinger’s trick𝛿𝐺(1,1′ )= −𝐺2 (1,1′ , 2,2+ ) + 𝑖𝐺(1,1′)𝐺(2,2+ ) = −𝐺2 (1,1′ , 2,2+ ) + 𝑖𝐺(1,1′ )〈𝑛̂(2)〉.
𝛿𝑢(2)With this, one can formally remove the two-body Green’s function from the equation of
motion:[i𝜕𝛿𝐺(1,1′ )− ℎ(1) − 𝑢(1) − ∫ 𝑑2𝑣(1,2)〈𝑛̂(2)〉] 𝐺(1,1′) − 𝑖 ∫ 𝑑2𝑣(1+ , 2)
= 𝛿(1,1′ ).
𝜕𝑡1
𝛿𝑢(2)

As we did above, we can define the self-energy operator
𝛴(1,1′ ) = ∫ 𝑑2𝑣(1,2)〈𝑛̂(2)〉 𝛿(1,1′ ) + 𝑖 ∬ 𝑑23 𝑣(1+ , 2)

𝛿𝐺(1,3) −1
𝐺 (3,1′ )𝛿𝑢(2)≡ 𝛴H (1,1′ ) + 𝛴xc (1,1′ )so that the equation of motion can be rewritten as[i𝜕− ℎ(1) − 𝑢(1)] 𝐺(1,1′) − 𝑖 ∫ 𝑑2𝛴(1,2)𝐺(2,1′ ) = 𝛿(1,1′ ).
𝜕𝑡1When unperturbed Green’s function is defined from

34

[i

𝜕
− ℎ(1)] 𝐺 (0) (1,1′ ) = 𝛿(1,1′ ),𝜕𝑡1

we have the equation of Dyson type𝐺(1,1′ ) = 𝐺 (0) (1,1′ ) + ∬ 𝑑23 𝐺 (0) (1,2)𝛴(2,3)𝐺(3,1′ ).

3. Self-energy
The self-energy operator has a form difficult to understand. Let us relate it physicallyunderstandable quantities. For this purpose, we define the total potential as a sum of theexternal potential and the Hartree-potential
𝑉(1) = 𝑢(1) + ∫ 𝑑2𝑣(1,2)〈𝑛̂(2)〉and consider the response of the Green’s function to the change in the total potential;−

δG−1 (1,2)↦ 𝛤(1,2; 3)𝛿𝑉(3)The vertex function 𝛤(1,2; 3) can be rewritten in a Dyson form asΓ(1,2; 3) = 𝛿(1,2)𝛿(2,3) + ∭ 𝑑4567

𝛿Σ xc (1,2)
G(4,6)G(7,5)Γ(6,7; 3),
𝛿𝐺(4,5)so that one may start the first order approximation from 𝛿(1,2)𝛿(2,3).
With the vertex, let us remove the unpleasant δ/δu operator from the definition of theself-energy. The exchange-correlation part can be rewritten as𝛴xc (1,2) = i ∬ 𝑑3456 𝑣(1+ , 3)

𝛿𝑉(4)𝐺(1,6)𝛤(6,2; 5).
𝛿𝑢(3)Using the fact that the change in the total potential induced by the external potential is thedielectric function
ϵ−1 (4,3) =

𝛿𝑉(4)
𝛿𝑢(3)

and the fact that the dielectric function operated to the Coulomb potential is the screened
Coulomb
𝑊(1+ , 4) = ∫ 𝑑3 𝑣(1+ , 3)ϵ−1 (4,3)Σ xc can be rewritten asΣ xc (1,2) = i ∬ 𝑑56 𝑊(1+ , 3)𝐺(1,6)Γ(6,2; 5).
To remove the δ/δu operator from the screened Coulomb, it is rewritten as35

𝑊(1+ , 2) = 𝑣(1+ , 2) + ∫ 𝑑34 𝑣(1+ , 4)𝛿𝐺(4,4+ ) 𝛿𝑉(5)𝑣(3,2)
𝛿𝑉(5) 𝛿𝑢(3)= 𝑣(1+ , 2) + ∫ 𝑑34 𝑣(1+ , 4)𝛿𝐺(4,4+ )𝑊(5,2)
𝛿𝑉(5)= 𝑣(1+ , 2) + ∫ 𝑑34 𝑣(1+ , 4)

𝛿〈𝜌̂(2)〉𝑊(5,2)
𝛿𝑉(5)and then, the derivative appearing in the last line, which is the polarization, is rewritten using
the vertex function as𝑃(1,2) = i ∫ 𝑑34 𝐺(1,3)

𝛿𝐺 −1 (3,4)
𝐺(4,1+ ) = −𝑖 ∫ 𝑑34 𝐺(1,3)𝐺(4,1+ )Γ(3,4; 2).
𝛿𝑉(2)Then, we can obtain a Dyson equation for the screened Coulomb as
𝑊(1+ , 2) = 𝑣(1+ , 2) + ∫ 𝑑34 𝑣(1+ , 4)𝑃(2,5)𝑊(5,2).
4. Hedin’s equation
Finally, we arrive at a closed set of equations called Hedin’s equation
𝐺(1,2) = 𝐺 0 (1,2) + ∫ 𝑑34𝐺 0 (1,3)Σ(3,4)𝐺(4,2)Σ xc (1,2) = 𝑖∫ 𝑑34𝑊(1+ , 3)𝐺(1,4)Γ(4,2; 3)𝑃(1,2) = −𝑖∫ 𝑑34 𝐺(1,3)𝐺(4,1+ )Γ(3,4; 2)𝑊(1,2) = 𝑉𝐶 (1,2) + ∫ 𝑑34𝑉𝐶 (1,3)𝑃(3,4)𝑊(4,2)Γ(1,2; 3) = 𝛿(1,2)𝛿(1,3) + ∫ 𝑑4567

𝛿Σ xc (1,2)
𝐺(4,6)𝐺(7,5)Γ(6,7; 3)
𝛿𝐺(4,5)

By self-consistently solving them, one can rigorously obtain the Green’s function. When thevertex function is approximated as the product of delta functions, it corresponds toapproximating as Σ xc (1,2) = 𝑖∫ 𝑑34𝑊(1+ , 3)𝐺(1,4). This is the reason why it is called the GW
approximation. This approximation also means to taking 𝑃(1,2) = −𝑖∫ 𝑑34 𝐺(1,3)𝐺(4,1+ ),which indicates two particles are propagating independently. As a higher order approximation,
one can use the first order approximation
Γ(1,2; 3) = 𝛿(1,2)𝛿(1,3) + ∫ 𝑑4567

𝛿Σ xc (1,2)
𝐺(4,6)𝐺(7,5)
𝛿𝐺(4,5)to obtain the polarization. The resulting set of equations is called the GW + Bethe Salpeter
equation.
To describe the equation, it is convenient to use the four-point response function definedabove𝜒(1,2; 1′ , 2′) = 𝐿(1,1′ , 2,2′ ) = i𝐺2 (1,1′ , 2,2′ ) − 𝐺(1,1′ )𝐺 ′ 2′ , 2).
It can be derived by Schwinger that this quantity is given as a response of the Green’s function36

𝐺(1,1′ ) perturbed by a nonlocal external potential 𝑢(2,2) as
𝐿(1,1′ , 2,2′ ) = −iδG(1,1′ ).
δu(2,2)Using the identity of the functional derivative,𝛿𝐹(1,1′ )𝛿𝐹 −1 (3,3′ )= − ∫ 𝐹(1,3)
𝐹(3′ , 1′ )𝑑3𝑑3′ ,
𝛿𝐺(3)
𝛿𝐺(2)
we obtain
δG(1,1′ )δG−1 (3,3′ )= − ∫ 𝐺(1,3)𝐺(3′ , 1′ )𝑑3𝑑3′ .
δu(2,2′)δu(2,2′)
When this is substituted into the Dyson equation
𝐺 −1(3,3′ ) = 𝐺0−1 (3,3′ ) − 𝑢(3,3′ ) − Σ(3,3′ ),We obtain
δG(1,1′)δΣ(3,3′ )= G(1,2)G(2′ , 1) + ∫ 𝐺(1,3)𝐺(3′ , 1′ )𝑑3𝑑3′
δu(2,2′)δu(2,2′)= G(1,2)G(2′ , 1) + ∫ 𝐺(1,3)𝐺(3′ , 1′ )δΣ(3,3′ ) δG(4,4′)𝑑3𝑑3′ 𝑑4𝑑4′.
δG(4,4′) δu(2,2′)

The first term in the left-hand side, or the product of the Green’s functions, is regarded as thefour-point response function within the independent pair approximation. Therefore, it will bedenoted as 𝐿0 (1,1′ , 2,2′ ). The second term contains the derivative of the self-energy, which
will be denoted as
𝐾(3,3′ , 4,4′) = i

δΣ(3,3′ ).
δG(4,4′)This quantity is the electron-hole interaction kernel because it relates the non-interacting
response function with the interacting one as
𝐿(1,1′ , 2,2′ ) = 𝐿0 (1,1′ , 2,2′) + ∫ 𝐿0 (1,1′ , 3,3′ )𝐾(3,3′ , 4,4′ )𝐿(4,4′ , 2,2′ )𝑑3𝑑3′ 𝑑4𝑑4′.
The kernel 𝐾 can be decomposed into the Hartree and exchange-correlation contributions.
The former
𝐾 x (3,3′ , 4′ ) ≡ 𝑖

δΣ H (3,3′ )δG(4,4′ )is called bare Coulomb exchange interaction and can be described as
𝐾 x (3,3′ , 4′ ) = 𝛿(3,3′ )𝛿(4,4′ )𝑣(3,4).
While the latter is given within the 𝐺𝑊 approximation as𝑖δΣ 𝐺𝑊 (3,3′)δW(3+ , 3′ )= −𝛿(3,4)𝛿(3′ , 4′ )𝑊(3+ , 3′ ) − 𝐺(3,3′ ).
′δG(4,4 )δG(4,4′ )The functional derivative, the second term in the right-hand side, can be removed by using37

the polarization, 𝑊 = 𝑣 + 𝑣𝑃𝑊, which provides the relation 𝛿𝑊 ⁄𝛿𝐺 = 𝑊 𝛿𝑃⁄𝛿𝐺 𝑊 as𝑖δΣ 𝐺𝑊 (3,3′ )δ𝑃(5,6)= −𝛿(3,4)𝛿(3′ , 4′ )𝑊(3+ , 3′ ) − ∫ 𝐺(3,3′ )𝑊(3+ , 5)𝑊(6,3′ )𝑑5𝑑6.
′
δG(4,4 )δ𝐺(4,4′ )Using the approximation 𝑃 = −𝐺𝐺,𝑖δΣ 𝐺𝑊 (3,3′ )= −𝛿(3,4)𝛿(3′ , 4′ )𝑊(3+ , 3′ )δG(4,4′ )+ i ∫ 𝐺(3,3′ )𝑊(3+ , 5)

δ(𝐺(5,6)𝐺(6,5+ ))𝑊(6,3′)𝑑5𝑑6,
δ𝐺(4,4′ )and further,𝑖δΣ 𝐺𝑊 (3,3′ )= −𝛿(3,4)𝛿(3′, 4′ )𝑊(3+ , 3′ ) + 𝑖𝑊(3+ , 4)𝐺(4′ , 4+ )𝑊(4′ , 3′ )δG(4,4′ )+ 𝑖𝑊(3+ , 4′ )𝐺(4′ , 4+ )𝑊(4, 3′ ).

The first term is called the direct term, and the second and third terms are called secondorder correction terms.
𝐾 x (3,3′ , 4,4′ ) = 𝛿(3,3′ )𝛿(4,4′)𝑣(3,4)
𝐾 d (3,3′ , 4,4′ ) = −𝛿(3,4)𝛿(3′ 4′ )𝑊(3+ , 3′ )′′𝐾 ′(3,3 ,4,4 ) = −𝐺(3,3′ )𝛿𝑊(3,3+ ).
𝛿𝐺(4,4′ )It is important that this GW + Bethe Salpeter equation can be rewritten using the Kohn-Shamorbitals as basis. The resulting equation is formally the same as the response function
appeared above in TD-DFT,
𝜒(𝑟𝑡, 𝑟 ′ 𝑡 ′ ) = 𝜒KS (𝑟𝑡, 𝑟 ′ 𝑡 ′ ) + ∭ 𝜒𝐾𝑆 (𝑟𝑡, 𝑟1 𝑡1 ) [𝛿(𝑡1 − 𝑡2 )+ 𝐾xc (𝑟1 𝑡1 , 𝑟2 𝑡2 )] 𝜒(𝑟2 𝑡2 , 𝑟 ′ 𝑡 ′ ).
|𝑟1 − 𝑟2 |When we take the standpoint that the response function is given as a functional of the electron
density through the KS orbitals, this is a density functional theory with advanced form for theexchange-correlation.
Direct random phase approximation
Let us go back to the density functional theory again. One of the crude approximations to
the response function was to take the direct term only, which is within the adiabatic
connection path
𝜆 (𝒓
𝜆
𝐾Hxc
1 , 𝒓2 ; 𝜔) ≃ 𝐾H (𝒓1 , 𝒓2 ) = 𝜆𝑣(𝒓1 , 𝒓2 ).

(4.28)

When this is used, the resulting equation for the response function is𝜒𝜆 (𝒓1, 𝒓2 ; 𝜔) = 𝜒0 (𝒓1 , 𝒓2 ; 𝜔) + 𝜆 ∭ 𝜒0 (𝒓1 , 𝒓3 ; 𝜔)𝑣(𝒓3 , 𝒓4 )𝜒𝜆 (𝒓4 , 𝒓2 ; 𝜔).

(4.29)

By this we can define the correlation functional within this approximation called direct RPA
38

as
∞
1 1
𝑑𝜔 −i𝜔0+𝜒0 (𝒓1 , 𝒓3 ; 𝜔)𝜒0 (𝒓4 , 𝒓2 ; 𝜔)
𝐸cdRPA = − ∫ 𝑑𝜆 ∫𝑒
[𝜆 ∭𝑑𝒓1 𝑑𝒓2 𝑑𝒓3 𝑑𝒓4|𝒓1 − 𝒓2 ||𝒓3 − 𝒓4 |
2 0
∞ 2𝜋i

+ 𝜆2 ∭

𝜒0 (𝒓1 , 𝒓3 ; 𝜔)𝜒0 (𝒓4 , 𝒓5 ; 𝜔)𝜒0 (𝒓6 , 𝒓2 ; 𝜔)𝑑𝒓1 𝑑𝒓2 𝑑𝒓3 𝑑𝒓4 𝑑𝒓5 𝑑𝒓6|𝒓1 − 𝒓2 ||𝒓3 − 𝒓4 ||𝒓5 − 𝒓6 |+ ⋯ ] (4.30)
Practical formulation of dRPA
Equation (4.29) can be expressed using the four-point response function asi𝜒𝜆 (𝒙𝟏 , 𝒙2 , 𝒙′𝟏 , 𝒙′2 ; 𝜏)= ⟨Ψ 𝜆 |𝑇[𝑛̂1𝜆 (𝒙1 , 𝒙1′ ; 𝑡1 )𝑛̂1𝜆 (𝒙2 , 𝒙′2 ; 𝑡2 )]|Ψ 𝜆 ⟩− ⟨Ψ 𝜆 |𝑛̂1𝜆 (𝒙1 , 𝒙1′ ; 𝑡1 )|Ψ 𝜆 ⟩⟨Ψ 𝜆 |𝑛̂1𝜆 (𝒙2 , 𝒙′2 ; 𝑡2 )|Ψ 𝜆 ⟩.

(4.34)

Within dRPA the Dyson equation for the response function is𝜒𝜆dRPA (𝒙1 , 𝒙2 , 𝒙1′ , 𝒙′2 ; 𝜔)
= 𝜒0 (𝒙1 , 𝒙2 , 𝒙1′ , 𝒙′2 ; 𝜔)+ ∭ 𝑑𝒙3 𝑑𝒙4 𝑑𝒙5 𝑑𝒙6 𝜒0 (𝒙1 , 𝒙4 , 𝒙1′ , 𝒙3 ; 𝜔)𝐾H𝜆 (𝒙3 , 𝒙6 ; 𝒙4 , 𝒙5 ) 𝜒𝜆dRPA (𝒙5 , 𝒙2 , 𝒙6 , 𝒙′2 ; 𝜔) (4.36)
with the kernel being given by
𝐾H𝜆 (𝒙1 , 𝒙2 , 𝒙1′ , 𝒙′2 ) = 𝜆𝑣(|𝒓1 − 𝒓2 |)𝛿(𝒙1 − 𝒙′1 )𝛿(𝒙2 − 𝒙′ 2 )Note that the non-interacting kernel can be given by the KS orbitals as𝜒0 (𝒙1 , 𝒙2 , 𝒙1′ , 𝒙′2 ; 𝜔) = ∑ [𝑖≤𝑁<𝑎𝜙𝑖∗ (𝒙1′ )𝜙𝑎 (𝒙1 )𝜙𝑎∗ (𝒙′2 )𝜙𝑖 (𝒙2) 𝜙𝑖∗ (𝒙′2 )𝜙𝑎 (𝒙2 )𝜙𝑎∗ (𝒙1′ )𝜙𝑖 (𝒙1 )−]𝜔 − (𝜀𝑎 − 𝜀𝑖 ) + i0+𝜔 + (𝜀𝑎 − 𝜀𝑖 ) − i0+

Now we take a basis representation. The basis is formed using the combined occupied andunoccupied orbitals as
𝑓𝑖𝑎 (𝒙1 , 𝒙1′ ) = 𝜙𝑖∗ (𝒙1′ )𝜙𝑎 (𝒙1 )𝑓𝑎𝑖 (𝒙1 , 𝒙1′ ) = 𝜙𝑎∗ (𝒙1′)𝜙𝑖 (𝒙1 ).
Then 𝜒0 can be expanded as𝜒0 (𝒙1, 𝒙2 , 𝒙1′ , 𝒙′2 ; 𝜔) ∑[𝜒0 ]𝑝𝑞 𝑓𝑝 (𝒙1, 𝒙1′ )𝑓𝑞∗ (𝒙2, 𝒙′2 )
𝑝𝑞

where
[𝜒0 (𝜔)]𝑖𝑎,𝑗𝑏 =𝛿𝑖𝑗 𝛿𝑎𝑏𝜔 − (𝜀𝑎 − 𝜀𝑖 ) + i0+[𝜒0 (𝜔)]𝑎𝑖,𝑏𝑗 =𝛿𝑖𝑗 𝛿𝑎𝑏𝜔 + (𝜀𝑎 − 𝜀𝑖 ) − i0+[𝜒0 (𝜔)]𝑖𝑎,𝑏𝑗 = [𝜒0 (𝜔)]𝑎𝑖,𝑗𝑏 = 0.
39

Therefore, in a matrix representation(𝜒0 )−1(𝜔) = − [(𝚫𝜺 𝟎 ) − 𝜔 (𝟏 𝟎 )]
𝟎 𝚫𝜺
𝟎 −𝟏

(4.39)

with
Δ𝜀𝑖𝑎,𝑗𝑏 ≡ (𝜀𝑎 − 𝜀𝑖 )𝛿𝑖𝑗 𝛿𝑎𝑏 .
Then
𝐀λ
−1
(𝜒𝜆dRPA ) (𝜔) = − [( ∗
𝑩𝜆

𝑩𝜆
𝟏 𝟎
)−𝜔(
)]
𝐀∗𝜆
𝟎 −𝟏

(4.40)

with(𝐴𝜆 )𝑖𝑎,𝑗𝑏 = Δ𝜀𝑖𝑎,𝑗𝑏 + 𝜆⟨𝜙𝑎 𝜙𝑗 |𝜙𝑖 𝜙𝑏 ⟩(𝐵)𝑖𝑎,𝑗𝑏 = 𝜆⟨𝜙𝑎 𝜙𝑏 |𝜙𝑖 𝜙𝑗 ⟩.

(4.41)

To take the inverse in Eq. (4.40), we solve the generalized eigenvalue problem
𝐀λ
( ∗
𝑩𝜆

𝑿𝑛,𝜆
𝑩𝜆 𝑿𝑛,𝜆𝟏 𝟎)() = 𝜔𝑛𝜆 ())(
𝐀∗𝜆 𝒀𝑛,𝜆𝟎 −𝟏 𝒀𝑛,𝜆(4.42)

Under the normalization condition for the eigenvectors (X ,Y), the response function is(𝜔) = ∑ [
𝝌dRPA
𝜆
𝑛

1

𝑿𝑛,𝜆
(
) (𝑿†𝑛,𝜆𝜆+𝜔 − 𝜔𝑛 + 𝑖0 𝒀𝑛,𝜆𝒀†𝑛,𝜆 ) −

𝑿∗𝑛,𝜆
∗†
(
∗ ) (𝑿𝑛,𝜆𝜔 + 𝜔𝑛𝜆 − 𝑖0+ 𝒀𝑛,𝜆1𝒀∗†𝑛,𝜆 )],and the correlation part of the two-body density matrix is then∞𝒀∗𝑛,𝜆 𝒀∗†𝑑𝜔 i𝜔0+𝜆,dRPA[𝝌𝜆 (𝜔) − 𝝌0 (𝜔)] = ∑ [( ∗ 𝑛,𝜆
𝒏2,c
= −∫
𝑒
2𝜋i
𝑿𝑛,𝜆 𝒀∗†∞𝑛,𝜆𝑛

𝒀∗𝑛,𝜆 𝑿∗†
𝑛,𝜆

𝟎 𝟎
)]
𝟏

)−(
𝟎
𝑿∗𝑛,𝜆 𝑿∗†𝑛,𝜆Therefore, the correlation functional is given by
1 1
𝐸cdRPA = ∫ 𝑑𝜆
2 0

∑

∗

∗

∑⟨𝜙𝑖 𝜙𝑏 |𝜙𝑎 𝜙𝑗 ⟩(𝑌𝑛,𝜆 )𝑖𝑎 (𝑌𝑛,𝜆 )𝑗𝑏 + ⟨𝜙𝑖 𝜙𝑗 |𝜙𝑎 𝜙𝑏 ⟩(𝑋𝑛,𝜆 )𝑖𝑎 (𝑋)𝑗𝑏

𝑖𝑗≤𝑁<𝑎𝑏 𝑛∗+ ⟨𝜙𝑎 𝜙𝑏 |𝜙𝑖 𝜙𝑗 ⟩(𝑋𝑛,𝜆 )𝑖𝑎 (𝑌𝑛,𝜆 )𝑗𝑏 + ⟨𝜙𝑖 𝜙𝑏 |𝜙𝑎 𝜙𝑗 ⟩ [(𝑋)∗𝑖𝑎 (𝑋𝑛,𝜆 )𝑗𝑏 − 𝛿𝑖𝑗 𝛿𝑎𝑏 ].
This is the formulation within the direct RPA. There are other formulations as well depending
on how the 𝜆 integration is done efficiently by introducing approximations/assumptions.
Most of the dRPA calculations are done non-self-consistently using a lower-level electron
density although self-consistent calculation can be done using the optimized effective
potential scheme. With this, we can incorporate the van der Waals interaction. Yet, thisscheme is not free-from the self-interaction, that will affect insulators more significantly. To
overcome it, one can add the exchange to the direct term as
λ (𝒓
λ
𝐾Hxc
1 , 𝒓2 ; 𝜔) ≃ 𝐾Hx (𝒓1 , 𝒓2 ; 𝜔) = 𝜆𝑣(𝒓1 , 𝒓2 ) + 𝜆𝐾𝑥 (𝒓1 , 𝒓2 ; 𝜔).

This approach is found to correct the problem only minorly. It is thus necessary to go beyond
the RPA approach by adding terms that include the vertex correction.

40

Linear response within classical statistical physics
Here I discuss a theory of density correlation for a classical many-body stem to becompared with DFT. The method is based on the standard statistical physics of a grandcanonical ensemble. We start from the Ornstein Zernike (OZ) equation. For detail, please
read the article written by Prof. Sato, (http://j-molsci.jp/archives/AC0004.pdf).
We now see how the correlation factors are described in a classical system(a) fluctuation and correlation
Let us define the density operator𝑁𝛿𝑛(𝑟⃗) = ∑ 𝛿(𝑟⃗ − 𝑟⃗𝑖 ) − 𝑛0 ,𝑖=1where 𝑛0 is the density averaged over the whole system. The average of 𝛿𝑛(𝑟⃗) over theensemble is zero, or𝑁〈∑ 𝛿(𝑟⃗ − 𝑟⃗𝑖 )〉 = 𝑛0 .
𝑖=1

The density fluctuation can be written as
𝑁𝑁𝑁〈𝛿𝑛(𝑟⃗)𝛿𝑛(𝑟⃗′)〉 = 〈∑ 𝛿(𝑟⃗ − 𝑟⃗𝑖 ) ∑ 𝛿(𝑟⃗′ − 𝑟⃗𝑗 )〉 − 2𝑛0 〈∑ 𝛿(𝑟⃗ − 𝑟⃗𝑖 )〉 + 𝑛02
𝑖=1

𝑗=1

𝑖=1

𝑁

𝑁

= 〈∑ 𝛿(𝑟⃗ − 𝑟⃗𝑖 ) 𝛿(𝑟⃗′ − 𝑟⃗𝑖 )〉 + 〈 ∑𝑖=1

𝛿(𝑟⃗ − 𝑟⃗𝑖 ) 𝛿(𝑟⃗′ − 𝑟⃗𝑗 )〉 − 𝑛02𝑖,𝑗=1,𝑖≠𝑗= 𝑛0 δ(𝑟⃗ − 𝑟⃗ ′ ) + 𝑛02 𝑔(𝑟⃗, 𝑟⃗ ′ ) − 𝑛02 = 𝑛0 δ(𝑟⃗ − 𝑟⃗ ′ ) + 𝑛02 ℎ(𝑟⃗, 𝑟⃗ ′ )where 𝑔 is the pair correlation function𝑁𝑔(𝑟⃗, 𝑟⃗′)𝛿(𝑟⃗ − 𝑟⃗𝑖 ) 𝛿(𝑟⃗′ − 𝑟⃗𝑗 )〉≡〈 ∑𝑖,𝑗=1,𝑖≠𝑗and ℎ is called as the total correlation function (two-body density matrix over square of theparticle density)ℎ(𝑟⃗, 𝑟⃗ ′ ) ≡ 𝑔(𝑟⃗, 𝑟⃗ ′ ) − 1.
Note that similar quantities have been discussed above although the naming is somewhat
different. (𝑛02 ℎ(𝑟⃗, 𝑟⃗ ′ ) ↔ 𝑛2 (𝒓, 𝒓′))To relate the correlation with the response function, it is convenient to introduce a functional
derivative scheme. In the first step, let us relate the correlation with the partition function.
(b) grand partition function and grand distribution function
41

Then, we focus on the grand partition function. For a system described by the Hamiltonian,
𝑁

𝐻=∑
𝑖

1 2
𝑝⃗ + 𝑈𝑁 (1, ⋯ , 𝑁),
2𝑚𝑖 𝑖

the partition function is written as, using the thermal de Broglie wave length 𝛬 =ℎ/√2𝜋𝑚𝑘𝐵 𝑇,
∞

Ξ= ∑∭
𝑁=0

1
Λ3𝑁 𝑁!

∞

= ∑∭
𝑁=0

exp[−𝛽(𝑈𝑁 (1, ⋯ , 𝑁) − 𝜇)] 𝑑(1 … 𝑁)𝑧𝑁
exp[−𝛽𝑈𝑁 (1, ⋯ , 𝑁)] 𝑑(1 … 𝑁),𝑁!

where 𝑧 is the activity𝑧 = Λ−3 exp(𝛽𝜇),and the 𝑚-particle distribution function (classical analogue of the 𝑚-body density matrix) is
given as
∞

𝑛

(𝑚)

1
𝑧𝑁
(1, ⋯ , 𝑚) = ∑ ∭ 𝑑(1 … 𝑁) exp[−𝛽𝑈𝑁 (1, ⋯ , 𝑁)]
Ξ
𝑁!
𝑁=0

× ∑ ∑ ⋯ ∑ 𝛿(1 − 𝑟⃗𝑖 )𝛿(2 − 𝑟⃗𝑗 ) ⋯ 𝛿(𝑛 − 𝑟⃗𝑘 ).
𝑖𝑗≠𝑖

𝑘≠𝑖,𝑗,⋯

This can be rewritten as a sum of multiples of the canonical partition function 𝑍𝑁 and the(𝑚)canonical distribution function 𝑛𝑁

∞

𝑛

(𝑚)

1
𝑧𝑁
(𝑚)
(1, ⋯ , 𝑚) = ∑𝑍 𝑛 (1, ⋯ , 𝑚)Ξ
𝑁! 𝑁 𝑁
𝑁=0

Note that(𝑚)∭ 𝑑(1 … 𝑁)𝑛𝑁 (1, ⋯ , 𝑚) =𝑁!
.
(𝑁 − 𝑚)!

Now that the correlation functions are given in terms of the partition function, let us consider
the functional derivative method.
(c) Functional derivative
When the potential energy consists of the external potential and the rest as𝑁

𝑈𝑁 (1, ⋯ , 𝑁) = ∑ 𝑉(𝑖) + 𝑈′𝑁 (1, ⋯ , 𝑁)
𝑖=1

the grand partition function can be written using generalized activity𝑧(𝑟⃗) = 𝑧 exp[−𝛽𝑉(𝑟⃗)]𝑧(𝑖) = 𝑧 exp[−𝛽𝑉(𝑟⃗𝑖 )]
as

42

∞

𝑁

𝑁=0

𝑖=1

1
Ξ = ∑ ∭ (∏ 𝑧(𝑖)) exp[−𝛽(𝑈 ′ 𝑁 (1, ⋯ , 𝑁))] 𝑑(1 … 𝑁).
𝑁!
When this is differentiated by 𝑧(𝑟⃗), we get
∞

𝑁

δΞ
1
δ𝑧(𝑖)
= ∑∭ ∑( ∏ 𝑧(𝑗)) exp[−𝛽(𝑈 ′ 𝑁 (1, ⋯ , 𝑁))] 𝑑(1 … 𝑁)δ𝑧(𝑟⃗)𝑁!
δ𝑧(𝑟⃗)
𝑁=0𝑖

𝑗=1,𝑗≠𝑖∞𝑁1= ∑ ∭ ∑ 𝛿(𝑟⃗𝑖 − 𝑟⃗) ( ∏ 𝑧(𝑗)) exp [−𝛽(𝑈′𝑁 (1, ⋯ , 𝑁))]𝑑(1 … 𝑁)
𝑁!
𝑁=0

𝑖

𝑗=1,𝑗≠𝑖

∞

𝑁

𝑁=0

𝑗=2

1
= ∑∭
(∏ 𝑧(𝑗)) exp[−𝛽(𝑈 ′ 𝑁 (𝑟⃗, 2, ⋯ , 𝑁))] 𝑑(2 … 𝑁).
(𝑁 − 1)!
Comparing with the above equations for the distribution functions, we can show𝑛(1) (𝑟⃗) =𝑛(2)𝑧(𝑟⃗) δΞδ ln Ξ
=Ξ δ𝑧(𝑟⃗) δ ln 𝑧(𝑟⃗)(𝑟⃗, 𝑟⃗ ′ ) = 𝑧(𝑟⃗′)δ𝑛(1) (𝑟⃗) δ𝑛(1) (𝑟⃗)=.
δ𝑧(𝑟⃗′)
δ ln 𝑧(𝑟⃗′)Note that δ ln 𝑧(𝑟⃗) = −𝛽δ𝑉(𝑟⃗), so that𝑛(2) (𝑟⃗, 𝑟⃗ ′ ) = −𝑘𝐵 𝑇δ𝑛(1) (𝑟⃗).
δ𝑉(𝑟⃗′)

This is a fluctuation dissipation theorem in the sense that the two-body density is related to
the response of the density to the local potential.
Let us then consider inverse of the response function
δ𝑛(1) (𝑟⃗)−𝛽.
δ𝑉(𝑟⃗′)
This should be the inverse of the non-interacting response function minus the interaction
kernel. When no interaction exists,(1)

𝑛0 (𝑟⃗) = exp[−𝛽𝑉(𝑟⃗)],so that𝜒 (0)−1 (𝑟⃗, 𝑟⃗ ′ ) = −𝛽𝛿𝑉(𝑟⃗)
𝛿(𝑟⃗ − 𝑟⃗ ′ )=.
(1)
𝛿𝑛(1) (𝑟⃗ ′ )𝑛 (𝑟⃗)
0Deviation from this is the result of the interaction, which will be denoted by 𝐾(𝑟⃗, 𝑟⃗ ′ ) as
𝜒 −1 (𝑟⃗, 𝑟⃗ ′ ) ≡δ ln 𝑧(𝑟⃗)𝛿𝑉(𝑟⃗)
𝛿(𝑟⃗ − 𝑟⃗ ′ )=−𝛽≡− 𝐾(𝑟⃗, 𝑟⃗ ′ ).
δ𝑛(1) (𝑟⃗′)𝛿𝑛(1) (𝑟⃗ ′ )𝑛(1) (𝑟⃗)When we use the chain rule43∫𝛿𝑛(1) (𝑟⃗) 𝛿 ln 𝑧(𝑟⃗ ′′ )= 𝛿(𝑟⃗ − 𝑟⃗′)𝛿 ln 𝑧(𝑟⃗ ′′ ) 𝛿𝑛(1) (𝑟⃗′)we get∫(𝑛(𝑟⃗)δ(𝑟⃗ − 𝑟⃗ ′′ ) + 𝑛(𝑟⃗)ℎ(𝑟⃗, 𝑟⃗ ′′ )𝑛(𝑟⃗′′)) (𝛿(𝑟⃗′′ − 𝑟⃗ ′ )− 𝐾(𝑟⃗′′, 𝑟⃗ ′ )) 𝑑𝑟⃗′′ = 𝛿(𝑟⃗ − 𝑟⃗ ′ ),𝑛(𝑟⃗′′)

and further𝑛(𝑟⃗)ℎ(𝑟⃗, 𝑟⃗ ′ ) − 𝑛(𝑟⃗)𝐾(𝑟⃗, 𝑟⃗ ′ ) − 𝑛(𝑟⃗) ∫ ℎ(𝑟⃗, 𝑟⃗ ′′ )𝑛(𝑟⃗′′)𝐾(𝑟⃗ ′′ , 𝑟⃗ ′ ) 𝑑𝑟⃗ ′′ = 0ℎ(𝑟⃗, 𝑟⃗ ′ ) − 𝐾(𝑟⃗, 𝑟⃗ ′ ) − ∫ ℎ(𝑟⃗, 𝑟⃗ ′′ )𝑛(𝑟⃗′′)𝐾(𝑟⃗ ′′ , 𝑟⃗ ′ ) 𝑑𝑟⃗ ′′ = 0
which is the OZ equation to relate the inverse of the response function with the density
correlation via a Dyson-like equation. Note that this derivation shows that OZ does not
require a pair-wise interaction. When the interaction kernel 𝐾 is given as a functional of thedensity 𝑛(𝒓) using the DFT for classical particles, the correlation ℎ is determined. In theclassical DFT, there are known approximate and explicit relation of 𝐾 and ℎ when theparticles are interacting via a two-body force. The relation is called as a closure and the famous
ones are the hypernetted chain closure and the Percus-Yevic closure.
Nonadiabatic coupling of electrons and nuclei (ions)
Here we consider the whole system, consisting of electrons and ions, quantum mechanicallyalthough ions have been treated classically so far. We will go beyond the Born-Oppenheimer(BO) approximation here. By the way, what was the BO approximation?
BO approximation is explained as originating from large difference in mass between an
electron and a nucleus; even for proton the ratio is 1/1700. For heavier elements, the ratio
amounts to that of Earth over Sun. It is no wonder to decouple the motion of nuclei from thatof electrons. Even though BO approximation will be generally valid to first orderapproximation, there are phenomena related to breaking down of BO approximation. Suchoccurs for chemical reactions, electronic friction, and so on.
Now let us discuss it quantitatively using the BO theory.
As the starting point, let us consider the system witĥ=−𝐻

ℏ2 2 ℏ2 2∇ −∇ + 𝑉(𝑅, 𝑟)
2𝑀 𝑅 2𝑚 𝑟
44

where an electron with coordinate 𝑟 and mass 𝑚 and a nucleus with coordinate 𝑅 and mass
𝑀 is interacting via 𝑉(𝑅, 𝑟). The time-dependent Schrödinger equation iŝ Ψ(𝑟, 𝑅, 𝑡).
𝑖ℏ𝜕𝑡 Ψ(𝑟, 𝑅, 𝑡) = 𝐻
When 𝑀 ≫ 𝑚, naïve approximation is to neglect the first term in the Hamiltonian and treat
the nuclear coordinate classically. It can be also assumed that the nuclear velocity is muchsmaller than the electronic one, by which the electrons will have enough time to deexcite to
the ground state via coupling to environmental degrees of freedom, like electromagnetic field,
phonon, and so on. In that case, one need to solve a static equation
̂𝑒𝑙 (𝑅(𝑡)) = −𝐻ℏ2 2∇ + 𝑉(𝑟; 𝑅(𝑡))2𝑚 𝑟̂el Ψel (𝑟; 𝑅(𝑡)) = 𝐸𝑒𝑙 (𝑅(𝑡))Ψel (𝑟; 𝑅(𝑡))𝐻where 𝐸𝑒𝑙 (𝑅(𝑡)) is the ground state energy of electrons for a given configuration of nuclei,
which is usually called the BO potential energy surface. The classical nuclei will evolve withtime according to the Newton equation𝑀𝜕𝑡2 𝑅(𝑡) = −∇R 𝐸𝑒𝑙 (𝑅(𝑡)).
Let us derive this more formally. To do this, we prepare a complete set of the electronic
Hamiltonian and expand the whole wave function asΨ(𝑟, 𝑅, 𝑡) = ∑ 𝜒𝑛 (𝑅, 𝑡)𝜓𝑛 (𝑟; 𝑅),𝑛wherêel 𝜓n (𝑟; 𝑅) = 𝐸𝑛 (𝑅)𝜓𝑛 (𝑟; 𝑅).
𝐻
Note that the nuclear coordinate is treated quantum mechanically through the nuclear wave
function 𝜒𝑛 (𝑅, 𝑡) , which is introduced as the expansion coefficient. 𝜓𝑛 (𝑟; 𝑅) is timeindependent and is given as a function of parameter 𝑅; time-dependence is considered only
through 𝜒.
The nuclear wave function follows an equation of motion, which can be straightforwardly
obtained by applying the whole wave function into the starting time-dependent Schrödinger
equation𝑖ℏ𝜕𝑡 𝜒𝑚 (𝑅, 𝑡) = [−−

ℏ2 2
∇ + 𝐸𝑚 (𝑅)] 𝜒𝑚 (𝑅, 𝑡)
2𝑀 𝑅

ℏ2
ℏ2
∑⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩∇𝑅 𝜒𝑛 (𝑅, 𝑡) −∑⟨𝜓𝑚 |∇2𝑅 |𝜓𝑛 ⟩𝜒𝑛 (𝑅, 𝑡)
𝑀
2𝑀
𝑛

𝑛

45

(1).

When the first term in the right-hand side is larger than the second and third terms, thenuclear wave function of 𝑚-th excited state 𝜒𝑚 (𝑅, 𝑡) follows quantum equation of motion
given by the corresponding potential surface 𝐸m (𝑅). Then 𝜒𝑚 (𝑅, 𝑡) can be regarded as thewave packet. The centroid of 𝜒𝑚 (𝑅, 𝑡) follows the Newton equation. This corresponds to the
BO approximation.
Because of the second and third terms, the nuclear wave packets are interfering with eachother. The BO approximation thus breaks down. The second term represents a couplingthrough the first-order nonadiabatic coupling (NAC) vector⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩and the third term represents the coupling through the second-order NAC ⟨𝜓𝑚 |∇2𝑅 |𝜓𝑛 ⟩.
The first-order NAC vector can be rewritten as⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩ =̂ )|𝜓𝑛 ⟩⟨𝜓𝑚 |(∇𝑅 𝐻.
𝐸n − 𝐸𝑚This indicates that NAC should be very large when there is a degeneracy of the BO surface.
When the numerator is non-zero, the coupling diverges at the degenerate point, which isusually called as crossing point. The BO surfaces generally looks like a Dirac cone of graphene
because of linear behavior of the BO surfaces. Therefore, the crossing point is sometimes
called as a conical intersection.
The breakdown of the BO surfaces can happen at the crossing point even when the nuclear
mass is large. For metals, the crossing should occur. For example, when a charged particle
approaches a metal surface, the ground state will be mixed with excited states though whichenergy transfer occurs from the particle to the surface. The particle is then decelerated. Thisphenomenon is called electronic friction. This contrasts the friction generally occurring forinsulators, where the energy transfer occurs from the particle to phonon.
On the other hand, the coupling ⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩ is generally large for insulators although thecrossing does not usually occur between the ground state and the first-excited state; this isreflected from the large electron-phonon coupling. The crossing, however, can occur between
two excited states. Therefore, the mixing occurs in excited states dynamics. Suppose there isa potential surface that may lead to a dissociation of a molecule and another one that may not.
Depending on the velocity of the wave packet ∇𝑅 𝜒𝑛 (𝑅, 𝑡), the second term in the Hamiltonian(Eq. (1)) can be large or small; the faster the velocity the larger the coupling is. This explains
behavior of fast chemical reaction occurring after photo-excitation.

46

How to access the NAC vector
To get the value of NAC vector, one needs to evaluate the change of the many-body wave
function (eigenstate of the electronic Hamiltonian) with respect to displacement of nucleus,⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩, which is usually too difficult to achieve. In this context, let us recall the equation
−𝑍𝑒 2
⟨𝜓
(
)
⟩ ⟨𝜓 |ℎ(𝑅)|𝜓 ⟩
|∇
̂
𝑚
𝑅
|(∇
𝐻
)|𝜓
⟨𝜓𝑚 𝑅
|𝑟 − 𝑅| |𝜓𝑛𝑛⟩𝑚𝑛⟨𝜓𝑚 |∇𝑅 |𝜓𝑛 ⟩ ==≡,𝐸𝑛 − 𝐸𝑚𝐸𝑛 − 𝐸𝑚𝐸𝑛 − 𝐸𝑚−𝑍𝑒 2

where ℎ(𝑟; 𝑅) ≡ ∇𝑅 (|𝑟−𝑅|). Let us then compare it with the density response induced by the
perturbation 𝑉ext (𝑟, 𝜔) = ℎ(𝑟; 𝑅)𝑒 𝑖𝜔𝑡 , that is,𝛿𝜌(𝑟, 𝜔) = ∫ 𝜒(𝑟, 𝑟 ′ , 𝜔)𝑉ext (𝑟 ′ , 𝜔)𝑑𝑟′
The response function is given by𝜒(𝑟, 𝑟 ′ , 𝜔) = ∑ [⟨𝜓0|𝜌̂(𝑟)|𝜓𝑛 ⟩⟨𝜓𝑛 |𝜌̂(𝑟 ′ )|𝜓0⟩ ⟨𝜓0|𝜌̂(𝑟 ′ )|𝜓𝑛 ⟩⟨𝜓𝑛 |𝜌̂(𝑟)|𝜓0⟩−],𝜔 − (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂
𝜔 + (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂𝑛𝑁
where 𝜌̂(𝑟) = ∑𝛼=1 𝛿(𝑟 − 𝑟𝛼 ). When taking ∫ 𝑑𝑟 𝛿𝜌(𝑟, 𝜔)𝑉ext (𝑟, 𝜔), the resulting quantity is

proportional to∑[𝑛⟨𝜓0 |𝑉ext |𝜓𝑛 ⟩⟨𝜓𝑛 |𝑉ext |𝜓0⟩ ⟨𝜓0|𝑉ext |𝜓𝑛 ⟩⟨𝜓𝑛 |𝑉ext |𝜓0 ⟩−]𝜔 − (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂
𝜔 + (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂
= ∑[𝑛⟨𝜓0|ℎ(𝑅)|𝜓𝑛 ⟩⟨𝜓𝑛 |ℎ(𝑅)|𝜓0⟩ ⟨𝜓0 |ℎ(𝑅)|𝜓𝑛 ⟩⟨𝜓𝑛 |ℎ(𝑅)|𝜓0 ⟩−]𝜔 − (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂
𝜔 + (𝐸𝑛 − 𝐸0 ) + 𝑖𝜂of which the numerator is square of the NAC vector. This indicates that the NAC vector canbe given by the susceptibility (correlation) of ℎ(𝑅).
According to the time-dependent DFT, this can be given via the Dyson equation
𝜒(𝑟, 𝑟 ′ , 𝜔) = 𝜒𝐾𝑆 (𝑟, 𝑟 ′ , 𝜔) + ∫ 𝑑𝑟1 𝑑𝑟2 𝜒𝐾𝑆 (𝑟, 𝑟1 , 𝜔)𝐾HXC (𝑟1 , 𝑟2 , 𝜔)𝜒(𝑟2 , 𝑟 ′ , 𝜔),
where
𝐾HXC (𝑟1 , 𝑟2 , 𝜔) =1
+ 𝑓xc (𝑟1 , 𝑟2 , 𝜔)|𝑟1 − 𝑟2 |and𝜒𝐾𝑆 (𝑟, 𝑟 ′ , 𝜔) = 2 ∑ [𝑖𝑎𝜑𝑖∗ (𝑟)𝜑𝑎 (𝑟)𝜑𝑎∗ (𝑟 ′ )𝜑𝑖 (𝑟 ′ ) 𝜑𝑎∗ (𝑟)𝜑𝑖 (𝑟)𝜑𝑖∗ (𝑟 ′ )𝜑𝑎 (𝑟 ′ )−],𝜔 − (𝜀𝑎 − 𝜀𝑖 ) + i𝜂
𝜔 + (𝜀𝑎 − 𝜀𝑖 ) − i𝜂

where 𝑖 (𝑎) runs over all the occupied (unoccupied) states. (Practically, one cannot take all
the unoccupied states into account. One needs to introduce a cutoff). Therefore, bycomputing the susceptibility using the response function of TD-DFT, one can obtain the47

value of NAC vector. There is a known compact form for the response function, known as
Casida equation, which can be used to conveniently. TD-DFT is accessible by Quantum
Espresso as well as may quantum chemical software packages like Gaussian.
Recall that the TD-DFT equation was derived from the response of the Kohn-Sham artificial
system to the change in the effective potential𝛿𝜌(𝑟, 𝑡) = ∫ 𝜒𝐾𝑆 (𝑟𝑡, 𝑟 ′ 𝑡′)𝛿𝑣eff (𝑟 ′ 𝑡′)𝑑𝑟′ = ∫ 𝜒𝐾𝑆 (𝑟𝑡, 𝑟 ′ 𝑡′)(𝛿𝑣ext (𝑟 ′ 𝑡′) + 𝛿𝑣HXC (𝑟 ′ 𝑡′))𝑑𝑟′= ∫ 𝜒𝐾𝑆 (𝑟𝑡, 𝑟 ′ 𝑡 ′ )(𝛿𝑣ext (𝑟 ′ 𝑡 ′ ) + 𝑓HXC (𝑟 ′ 𝑡 ′ , 𝑟 ′′ 𝑡 ′′ )𝛿𝑣ext (𝑟 ′ ′𝑡′′))𝑑𝑟′,
which is Fourier transformed to the previous equation. The quantity 𝑓HXC , or the Hartree
exchange-correlation kernel, is a functional of the time-dependent density. Note that the
adiabatic approximations (or ω-independent approximations) are usually adopted. Note alsothat the kernel formally plays a role of describing the interacting electrons and holes as does
the exchange-correlation self-energy in the theory of Hedin.
It is worth mentioning that the TD-DFT has been shown within the linear response scheme.
But the importance of TD-DFT is in the fact that non-linear response can be accessed despite
its limited accuracy originated from using approximate functionals. The linear responsemeans a linear deviation from the ground state, which is generally invalid for representing theexcited state density. The excited-state of a small molecules, for example, has a density that isvery different from that of the ground state. The excitation is no more a small perturbation to
the system. In that case, one should alternatively follow a real-time evolution of electrons,
which can be achieved by integrating the time-dependent Kohn-Sham equation with respect
to time.
The simplest expansion formula likeiℏ

φ𝑛 (𝑡 + Δ𝑡) − φ𝑛 (𝑡 − Δ𝑡)1= (− ∇2 + 𝑉eff (𝑟𝑡)) 𝜑𝑛 (𝑡)2Δ𝑡2

is not very appropriate in this case because this does not guarantee orthonormalization of theKohn-Sham orbitals. In this context, it is convenient to use
φ𝑛 (𝑡 + Δ𝑡) = exp [−1
1(− 𝛻 2 + 𝑉eff (𝑟𝑡)) Δ𝑡] 𝜑𝑛 (𝑡)iℏ2≃ exp [−1
1
1
1
1(− 𝛻 2 ) Δ𝑡] exp [− 𝑉eff (𝑟𝑡)Δ𝑡] exp [−(− 𝛻 2 ) Δ𝑡] 𝜑𝑛 (𝑡)
2iℏ
2
iℏ
2iℏ
2

which guarantees the orthonormalization and thus advantageous for parallel computation.
48

When applying the classical approximation to nuclei aŝel Ψel (𝑟; 𝑅(𝑡))iℏ ∂t Ψel (𝑟; 𝑅(𝑡)) = 𝐻
𝑀𝜕𝑡2 𝑅(𝑡) = −∇R 𝐸𝑒𝑙 (𝑅(𝑡))this coupled electronic and nuclear equation is called Ehrenfest equation of motion. This is akind of mean-field approximation to the coupling of electrons and ions. In integrating thiscoupled electron-nucleus equation, one will encounter a problem of very different time-scale
for the electronic and nuclear motion. The large difference in the time-scale makes thenumerical approach very unstable.
It is worth mentioning that the mixed classical-quantum approach like the Ehrenfest is known
to introduce a fundamental problem; as note in many textbook of quantum mechanics, nuclear
force derived thereby can be meaningless. Nuclear motion will be driven by the “force fromelectron” averaged over quantum states of electrons, which needs to be handled carefully.
When, for example, there are two electronic states for which a nucleus will feel a dissociative
force and a binding force, respectively. The nucleus will feel just an averaged force although,in proper quantum mechanical approach, the nuclear wave functions of dissociative andbound character should be used to describe the nuclear state. For a system of very small
degrees of freedom, there are ways to avoid such problem but is not practical for large systems.
There are movies on the Ehrenfest TD-DFT simulations, for example, on the internet.
https://figshare.com/articles/Physical_Factors_Affecting_Charge_Transfer_at_the_Pe_CO
OH_TiO_sub_2_sub_Anatase_Interface/2238757
One can observe a charge-transfer phenomenon that occurs between the adsorbate molecule
Pe-COOH to a titania surface. The authors investigated how the charge transfer dynamics is
affected by nuclear motion.
Quantum mechanical behavior of nuclei
Although nuclei are heavy, they behave as a classical particle at high temperatures. However,this is not always the case. Especially for hydrogen, or proton, the quantum effect is negligibly
small. To take the quantum effect into account, at least approximately, one can use the path
integral method.
For this purpose, let us start from the Hamiltonian𝑃̂2
̂ = ∑ 𝐼 + 𝐻𝑒 (𝑟̂ , 𝑝̂ ; 𝑅̂ )
𝐻
2𝑀
𝐼

49

and adiabatic states |𝜓𝑖 (𝑅)⟩. By decoupling the whole system into electronic and nuclear
states as|𝜓𝑖 (𝑅)⟩|𝑅⟩ ≡ |𝑖, 𝑅⟩
the following equation is assumed to hold𝐻𝑒 (𝑟̂ , 𝑝̂ ; 𝑅̂ )|𝜓𝑖 (𝑅)⟩|𝑅⟩ = 𝐸𝑖 (𝑅)|𝜓𝑖 (𝑅)⟩|𝑅⟩.
We are interested in an equilibrium density at temperature 𝑇. Then we will derive a formulafor the partition function̂ )] = ∫ 𝑑𝑅 ∑⟨𝑖, 𝑅| exp(−𝛽𝐻̂ ) |𝑖, 𝑅⟩
𝑍 = Tr[exp(−𝛽𝐻𝑖When the exponential is partitioned into 𝑃 parts by inserting the identity∫ 𝑑𝑅𝑖 ∑|𝑖𝑅⟩⟨𝑖𝑅|,𝑖

we get𝛽̂𝑍 = ∫ 𝑑𝑅1 ⋯ 𝑑𝑅𝑃 ∑ ⟨𝑖, 𝑅| exp (− 𝐻
) |𝑖𝑃 , 𝑅𝑃 ⟩𝑃
𝑖1 ,⋯,𝑖𝑃𝛽̂𝛽̂
× ⟨𝑖𝑃 , 𝑅𝑃 | exp (− 𝑃 𝐻) |𝑖𝑃−1 , 𝑅𝑃−1 ⟩ ⋯ ⟨𝑖1 , 𝑅1 | exp (− 𝑃 𝐻) |𝑖, 𝑅⟩
Then we need to evaluate terms like𝛽̂ρ(s, s − 1; τ) ≡ ⟨𝑖𝑠 , 𝑅𝑠 | exp (− 𝐻
) |𝑖𝑠−1 , 𝑅𝑠−1⟩
𝑃
This can be approximated as𝑃̂2𝛽ρ(s, s − 1; τ) ≃ ⟨𝜓0 (𝑅𝑠 )|𝜓0 (𝑅𝑠−1 )⟩ ⟨𝑅𝑠 | exp (− ∑𝐼 𝐼 ) |𝑅𝑠−1 ⟩𝑃 2𝑀𝐼× exp (−

𝛽 𝐸0 (𝑅𝑠 ) + 𝐸0 (𝑅𝑠−1 )) 𝛿𝑖𝑠 ,0 𝛿𝑖𝑠−1 ,0𝑃
2from the adiabatic sense. The overlap appearing above can be approximated by using Taylor
expansion as
𝑆𝑠,𝑠−1 ≡ ⟨𝜓0 (𝑅𝑠 )|𝜓0 (𝑅𝑠−1 )⟩ = ∫ 𝑑𝑟𝜓0∗ (𝑟; 𝑅𝑠 )𝜓0(𝑟; 𝑅𝑠−1 )1= 1 + ∑(𝑅𝐼,𝑠−1 , 𝑅𝐼,𝑠 )(𝑅𝐽,𝑠−1 , 𝑅𝐽,𝑠 ) ∫ 𝑑𝑟𝜓0∗ (𝑟; 𝑅𝑠 )∇𝐼,𝑠 ∇𝐼,𝑠−1 𝜓0 (𝑟; 𝑅𝑠−1 )2
𝐼,𝐽where we have assumed that the wave function is real, so that we can apply∫ 𝑑𝑟𝜓0∗ (𝑟; 𝑅)∇𝐼 𝜓0 (𝑟; 𝑅) = 0.
Since the derivative terms are small, we can reasonably assume that the overlap is small. Then
50

2

𝛽
𝑀𝐼 √𝑃
2
ρ(s, s − 1; τ) ≃ exp (− ∑ ( ) (𝑅𝐼,𝑠 − 𝑅𝐼,𝑠−1 ) )
𝑃
2 𝛽ℏ
𝐼

× exp (−

𝛽 𝐸0 (𝑅𝑠 ) + 𝐸0 (𝑅𝑠−1 )) 𝛿𝑖𝑠 ,0 𝛿𝑖𝑠−1 ,0𝑃
2Therefore, we have
3𝑃

2

𝑀𝐼 𝑃 2𝑀𝐼 √𝑃
1
2𝑍 = ∏ [() ] ∫ 𝑑𝑅1 ⋯ 𝑑𝑅𝑃 exp (−𝛽 ∑ ∑ ( ) (𝑅𝐼,𝑠 − 𝑅𝐼,𝑠−1 ) + 𝐸0 (𝑅𝑠 ))
2𝜋𝛽ℏ
2 𝛽ℏ
𝑃
𝐼

𝑠

𝐼

This has the same form as the partition function of classical systems having 𝑁 × 𝑃 particles
interacting via
2

𝑀𝐼 √𝑃
1
2
Veff = ∑ ∑ ( ) (𝑅𝐼,𝑠 − 𝑅𝐼,𝑠−1 ) + 𝐸0 (𝑅𝑠 )
2 𝛽ℏ
𝑃
𝑠

𝐼

This means that, if we perform a molecular dynamics (MD) simulation using 𝑁𝑃 particles
under a constant temperature 𝑇, we get canonical ensemble of 𝑁 quantum particles. So, theproblem is then how to do the DFT MD simulation. In fact, the method is well established.
The quantum simulation can be done conveniently using, for example, PIMD developed by
Dr. Shiga.
http://ccse.jaea.go.jp/ja/download/pimd/index.jp.html
Maximally localized Wannier orbitals
In solid state physics, Wannier orbitals are constructed from the Bloch orbitals by taking a“Fourier transformation”. Wannier orbital is localized spatially with the decay rate dependingon the electronic structure; it decays exponentially for insulators and polynomially for metals.
The decay rate can be enhanced by adding different bands to the transformation. When thisis done maximally, the orbitals are localized in atomic scale and can be taken as atomic orbitals
under the crystalline environment. The resulting orbitals can be conveniently used as a basisof the strongly correlated electron models. In addition, it has a peculiar property surprisingly
suitable for the study of ferroelectrics.
Formulation
We begin by defining the Wannier orbitals taking the phase factor 𝜑 explicitly taken into
account.

|𝑅𝑛 ⟩ =

Ω∫ |𝜓𝑛𝑘 ⟩𝑒 𝑖𝜑(𝑘)−𝑖𝑘𝑅 𝑑𝑘
(2𝜋)3
𝐵𝑍

51

Here, 𝜓𝑛𝑘 (𝑟) is the Bloch orbital,
𝜓𝑛𝑘 (𝑟) = 𝑢𝑛𝑘 (𝑟)𝑒 𝑖𝑘𝑟 :𝑢𝑛𝑘 (𝑟) = 𝑢𝑛𝑘 (𝑟 + 𝑅),

𝑅 is the lattice vector and 𝜑(𝑘) is the phase factor that is periodic in BZ; 𝜑(𝑘) = 𝜑(𝑘 + 𝐺)
for reciprocal lattice vectors 𝐺. We can localize the orbitals more by including a unitarytransformation among bands:|𝑅𝑛 ⟩ =

Ω𝑘 |𝜓𝑖𝜑(𝑘)−𝑖𝑘𝑅𝑑𝑘∫ ∑ 𝑈𝑚𝑛
𝑚𝑘 ⟩𝑒
(2𝜋)3
𝐵𝑍

𝑚

The transformation matrix 𝑈 can be determined by minimizing a variance𝐿 = ∑(⟨0𝑛 |𝑟 2 |0𝑛 ⟩ − ⟨0𝑛 |𝑟|0𝑛 ⟩2 )𝑛

One can show that the matrix elements can be written using the periodic part of the Blochfunction and derivatives in the reciprocal space.

⟨𝑅𝑛 |𝑟|0𝑚 ⟩ = 𝑖Ω∫ ⟨𝑢𝑛𝑘 |∇𝑘 |𝑢𝑚𝑘 ⟩𝑒 𝑖𝑘𝑅 𝑑𝑘
(2𝜋)3
𝐵𝑍

⟨𝑅𝑛 |𝑟 2|0𝑚 ⟩ = 𝑖Ω∫ ⟨𝑢𝑛𝑘 |∇2𝑘 |𝑢𝑚𝑘 ⟩𝑒 𝑖𝑘𝑅 𝑑𝑘
(2𝜋)3
𝐵𝑍

The derivatives are in practice estimated using the finite mesh points spanning the Brillouinezone; this is the computational physics. In this case, we approximate the derivative by a finite
difference as𝛻𝑓(𝑘) = ∑ 𝑤𝑏 𝑏[𝑓(𝑘 + 𝑏) − 𝑓(𝑘)]𝑏where b’s are vectors chosen so that 𝑘 and 𝑘 + 𝑏 are neighboring mesh points. Similarly,𝛻 2 𝑓(𝑘) = ∑ 𝑤𝑏 𝑏[𝑓(𝑘 + 𝑏) − 𝑓(𝑘)]2𝑏This indicates that ⟨𝑢𝑛𝑘 |∇𝑘 |𝑢𝑚𝑘 ⟩ is approximated using overlaps of the form ⟨𝑢𝑛𝑘 |𝑢𝑚𝑘+𝑏 ⟩. It
was shown that⟨0𝑛 |𝑟|0𝑛 ⟩ = −1∑ 𝑤𝑏 𝑏 Im(ln⟨𝑢𝑛𝑘 |𝑢𝑛𝑘+𝑏 ⟩):𝑁

Wannier center𝑘,𝑏⟨0𝑛 |𝑟 2|0𝑛 ⟩ = −1∑ 𝑤𝑏 {1 − |⟨𝑢𝑛𝑘 |𝑢𝑛𝑘+𝑏 ⟩|2 + |Im(ln⟨𝑢𝑛𝑘 |𝑢𝑛𝑘+𝑏 ⟩)|2 }
𝑁
𝑘,𝑏

(𝑘,𝑏)

hold. Therefore, the overlaps, Mnm ≡ ⟨𝑢𝑛𝑘 |𝑢𝑚𝑘+𝑏 ⟩, plays the central role in calculating the
localized Wannier orbitals.

52Modern theory of polarization
1. Motivation: Intuitive definition of polarization will be
1
𝑷 = ∫ 𝑑𝑟𝒓𝜌(𝒓),𝛺 𝑐𝑒𝑙𝑙where 𝛺 is the volume of the cell and 𝜌(𝑟) is the charge density, but this is ambiguous forsolids. The result depends on how we take the origin and also how we take the unit cell. Thisis the well-known problem of the polarization. You may consider that the real materials are
finite, so that the polarization should be well-defined after integrating over the whole sample.
In that case, the integration can be transformed to the surface integration using the Green’s
theorem. Thus, the theory requires argument of the surface. But, there is a way to define thepolarization using the unit cell information only.
2. New approach: One can solve the problem of the ill-defined “crystalline polarization” byconsidering “flow of charge” associated with displacement of ions. Phenomenologically, thepolarization occurs associated with the ionic displacement, such that the polarization 𝑃 iszero at an initial configuration of ions but is nonzero at a displaced configuration. This fact
makes it possible to define the polarization by connecting the initial unpolarized state to thefinal polarized state continuously and integrate the associated flow of charge as
1
𝑃 = ∫ 𝑑𝑡 ∫ 𝑑𝑟𝒋(𝒓, 𝑡).
𝛺 𝑐𝑒𝑙𝑙This “adiabatic” definition has made it possible to modernize the theory of polarization.
2.1 (DFPT formulation) Before introducing the modern theory, let us apply the DFPT
explaining associated problems. To make the polarization to occur, we apply electric field withlarge wave length; because of the ill-defined nature as explained above, the wave length is let
to be finite and will be extrapolated to infinity after the calculation. Because of the electricfield, the atoms will be displaced in proportional to the field strength thereby generating apolarization 𝑃0 = 𝑍𝜏 𝑅. The nuclear charge 𝑍𝜏 can be defined alternatively as the derivative of
the polarization with respect to the atomic displacement 𝛿𝑅. This can be also described as
the second derivative of the total energy with respect to the electric field and the atomic
displacement
𝜕𝜌
𝜕𝑅

1

𝜕2𝜌

𝜕 2 𝐸tot
𝜕𝐸𝜕𝑅

. The electron density will also be changed by the electric filed by

𝜕𝜌
𝜕𝐸

𝛿𝐸 +

𝛿𝑅 + 2 𝜕𝐸𝜕𝑅 𝛿𝐸𝛿𝑅 + ⋯ : Because the first order changes are zero, the second order isimportant. Since the uniform electric field is also ill-defined, we practically apply the field offinite wave length and then extrapolated to infinity after calculation. From the second
derivative modifies the value of Z, as 𝑍𝜏 + 𝑍eff , which is called Born effective charge. This isa DFPT approach to the Born effective charge. The modern approach is different.
3. Formulation: Let us introduce a parameter λ connecting the unpolarized and the polarized53

states; you may regard λ as a parameter for the displacement. By applying a first-order
perturbation theory, change of the wave function can be written as⟨𝜓𝑚𝑘 |𝜕𝜆 |𝜓𝑛𝑙 ⟩|𝛿𝜓𝑛𝑘 ⟩ = −𝑖ℏ𝛿𝜆 ∑|𝜓𝑚𝑘 ⟩
𝐸𝑛𝑘 − 𝐸𝑚𝑘𝑚≠𝑛

Corresponding current is obtained by calculating the expectation value of the current operator
as,
𝑗𝑛 =

⟨𝜓𝑛𝑘 |𝑝|𝜓𝑚𝑙 ⟩⟨𝜓𝑚𝑘 |𝜕𝜆 |𝜓𝑛𝑙 ⟩
𝑑𝑃𝑛
𝑖ℏ𝑒𝜆̇
=
∑ ∫𝑑𝑘
+ 𝑐. 𝑐.
3(2𝜋)𝑑𝑡
𝐸𝑛𝑘 − 𝐸𝑚𝑘𝑚≠𝑛This can be rewritten using the periodic part of the Bloch wave function. The result is
𝑑𝑃𝑛
𝑖𝑒𝜆̇
=
∫𝑑𝑘⟨𝛻𝑘 𝑢𝑛𝑘 |𝜕𝜆 𝑢𝑛𝑘 ⟩ + 𝑐. 𝑐.
(2𝜋)3𝑑𝑡This is appealing in that the contribution from the unoccupied states disappears and that thederivative with respect to 𝑘 is conveniently accessed by the maximally localized Wannierorbital method.
When summing over all the occupied states 𝑛 and integration over 𝑡 (to connect theinitial and final states) we finally arrive at𝑒𝑃(𝜆) =Im ∑ ∫𝑑𝑘⟨𝑢𝑛𝑘 |𝛻𝑘 |𝑢𝑛𝑘 ⟩(2𝜋)3𝑛

𝛥𝑃 = 𝑃(1) − 𝑃(0) = 𝑃(1)This is equivalent to the quantity appearing in the Berry phase theory and thus we will regard𝑖⟨𝑢𝑛𝑘 |𝛻𝑘 |𝑢𝑛𝑘 ⟩ as a Berry connection or gauge potential. Note that this is the electronic
contribution to the polarization. We need to add the ionic contribution𝑒𝑃ion = ∑ 𝑍𝜏 𝑅𝜏 ,
Ω
𝜏

as
𝑃=

𝑒
𝑒
Im ∑ ∫𝑑𝑘⟨𝑢𝑛𝑘 |𝛻𝑘 |𝑢𝑛𝑘 ⟩ + ∑ 𝑍𝜏 𝑅𝜏
3
(2𝜋)
Ω
𝑛

𝜏

This way we arrive at this very useful formula.
4. Quantum of polarization: Suppose we are integrating over a path within a one-dimensional
Brillouine zone;
𝑃𝑛 =

π/L
𝑒
𝑒
𝜙𝑛 =
Im ∫ 𝑑𝑘⟨𝑢𝑛𝑘 |𝜕𝑘 |𝑢𝑛𝑘 ⟩
2𝜋
2𝜋
−π/L

In this case, the initial and final points are equivalent because of the periodicity. So, the above
integration path corresponds to a closed loop.
The phase of the wave function can be taken arbitrarily like𝑢̃𝑛𝑘 = 𝑒 −𝑖𝛽(𝑘) 𝑢𝑛𝑘
except that the difference at the boundary 𝛽(𝜋/𝐿) − 𝛽(−𝜋/𝐿) is a multiple of 2π. When thepolarization is calculated using this new wave function, we get
𝜙̃𝑛 = 𝜙𝑛 + 2𝜋𝑚:𝑚 ∈ Integer54indicating that the polarization can be defined within the modulo of e. In other words, wecannot definite uniquely the crystalline polarization, which may be understood as different
paths can be used to connect the equivalent points in the Brillouine zone.
When this discussion is extended to a 3D system, it can be shown that the polarization iswell defined within the modulo of𝑒𝑅/𝛺
where 𝑅 is a lattice vector
𝑅 = ∑ 𝑚𝑗 𝑅𝑗 .
𝑗

Please see a literature on modern theory of polarization for the proof.
The polarization can be rewritten using the relation 𝑟 = −π
L 𝑑𝑘

𝑃𝑛 = −𝑒 ∫π− L 2𝜋𝑖𝜕𝜕𝑘, as the real part ofπ/L𝑑𝑘⟨𝑢𝑛𝑘 |𝑟|𝑢𝑛𝑘 ⟩−π/L 2𝜋⟨𝑢𝑛𝑘 |−𝑖𝜕𝑘 |𝑢𝑛𝑘 ⟩ = −𝑒 ∫5. Wannier center: We have learned that the Wannier orbitals are
𝑑 3 𝑘 −𝑖𝑘𝑅
𝑑 3 𝑘 𝑖𝑘(𝑟−𝑅)
(𝑟)
𝑒
𝜓
=
𝛺
𝑒
𝑢𝑛𝑘 (𝑟).
∫𝑛𝑘
3
3
𝐵𝑍 (2𝜋)𝐵𝑍 (2𝜋)𝑤𝑛 (𝑟 − 𝑅) = 𝛺 ∫Then, it is possible to show that the polarization can be described by the position of the
Wannier center (after adding the contribution from the nuclear charge)𝑒𝑒𝑃𝑛 = (𝑟𝑛𝑅 − 𝑅) = (⟨𝑤𝑛𝑅 |𝑟|𝑤𝑛𝑅 ⟩ − 𝑅).
𝛺𝛺Note that the Wannier center represents the spatially distributed charge density as a point
charge (see figure). We can thus regard the polarization effectively as the charge assigned to
the Wannier center. Note that the calculation can be done efficiently when using the most
localized Wannier orbitals, instead of the standard one shown here, but that formulation isskipped here.

There are various ways to relate it to the flow of charge differently like55

This difference corresponds to the existence of modulus.
Density functional perturbation theory(http://iffwww.iff.kfa-juelich.de/~wortmann/IFF-Springschool45/C2.pdf)
Here we study the perturbation theory within DFT, namely DFPT. DFPT is different fromthe perturbation theory of noninteracting systems in that the density response should also betaken into account. I will introduce DFPT in relation to the polarization of materials.
1. Sternheimer equation for noninteracting systems(1)

The first-order change in the wave function 𝜓𝑖can be obtained by solving the Sternheimer

equation(1)(0)(𝐻 (0) − 𝜀 (0) ) |𝜓𝑖 ⟩ = −(𝐻(1) − 𝜀 (1) ) |𝜓𝑖 ⟩where 𝐻 0 and 𝜀 0 are unperturbed ones and 𝐻 (1) and 𝜀 (1) are the perturbed ones and
(0)

(0)

𝜀 (1) = ⟨𝜓𝑖 |𝐻 (1) |𝜓𝑖 ⟩.
This equation is solved by superimposing the unperturbed states as
(1)

(1)

(0)

(1)

(0)

(1)

(0)

|𝜓𝑖 ⟩ = ∑ 𝑐𝑖𝑗 |𝜓𝑗 ⟩ = ∑ 𝑐𝑖𝑗 |𝜓𝑗 ⟩ + ∑ 𝑐𝑖𝑗 |𝜓𝑗 ⟩𝑗where 𝐼 and 𝐼

⊥

𝑗∈𝐼 ⊥

𝑗∈𝐼

are states degenerate to(0)𝜓𝑖and other states, respectively. By thisdefinition,
(0)

(1)

(1)

(0)

(𝐻 (0) − 𝜀𝑖 ) |𝜓𝑖 ⟩ = ∑ 𝑐𝑖𝑗 (𝜀𝑗

(0)

(0)

− 𝜀𝑖 ) |𝜓𝑗 ⟩.

𝑗∈𝐼 ⊥Therefore,(0)(1)
𝑐𝑖𝑗 = −

(0)⟨𝜓𝑗 |𝐻 (1) |𝜓𝑖 ⟩
(0)

𝜀𝑗

(0)

− 𝜀𝑖
(1)

for 𝑗 ∈ 𝐼 ⊥ .

Using the gauge degrees of freedom to set 𝑐𝑖𝑗 = 0 for 𝑗 ∈ 𝐼. Therefore,
56

(0)

(1)
(0)
|𝜓𝑖 ⟩ = − ∑ |𝜓𝑗 ⟩
𝑗∈𝐼 ⊥

(0)

⟨𝜓𝑗 |𝐻 (1) |𝜓𝑖 ⟩
(0)

𝜀𝑗

(0)

− 𝜀𝑖

.

This can be alternatively written as
(0)

(1)

(0)

(0)

(0)

𝑃𝐼⊥ (𝐻 (0) − 𝜀𝑖 ) 𝑃𝐼⊥ |𝜓𝑖 ⟩ = −𝑃𝐼⊥ 𝐻(1) |𝜓𝑖 ⟩ with 𝑃𝐼⊥ = ∑ |𝜓𝑗 ⟩ ⟨𝜓𝑗 |
𝑗∈𝐼 ⊥

or
−1

(1)

𝑃𝐼⊥ |𝜓𝑖 ⟩ = 𝐺𝐼⊥ (𝜀𝑖0 )𝐻(1) |𝜓𝑖0 ⟩ with 𝐺𝐼⊥ (𝜀) = [𝑃𝐼⊥ (𝜀 − 𝐻 (0) )𝑃𝐼⊥ ] .

2. DFT
In DFT, one has to take the density variation into account. The first-order change in theenergy is
𝑑
𝑑
(0)
(0)
∑⟨𝜓𝑖 |𝑇 + 𝑉ext + 𝑉H + 𝑉xc |𝜓𝑖 ⟩ = ⟨𝜓𝑖 |(𝑇 + 𝑉ext )(1)|𝜓𝑖 ⟩ + 𝐸Hxc [𝜌(0) ].
𝑑𝜆𝑑𝜆𝑖

The second order change is
(0)

(0)

(1)

(0)

(0)

(1)

⟨𝜓𝑖 |(𝑇 + 𝑉ext )(2) |𝜓𝑖 ⟩ + ⟨𝜓𝑖 |(𝑇 + 𝑉ext )(1) |𝜓𝑖 ⟩ + ⟨𝜓𝑖 |(𝑇 + 𝑉ext )(1) |𝜓𝑖 ⟩
(1)

(1)

+ ⟨𝜓𝑖 |(𝑇 + 𝑉ext )(0) |𝜓𝑖 ⟩1
𝛿 2 𝐸Hxc [𝜌(0) ] (1)𝑑 𝛿𝐸Hxc [𝜌(0) ] (1)
(1) ′
′
(𝑟)𝜌
(𝑟
)𝑑𝑟𝑑𝑟
+ ∬
𝜌
+
(
) 𝜌 (𝑟)𝑑𝑟∫2
𝛿𝜌(𝑟)𝛿𝜌(𝑟 ′ )
𝑑𝜆
𝛿𝜌(𝑟)
+

1 𝑑2
𝐸 [𝜌(0) ],2 𝑑𝜆2 Hxc

with
(1)∗

(𝑟)𝜓𝑖 (𝑟) + 𝜓𝑖

(1)

are obtained by solving the Sternheimer equation:

𝜌(1) (𝑟) = ∑ 𝜓𝑖

(0)

(0)∗

(1)

(𝑟)𝜓𝑖 (𝑟).

𝑖

The perturbed Kohn-Sham orbitals 𝜓𝑖
(0)

(1)

(1)

(0)

(𝐻 (0) − 𝜀𝑖 ) |𝜓𝑖 ⟩ = − (𝐻 (1) − 𝜀𝑖 ) |𝜓𝑖 ⟩,
with
𝐻

(1)

= (𝑇 + 𝑉ext )(1) + ∫𝛿 2 𝐸Hxc [𝜌(0) ] (1) ′𝜌 (𝑟 )𝑑𝑟 ′𝛿𝜌(𝑟)𝛿𝜌(𝑟 ′ )

and
(1)

𝜀𝑖

(0)

(0)

= ⟨𝜓𝑖 |𝐻 (1) |𝜓𝑖 ⟩.
57

3. DFPT for polarization
Let us study the polarization, which can be characterized by change in the polarization 𝑃𝛽caused by a uniform displacement of an atom 𝜅 in the direction 𝑅𝜅𝛼 : here uniform meansthat the atom in a cell is displaced together with its periodic images. The relevant quantity is
the Born effective charge defined by
∗
𝑍𝜅,𝛽𝛼
=Ω

𝜕𝑃𝛽
.
𝜕𝑅𝜅𝛼

Since the total energy can be described as𝐸[𝜌] − ΩP ⋅ ℰ
under the electric field ℰ induced by the polarization,𝑃𝛽 = −

1 𝜕𝐸[𝜌].
Ω 𝜕ℰ𝛽Therefore,
∗
𝑍𝜅,𝛽𝛼
=−

𝜕 2 𝐸[𝜌].
𝜕𝑅𝜅𝛼 𝜕ℰ𝛽This can be calculated by applying the DFPT. Since the 𝐸Hxc has not dependence on theatomic displacement and electric field, 𝑇 + 𝑉ext does not depend on the electric field, and the
kinetic energy operator is not affected by a uniform atomic displacement, the equation is
simplified as𝜕 2 𝐸[𝜌]
(ℰ )
(ℰ )
(ℰ )
(𝑅 )
(0)
(0)
= ⟨𝜓𝑖 𝛽 |(𝑉ext )(𝑅𝜅𝛼) |𝜓𝑖 ⟩ + ⟨𝜓𝑖 |(𝑉ext )(𝑅𝜅𝛼) |𝜓𝑖 𝛽 ⟩ + 2 ⟨𝜓𝑖 𝜅𝛼 |(𝑇 + 𝑉ext )(0) |𝜓𝑖 𝛽 ⟩𝜕𝑅𝜅𝛼 𝜕ℰ𝛽+∬

𝛿 2 𝐸Hxc [𝜌(0) ] (𝑅 )𝜌 𝜅𝛼 (𝑟)𝜌(ℰ𝛽 ) (𝑟 ′ )𝑑𝑟𝑑𝑟 ′ ,𝛿𝜌(𝑟)𝛿𝜌(𝑟 ′ )
(𝑅

)

(ℰ )

So, we need to solve the Sternheimer equation twice to get 𝜓𝑖 𝜅𝛼 and 𝜓𝑖 𝛽 , corresponding
to change in the Kohn-Sham orbitals with respect to the atomic displacement and thegeneration of the electric field.
About quantum espresso
Now let us use the environment
https://www.materialscloud.org/work/quantum-mobile
which will install a virtual machine (Virtual Box 5.2.8 + Ubuntu 16.04.4LTS) on which you
can run quantum espresso (v6.2.1) and wannier90 (2.1).
<First assignment: transfer of the benzene >58

Change your directory to wannier90-2.1 (depending on the version you installed) and then to
examples/example12. Then you will find input files for the calculation of benzene.
First, do an SCF calculation: In the SCF calculation, you get the effective Kohn-Shampotential, electron density and KS orbitals.
pw.x < benzene.scf > scf.out
Second, run wannier90 to generate a list of the required overlaps (written into thebenzene.nnkp file).
wannier90.x -pp benzene
Third, run pw2wannier90 to compute the overlap between Bloch states and the projections
for the starting guess (written in the benzene.mmn and benzene.amn files).
pw2wannier90.x < benzene.pw2wan > pw2wan.out
Finally, run wannier90 to compute the MLWFs.
wannier90.x benzene
Do this again by adding the lines to the benzene.win
restart = plotwannier_plot = true
wannier_plot_format = cube
wannier_plot_list = 2-4
Then you have files for the second to fourth MLWFs. By changing the digits, you can draw all
the orbitals.
1. Use xcrysden to draw all the MLWFs.
Do the same thing after you have inserted a line at the end of benzene.win
write_hr = true
Then you will get benzene_hr.dat
Focus on the third to fifth columns, say a, b, H. It indicates the transfer integral between theatoms a and b is H.
2. Show the value of transfer on a picture describing the structure of benzene. You need toshow only the largest ones, for example, only those whose magnitude is larger than 1 eV.
59

3. Describe scientific significance of the results.
<Second assignment: Born effective charge of Pb in PbTiO3 >Following the example04 of the quantum espresso, calculate the Born-effective charge of
Pb. Explain briefly procedures for the calculation, important parameters in the input file,result of the calculation. For this, use the ball and stick model obtained using, for example,
xcrysden. Also, show scientific significance of the result.

60