The 35th CMD workshop”Practice of ESopt”
Basic theory
Koichi KusakabeGrad. Sch. Eng. Sci, Osaka University1Born-Oppenheimer approximation◼ Based on difference in the mass, separation of the quantum

mechanical motion is done. (the error scales as (m/M)6/4.)
Cf. S. Takahashi, K. Takatsuka, J. Chem. Phys. 124, 144101 (2006).

Hˆ + Hˆelel −ion+ Hˆ ion −ion    ion (R I )= E + Eion −ion     ion (R I ) .

m: electron mass

M:nuclear mass

◼ Existence of stable materials allows us to fix positions for the masscenters of nuclei. Then,
Multiply R I  to the Schrödinger equation from the left,
N ion
Z I e2and define the external potential vext (r ) = 
for the electrons, we have,
r−R
I =1
I
22


ˆ (r )nˆ (r ') :3    †2
33 e :n3ˆ  (r ) ˆ  (r ) +  d rd r '+  d rvext (r )nˆ (r ) 
−   d r r − r'
 2m  =E.

The Coulomb gauge is fixed as the photon propagator in a steady state.
2Hamiltonian◼ Static state of electron system:◼

Represented by potential:

vext (ri )◼

e2V (ri , r j ) =.
ri − r j

Cf: Ehrenfest’s theoremrjrid p  =  − vext (r ) dt◼ The Hamiltonian◼ Velocities of electrons < the velocity of light,◼

In a steady material, internal gauge field behaves as theCoulomb potential.
N 2  2e2
 ri + H = − +  vext (ri ).
i =1  2m i , j ri − r ji =1N3

A practice◼ Let’s have an expression of the energy by vext.
NEext =   d rl   (rl ,  l ) vext (ri ) (rl ,  l )*3l

 l = 1i =1= N   d 3 rl  vext (r1 )  (rl ,  l )l2

 l = 1=  d 3 rvext (r ) N   d 3 rl   (r − r1 )  (rl ,  l )l2

 l = 1

N=  d rvext (r )   d rl   (rl ,  l )  (r − ri ) (rl ,  l )*33l

 l = 1i =1=  d 3 rvext (r )  nˆ (r )  =  d 3 rvext (r )n(r ).

◼ Variation of the energy w.r.t. vext, we have theexpression in the next page.
◼ [Problem] Derive the above expression in a different
smart manner.

4The density determining “order”◼ Energy as a functional of vext: E[vext].

◼ A functional derivative:E vext E vext + vext  − E vext  =  d rvext (r ) + O( 2 )vext (r )3( )=  d 3 rvext (r )n(r ) + O  2 .

(1)◼ The order parameter : n(r).

K.Kusakabe and I. Maruyama, J. Phys. A: Math. Theor. 44 (2011) 135305.
5

A theorem in quantum mechnics◼ Consider two potentials, v1(r) and v2(r), which

are different more than a constant from eachother.
◼ Let |1> and |2> be an eigen state of asystem by v1(r) and that by v2(r), respectively.
◼ Then, |1> cannot be an eigen state of v2(r)and vice versa.
Proof: The unique continuation theorem for the Poisson equation tells usthat the statement above is correct. See the next page.

6

A proof◼ When our system is treated by the Born-Oppenheimer

approximation, nuclei are treated as fixed point charges.
◼ v1 (and v2 also) follows the Poisson equation.
◼ Suppose that  is a solution of both of two Schrödinger equations
given by v1 and v2, where v1 and v2 are different from each othermore than a constant.
◼ We have v1 =(v2 +const.) .
◼ Let N be the number of electrons.
◼ We may find an open set S in R3 , where (r1+r,1, r2, 2,・ ・ ・ , rN, N) ≠ 0for r ∊ S.
◼ In S, we have v1(r) = v2(r) + const.
◼ By the unique continuation theorem, we conclude that v1(r) = v2(r)+ const in R3.
◼ This conclusion contradicts to the assumption.
7

Hohenberg-Kohn’s theorem◼ By the last theorem, we have a next inequality which is a contradiction.

E1 = 1 T + Vee 1 +  d 3 rnv1
= 1 H1 1  2 H1 2= 2 T + Vee 2 +  d 3 rnv1= 2 T + Vee 2 +  d 3 rnv2 +  d 3 rn(v1 − v2 )= 2 H 2 2 +  d 3 rn(v1 − v2 ) 1 H 2 1 +  d 3 rn(v1 − v2 )= 1 T + Vee 1 +  d 3 rnv2 +  d 3 rn(v1 − v2 )= 1 H1 1 = E1.

◼ This implies existence of “one-to-one correspondence between thepotential and the electron density”. (by Hohenberg and Kohn (1964))8A logic tells physics◼ Fixed positions of nuclei give a scalar potential vi.
◼ When v2 ≠ v1 + const., 1≠2 , where i is a ground state for a

problem by vi. (by the unique continuation theorem.)◼ In this case, n1 ≠ n2. (Proof by Hohenberg & Kohn.) Namely,◼ If “n1=n2=n” and “1≠2”, thenE1 = 1 T + Vee 1 +  d 3 rnv1 = 1 H1 1  2 H1 2 = 2 T + Vee 2 +  d 3 rnv1= 2 T + Vee 2 +  d 3 rnv2 +  d 3 rn(v1 − v2 ) = 2 H 2 2 +  d 3 rn(v1 − v2 ) 1 H 2 1 +  d 3 rn(v1 − v2 ) = 1 T + Vee 1 +  d 3 rnv2 +  d 3 rn(v1 − v2 ) = 1 H1 1 = E1.
Thus, E1  E1 , which is a contradiction.

◼ Thus, we have a one to one correspondence between v and n.
◼ Namely, n gives v, the Hamiltonian, and the wave function of

the ground state  with energy E[n].
9Variational energyE0 = min F [n] +  n(r )  vext (r )drn= min min  Tˆ + Vˆee  +  n(r )  vext (r )dr = min Evext [n].
n →n

n1 (r )n2 (r )n(2)A subset of w.f.

Evext1 [n]Evext 2 [n]n3 (r )This definition of the density functional theory does not rely on HK.
10

’’’

Definition of Exc[n].

11Derive differential equations◼ First, set up a differentiable energy functional.

Exc nE xc n

GT  GT  

◼ Derive a determining equations, which are to be

self-consistent equations.
 G   − E (   − 1) = 0,T

Tˆ +  d rv (r )nˆ(r )  = E  ,3effn(r ') E xc nveff (r ) =  d r '++ vext (r ),r − r ' n(r )3n(r ) =   nˆ (r )  .
12The local density approximation I.

13

The local density approximation II.

dddd

14

The local density approximation III.

15

Correlation energy density for LDA
Let’s introduce rs byNe 4 3 = 1 /
rs  We may divide exc as, e xc (rs ) = e x (rs ) + e c (rs ).
V 3 Vosko, Wilk & Nusair (1980)ec

Perdew & Zunger (1981)

Hedin & Lundqvist (1971)

Nozières & Pines (1958)
Wigner (1934)

16

The plane-wave expansion I.

17

The plane-wave expansion II.

‘• Dimension of the matrix can be O(103) or O(104).
• The diagonalization is often performed using the Housholder method or theconjugate gradient method (CG).
However,• the process to find the ground state is regarded as an optimization process forthe energy functional in the function space.
The Car-Parrinello method (Conceptually different idea)
18

The pseudopotential I.

Ca 4p orbital
All electron calc.

Ca 4s orbitalPseudization calc.

All electron calc.
Pseudization calc.

19

The pseudopotential II.

20

The pseudopotential III.

21

The pseudopotential IV.

22

The pseudopotential V.

23

Application for a metal◼ The band structure of Al (FCC)
The plane-wave expansion with a pseudo-potential◼ GGAPBE96◼

24

Covalent crystals I.
◼ Diamond and Graphite as covalent crystals.

Diamond structure

Graphite structure25Band structure of cubic diamond
* bands(anti-bondingbands)

 bands(bonding bands)

An energy gap appears and the system is a wide-gap
semiconductor.
26

Bonding charge in hex-diamond◼ In a covalent crystal,

we can see charge
density of electrons
at each bondconnection.
◼ Yellow object
represents charge
density and whitespheres are carbons.

27Band structure of graphite* bands(anti-bondingbands)
* bands(anti-bonding) bands(bonding bands)
 bands(bonding bands)

The -band is half-filled and there are small
Fermi pockets both for electrons and holes. (Semimetal)
28

Bonding charge in graphite◼ Bonding charge

comes from electrons.
◼ This system is asemimetal wherethe Fermi surface
is made of bands.

29

Techniques in the molecular dynamics◼ Constant-Temperature Scheme (Nose` & Hoover)◼ Average of the kinetic energy is made constant byintroduction of interaction between the system and an
imaginative heat bath.
◼ An equation of motion with dissipative term written by an
artificial variable is solved.
◼ Constant-Pressure Scheme (Parrinello & Rahman)◼ The volume (cell parameters) of the system is made avariable.
◼ MD and structural optimization is performed using theenthalpy, H=E+PV.

Realization of realistic situation

30

The Hellmann-Feynman forcei31

The internal stress I

32

The internal stress II

33

The internal stress III34

The calculation scheme of FPMD at a constant pressure

Structuraloptimization

Input: Unit Cell{a0}
Atoms
{RI}
External Pressure Pexta = (1 +e ) a0

a = (1 + e)a0

Wave Functions{i }HF Force
Quantum Stress{F}{ }Convergence check{RI} {e} update

Moleculardynamics{Kohn-Sham eq.
by CG method

Molecular Dynamics
using {F} { }{RI} {e} update

Optimized { i} {RI} {e}35Nose-Hoovermethod (Introduction)´◼ For constant temperature calculations, we want

to reproduce the canonical distribution.

Microcanonical ensemble(E=const. N=const.)

Contact withartificial heat-bathCanonical ensemble(T=const. N=const.)

The equation for continuity

Classical dynamics

Liouville’s equation

The Liouville theorem

Nose-Hooverdynamics´

Canonical ensemble simulated by Nose´
Hoover dynamics

Nose-Hoover
chain method´Heat-bath is represented by a chain of variables h.

Parrinello-Rahman method(Introduction)◼ To determine lattice structure automatically, or tosimulate structural transformation, we want to simulatea system with a variable unit cell.

bCompressionb’
ri’=xia’+hib’+zic’

ri=xia+hib+zica’aLow pressure

High pressure

Definitions

Parrinello-Rahman LagrangianEquations of motion for CPMD

This set of equations allows us to perform the molecular dynamics
at a constant pressure.

The equations of motion for CP-FPMD

tot

48

Graphite-diamond transformation in a FPMD
simulation
At ambient pressure

At high pressure49Formation of new bondings at the transformation• White objects are carbon
atoms and yellow isosurfaces represent charge
density of electrons.
• We see new bonding
represented by bonding
charge between graphitelayers.
• Sliding of layers
occurs due to formationof sp3 bond connections.

50Graphitic structure stable in 40GPa range• G-ball has a concentric
structure made of graphitic
sheets.
• The ball does not show
any structuraltransformation up to ～40GPa.
• The conductance of theball becomes high at highpressure.
TEM image of the graphite ballby F. Kokai et al. (2001).
51

Lattice constants of graphitic system in high pressure

1. Calculation has been done
for hexagonal graphite withAB stacking of graphenelayers.
2. Since we assume thegraphite structure in calc.,
no phase transformationis seen and the structureis meta-stable up to 60GPa.

Exp: by Nakayama et al.
52What is ESopt?
◼ Based on “opt” by ISSP, Univ. of Tokyo, werevised the program to have Esopt.
◼ Plane wave expansion method.
◼ Pseudo potential (Norm-conserving P.P.)◼ Exchange correlation: PW91

◼ Some merits of “opt” succeed to ESopt.
◼ Readability of the source code
◼ Optimization ofthe wavefunction by the original

CG method53

Characteristics of basis sets for DFT calculation••

Plane-wave expansion method with pseudo-potentials• Since plane waves are independent of position of atoms, the result isaccurate with respect to the valence electrons.
• Accuracy of the calculation is determined by the maximum energy of planewaves.
• The kinetic energy is diagonal in the Fourier space, while the potential
energy is diagonal in the real space. FFT is used to connect two spaces.
• The Hellmann-Feynman force and the quantum stress are easily obtained.
FLAPW (Full-potential linearlized augmented plane wave)• The wave functions in an atomic sphere are expanded in spherical waves.
Otherwise, they are written in the plane waves.
• Accuracy is determined by number of spherical waves and the maximum
energy of plane waves.
• Less ambiguity compared to the pseudo-potential method.
• Pulay force has to be evaluated.
54