DFT Module Review - Summer School for Integrated Computational

Integrated Computational Materials
Engineering Education
Calculation of Equation of State Using
Density Functional Theory
Mark Asta1, Larry Aagesen2 and Katsuyo Thornton2
1Department
of Materials Science and Engineering, University of California, Berkeley
2Department of Materials Science and Engineering, University of Michigan, Ann Arbor
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Equation of State
A Probe of Interatomic Interactions
Schematic Energy vs.
Volume Relation
Diamond Cubic
Structure of Si
Energy
per atom
a
a
a
http://www.e6cvd.com/cvd/page.jsp?pagei
d=361
Volume per atom (=a3/8 for Si)
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Equation of State
What Properties Can we Learn from It?
Pressure versus Volume Relation
P=-
¶E
¶V
Given E(V) one can compute P(V) by taking derivative
Recall 1st Law of Thermo: dE = T dS - P dV and consider T = 0 K
Equilibrium Volume (or Lattice Constant)
Volume corresponding to zero pressure = Volume where slope of E(V) is zero
≈ Volume measured experimentally at P = 1 atm
Bulk Modulus
¶P
¶2 E
B = -V
=V 2
¶V
¶V
B related to curvature of E(V) Function
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Generalize to Non-Hydrostatic Deformation
Example of Uniaxial Deformation
Lz
Lz(1+e)
Ly
Lx
Ly
Lx
Definition of Deformation In Terms of Strain:
e33 = e
(All other strains are zero)
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Linear-Elasticity for Single Crystals
3
3
General form of Hook’s Law s = C e =
Cijklekl
å
å
ij
ijkl kl
(Linear Elasticity):
k=1 l=1
Stress
Tensor
Elastic
Constant
Tensor
Strain
Tensor
Voigt Notation: 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6
6
s i = Cije j = åCije j
j=1
Elastic Energy:
In example from
previous slide:
1
E = Cijeie j
2
e33 = e
(All other strains are zero)
s 33 = C3333e33 º C11e
s 22 = C2233e33 º C12e
Note: for cubic crystal C11=C22=C33, C12=C13=C23
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Equation of State
How to Calculate from Density Functional Theory
Formulation: for a given arrangement of nuclei defined by the
lattice constant, crystal structure, and non-hydrostatic strains,
compute the total energy corresponding to the optimal
arrangement of the electron density
Theoretical Framework: Quantum mechanical calculation of
energy of electrons and nuclei interacting through Coulomb
potential
Practical Implementation: Density functional theory
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Total Energy in Density Functional Theory
Ne
å f (r)
Electron Density n(r) = -e
Electron Wavefunctions
Potential Electrons Feel from Nuclei
Exchange-Correlation Energy
fi (r)
2
i
i=1
Vext (r)
Exc [n(r)]
Form depends on whether you use Local Density Approximation (LDA)
or Generalized Gradient Approximation (GGA)
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Kohn-Sham Equations
Schrödinger Equation for Electron Wavefunctions
2
é
ù
n(r
')
Ñi2 +Vext (r) + ò
d 3r '+Vxc (r)úfi (r) = eifi (r)
êr - r'
êë 2me
úû
Exchange-Correlation Potential
dE xc [n(r)]
Vxc (r) º
dn(r)
Ne
å f (r)
Electron Density n(r) = -e
2
i
i=1
Note: fi depends on n(r) which depends on fi 
Solution of Kohn-Sham equations must be done iteratively
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Self-Consistent Solution to DFT Equations
1.
Set up atom positions
2.
Make initial guess of “input” charge density
(often overlapping atomic charge densities)
guess charge density
3.
Solve Kohn-Sham equations with this input
charge density
compute effective
potential
4.
Compute “output” charge density from
resulting wavefunctions
compute Kohn-Sham
orbitals and density
5.
If energy from input and output densities
differ by amount greater than a chosen
threshold, mix output and input density and
go to step 2
6.
Quit when energy from input and output
densities agree to within prescribed
tolerance (e.g., 10-5 eV)
Input Positions of Atoms for a Given
Unit Cell and Lattice Constant
different
compare output and
input charge densities
Energy for Given
Lattice Constant
same
Note: In this module the positions of atoms are dictated by symmetry. If this is not the
case another loop must be added to minimize energy with respect to atomic positions.
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Implementation of DFT for a Single Crystal
Crystal Structure Defined by Unit Cell Vectors and
Positions of Basis Atoms
Example: Diamond Cubic Structure of Si
Unit Cell Vectors
a1 = a (-1/2, 1/2 , 0)
a2 = a (-1/2, 0, 1/2)
a3 = a (0, 1/2, 1/2)
a
a
a
Basis Atom Positions
000
¼¼¼
All atoms in the crystal can be obtained by adding integer
multiples of unit cell vectors to basis atom positions
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Electron Density in Crystal Lattice
Unit-Cell Vectors
a1 = a (-1/2, 1/2 , 0)
a2 = a (-1/2, 0, 1/2)
a3 = a (0, 1/2, 1/2)
a
a
a
Electron density is periodic with periodicity given by Ruvw
n ( r) = n ( r + Ruvw )
Translation Vectors: Ruvw = ua1 + va 2 + wa3
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Electronic Bandstructure
Example for Si
Brillouin Zone
Bandstructure
http://en.wikipedia.org/wiki/Brillouin_zone
http://de.wikipedia.org/wiki/Datei:Band_structure_Si_schematic.svg
Electronic wavefunctions in a crystal can be indexed by
point in reciprocal space (k) and a band index (b)
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Why?
Wavefunctions in a Crystal Obey Bloch’s Theorem
For a given band b
f ( r) = exp (ik × r) u ( r)
b
k
b
k
b
b
b
u
Where k ( r) is periodic in real space: uk ( r) = uk ( r + Ruvw )
Translation Vectors: Ruvw = ua1 + va 2 + wa3
The envelope function represents delocalized distribution of electron density
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Representation of Electron Density
f ( r) = exp (ik × r) u ( r)
b
k
Ne
b
k
n(r) = -eå fi (r)
occ
n(r) = -eå ò
2
b
i=1
WBZ
fkb ( r)
2
d 3k
WBZ
Integral over k-points in first Brillouin zone
In practice the integral over the Brillouin zone is replaced
with a sum over a finite number of k-points (Nkpt)
N kpt
n(r) » -eå å w j fk j ( r)
b
b
2
j=1
Band occupation (e.g., the Fermi function)
One parameter that needs to be checked for numerical
convergence is number of k-points
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Representation of Wavefunctions
Fourier-Expansion as Series of Plane Waves
For a given band:
fkb ( r) = exp (ik × r) ukb ( r)
Recall that uk ( r) is periodic in real space: uk ( r) = uk ( r + Ruvw )
b
b
b
ukb ( r) can be written as a 3D Fourier Series:
ukb ( r) = åukb ( Glmn ) exp (iGlmn × r)
lmn
Glmn = la1* + ma*2 + na*3
where the a *i are primitive reciprocal lattice vectors
a1* × a1 = 2p
a1* × a 2 = 0
a1* × a 3 = 0
a*2 × a1 = 0
a*2 × a 2 = 2p
a*2 × a 3 = 0
a*3 × a1 = 0
a*3 × a 2 = 0
a*3 × a 3 = 2p
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Recall Properties of Fourier Series
Black line = (exact) triangular wave
Colored lines = Fourier series
truncated at different orders
http://mathworld.wolfram.com/FourierSeriesTriangleWave.html
General Form of Fourier Series:
For Triangular Wave:
Typically we expect the accuracy of a truncated Fourier series to
improve as we increase the number of terms
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Representation of Wavefunctions
Plane-Wave Basis Set
For a given band
fkb ( r) = exp (ik × r) ukb ( r)
Use Fourier Expansion
fkb ( r) = åuˆkb ( G) expéëi ( G + k) × rùû
G
In practice the Fourier series is truncated to include all G for which:
2
2m
( G + k)
2
< E cut
Another parameter that needs to be checked for convergence is
the “plane-wave cutoff energy” Ecut
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
Examples of Convergence Checks
Effect of Ecut
Effect of Number of k Points
Note: the different values of kTel
corresponds to different choices for
occupation function (wj in slide 14)
http://www.fhi-berlin.mpg.de/th/Meetings/FHImd2001/pehlke1.pdf
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
DFT Module
• Problem 1: Calculate equilibrium volume and bulk
modulus of diamond cubic Si using Quantum Espresso
on Nanohub (http://nanohub.org/)
o Outcome 1: Understand effect of numerical parameters on
calculated results by testing convergence with respect to
number of k-points and plane-wave cutoff
o Outcome 2: Understand the effect of the theoretical model
for exchange-correlation potential on the accuracy of the
calculations by comparing results from Local Density
Approximation (LDA) and Generalized Gradient
Approximation (GGA) with experimental measurements
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
DFT Module
• Problem 2: Calculate the single-crystal elastic constants
C11 and C12
o Outcome 1: Understand how to impose homogeneous
elastic deformations in a DFT calculation
o Outcome 2: Understand the effect of the theoretical model
for exchange-correlation potential on the accuracy of the
calculations by comparing results from Local Density
Approximation (LDA) and Generalized Gradient
Approximation (GGA) with experimental measurements
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
DFT Module
• For problem 1 you will make use of the unit cell for
diamond-cubic Si shown below. You will vary only the
lattice constant a.
Unit Cell Vectors
a1 = a (-1/2, 1/2 , 0)
a2 = a (-1/2, 0, 1/2)
a3 = a (0, 1/2, 1/2)
a
a
a
Basis Atom Positions
(Fractional Coordinates)
000
¼¼¼
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education
DFT Module
• For problem 2 you will impose a homogeneous tensile
strain (e) along the  axis (see slide 4)
• Such a strain results in the x3 coordinate of all atoms
changing to x3*(1+e)
• This homogeneous deformation can be represented by
changing the unit cell vectors as follows:
Unit Cell Vectors
a1 = a (-1/2, 1/2 , 0)
a2 = a (-1/2, 0, (1+e)/2)
a3 = a (0, 1/2, (1+e)/2)
a
a
a
Basis Atom Positions
(Fractional Coordinates)
000
¼¼¼
DFT Module Review, The 3rd Summer School for Integrated Computational Materials Education