Report

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 [001] 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