Integrated Computational Materials

Engineering Education

Calculation of Equation of State Using

Density Functional Theory

Mark Asta1, Katsuyo Thornton2, and Larry Aagesen3

1
Department of Materials Science and Engineering, University of California, Berkeley
2
Department of Materials Science and Engineering, University of Michigan, Ann Arbor
3
Idaho National Laboratory

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Purposes of Density Functional Theory Module

• Understand fundamentals of Density Functional

Theory (DFT)
• Apply DFT to calculate:
– Equilibrium lattice constant
– Bulk Modulus
– Components of elastic constant tensor
• Understand how to check for convergence of
results

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Equation of State
A Probe of Interatomic Interactions

Schematic Energy vs. Diamond Cubic

Volume Relation Structure of Si
Energy
per atom

a
a
http://www.e6cvd.com/cvd/page.jsp?
pageid=361
Volume per atom (=a3/8 for Si)

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Equation of State
What Properties Can we Learn from It?

Pressure versus Volume Relation

¶E Given E(V) one can compute P(V) by taking derivative
P =-
¶V
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 B related to curvature of E(V) Function
¶V ¶V

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Generalize to Non-Hydrostatic Deformation
Example of Uniaxial Deformation

Lz Lz(1+)

Ly Ly
Lx Lx

Definition of Deformation In Terms of Strain:

e33 =e (All other strains are zero)

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Linear-Elasticity for Single Crystals
3 3
General form of Hook’s Law s =C e =
(Linear Elasticity):
ij ijkl kl å å Cijklekl
k=1 l=1

Stress Elastic Strain

Tensor Constant Tensor
Tensor
Voigt 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5, 12 → 6
6
Notation: s i =Cije j =å Cijej
j=1
1
Elastic Energy: E = Cijeiej
2
In example from e33 =e (All other strains are zero)
previous slide: s 33 =C3333e33 º C11e
s 22 =C2233e33 º C12e
Note: for cubic crystal C11=C22=C33, C12=C13=C23
DFT Module Review, The 7th 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 7th Summer School for Integrated Computational Materials Education
 Total Energy in Density Functional Theory

Ne

Electron Density n(r) =- e å f i (r)

2

i=1
Electron Wavefunctions fi (r)
Potential Electrons Feel from Nuclei Vext (r)
  Energy
Exchange-Correlation Exc [n(r)]
Form depends on whether you use Local Density Approximation (LDA)
or Generalized Gradient Approximation (GGA)
DFT Module Review, The 7th 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)úf i (r) =eifi (r)
ê
ë 2me r - r' ú
û
dE xc [n(r)]
Exchange-Correlation Potential Vxc (r) º
dn(r)
Ne

Electron Density n(r) =- e å f (r)i

2

i=1
 
Note: i depends on n(r) which depends on i 
Solution of Kohn-Sham equations must be done iteratively
 
DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Self-Consistent Solution to DFT Equations
Input Positions of Atoms for a Given 1. Set up atom positions
Unit Cell and Lattice Constant
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
different compare output and threshold, mix output and input density and
input charge densities go to step 2
6. Quit when energy from input and output
Energy for Given same
Lattice Constant densities agree to within prescribed
tolerance (e.g., 10-5 eV)

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 7th 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)
a a3 = a (0, 1/2, 1/2)
Basis Atom Positions
a 000
a ¼¼¼

All atoms in the crystal can be obtained by adding integer

multiples of unit cell vectors to basis atom positions

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Electron Density in Crystal Lattice

Unit-Cell Vectors
a1 = a (-1/2, 1/2 , 0)
a
a2 = a (-1/2, 0, 1/2)
a3 = a (0, 1/2, 1/2)
a
a
Electron density is periodic with periodicity given by R uvw
n ( r) =n ( r + R uvw )
Translation Vectors: R uvw =ua1 + va 2 + wa 3

DFT Module Review, The 7th 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 ()
DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Why?
Wavefunctions in a Crystal Obey Bloch’s Theorem

For a given band 

b b
f ( r) =exp ( ik ×r) u ( r)
k k

b b b
Where k ( r ) is periodic in real space: uk ( r ) =uk ( r + R uvw )
u

Translation Vectors: R uvw =ua1 + va 2 + wa 3

The envelope function represents delocalized distribution of electron density

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 Representation of Electron Density
b b
f ( r) =exp ( ik ×r) u ( r )
k k
occ
Ne 2 d 3k
n(r) =- eå f i (r)
2
n(r) =- eå ò f kb ( r )
i=1 b
W BZ
W BZ
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 f
2
b
kj ( r)
b 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 7th Summer School for Integrated Computational Materials Education
 Representation of Wavefunctions
Fourier-Expansion as Series of Plane Waves

For a given band: f kb ( r) =exp ( ik ×r) ukb ( r)

b b b
Recall that uk ( r ) is periodic in real space: uk ( r ) =uk ( r + Ruvw )

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

f kb ( r) =å ûkb ( G) exp éëi ( G + k) ×rùû

G

In practice the Fourier series is truncated to include all G for which:

2 2
( G + k) < E cut
2m

Another parameter that needs to be checked for convergence is

the “plane-wave cutoff energy” Ecut
 
DFT Module Review, The 7th 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 7th 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 7th 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 7th 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)
a a3 = a (0, 1/2, 1/2)
Basis Atom Positions
a (Fractional Coordinates)
a 000
¼¼¼

DFT Module Review, The 7th Summer School for Integrated Computational Materials Education
 DFT Module
• For problem 2 you will impose a homogeneous tensile
strain () along the [001] axis (see slide 4)
• Such a strain results in the x3 coordinate of all atoms
changing to x3*(1+)
• 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+)/2)
a a3 = a (0, 1/2, (1+)/2)
Basis Atom Positions
a (Fractional Coordinates)
a 000
¼¼¼
DFT Module Review, The 7th Summer School for Integrated Computational Materials Education