CHEM 938: Density Functional Theory

CHEM 938: Density Functional Theory condensed-phase systems February 11, 2010

Pseudopotentials a huge number of planewaves are needed to describe core states, which have a very high curvature • consider the all-electron wavefunction (AE) in the plot • high curvature in core region requires high values of G in planewave expansion (lots of basis functions) • core electrons (usually) don’t matter for chemical applications, so let’s replace them with a pseudopotential (PS) • wavefunction and potential in valence region match, but core region is less curved so that lower values of G can be used distance from nucleus

Pseudopotentials we have to consider how core electrons affect valence electrons Fock operator: these terms describe: • Coulomb attraction between valence electron a and the nuclei • Coulomb repulsion between valence electron a and the core electrons • exchange repulsion between valence electron a and the core electrons

Pseudopotentials we have to consider how core electrons affect valence electrons effective core potentials combine the effects of the core electrons plus nuclear charges into a single potential energy term each atom gets an ECP, which accounts for core electrons and nuclear charges eliminates an explicit treatment of core electrons with basis functions

Pseudopotentials Advantages of Pseudopotentials: • reduce significantly the number of basis functions used in a calculation • can easily introduce relativistic effects for core electrons • core electrons of heavy elements move at near the speed of light, where relativistic effects matter • with ECPs we can include relativistic effects without using a relativistic Hamiltonian Disadvantages of Pseudootentials: • sometimes core electrons matter • defining core states can be ambiguous • are 3p electrons included in valence or core for first row transition metals? can be hard to construct (at least good ones can be) transferability is not guaranteed NOTE:pseudopotentialsare also called effective core potentials

Projector Augmented Wavefunction Methods what if you want to keep the core electrons? • transform from pseudo wavefunction to all-electron wavefunction pseudo orbital all electron orbital transformation operator • transformation operator is defined as: projector function pseudo atomic orbital i on atom a all electron atomic orbital i on atom a

Projector Augmented Wavefunction Methods what if you want to keep the core electrons? • transform from pseudo wavefunction to all-electron wavefunction difference between real and pseudo atomic orbitals + - = difference between real and pseudo atomic orbitals pseudo wavefunction all electron wavefunction

Projector Augmented Wavefunction Methods PAWs employ some approximations and have some benefits approximations: • core orbitals are not optimized (frozen-core) • truncated plane-wave expansion in basis set • truncated wave expansion in core (augmentation) benefits: • incorporate explicitly effects of core electrons • retain kinetic energy cutoffs used with ‘ultrasoft’ pseudopotentials (Ecut ~ 30Ry) • only a few projectors per atoms and angular momentum channel • better transferability I would advise the use of PAWs (or other augmentation schemes) if they are available in the code you are using

The Ingredients periodic boundary conditions • need a simulation cell defined by lattice vectors a, b, c • need to place the atoms at appropriate positions within that cell pseudopotentials or PAWs • these may come with the program you are using (but don’t trust the at face value) • you may have to construct them yourself (tough work, but sometimes necessary) planewave basis sets • you simply define a value of Ecut • this will generally depend on the atoms in the system, as well as the pseudopotentials you are using • this is something you just have to test out k-point grid • you define the number of divisions along each reciprocal lattice vector • generally dependent on the system (but requires testing) • for non-metals ~0.05 Å-1 spacing is usually sufficient, denser grids are needed for metals

Methods and Accuracy systematic studies of the performance of DFT methods for many material properties are lacking methods: • virtually all electronic structure calculations of materials are performed with DFT at the LDA or GGA levels • hybrid functionals, Hartree-Fock and ab initio methods require the evaluate of the exact exchange integral • evaluating exact exchange scales as K4 and with planewave basis sets, K ~ two to three orders of magnitude greater than with atom-centered basis functions accuracy: • results are usually quite good, but some problematic systems have been identified • unlike molecular systems, LDA functionals can do a good job for some materials • this is because LDA functionals are based on a homogeneous electron gas, which is very much like a metal • however, GGAs usually perform better and should be used when possible

Quick Tests consider convergence on total energies of chromia (Cr2O3) with planewave cutoffs and k-point grid converged to ~1 meV/atom converged to ~1 meV/atom numbers increase instead of decrease because of peculiarities of VASP software Ecut = 450 eV is sufficient 2x2x1 k-point grid is sufficient

Geometry Optimization just as with molecules, getting optimized structure is important for materials geometry optimization to minima: • works just like with molecules • also called relaxation • only advantage is there are no Pulay forces with planewave basis sets • still need derviatives with respect to grids used to evaluate Exc: 0

Transition State Optimization just as with molecules, getting optimized structure is important for materials geometry optimization to transition states: • usually perform nudged-elastic band calculations • guess a series of structures leading from reactants to products (usually a linear interpolation between these states) • connect these structures with artificial springs • optimize each structure while minimizing the action along the total path to obtain the minimum energy path • variants exist to allow points to move to locate the maximum (TS) along the reaction pathway

Cell Optimization cell optimization is also important for materials we also need to optimize the cell size and shape: • two common ways to do this: 1. vary lattice constants to minimize energy • approach works well when crystal structure is well-defined and affected by only one parameter • e.g. cubic cell can be varied only changing length, because all vectors have same length and angles between vectors are all 90° lattice vector length of Si in fcc cell • more complicated structures, involve more combinations of parameters to optimize and are not best treated through manual optimization • optimal length ~5.47 Å, expt. = 5.43 Å

Cell Optimization cell optimization is also important for materials we also need to optimize the cell size and shape: • two common ways to do this: 2. vary lattice vectors to minimize stress • stress is force over area • in materials, ∂x corresponds to a change in the lattice vector (strain), and the area A is defined in terms of lattice vectors • area = |a X b| • direction of a X b is normal to plane spanned by a and b • different combinations of direction along which F is applied and the area over which it is applied give rise to the stress tensor, σij

Cell Optimization cell optimization is also important for materials stress tensor: • consider a cubic cell with lattice vectors a, b, and c aligned with the x, y, and z axes, respectively i = direction of vector normal to area over which force is applied j = direction along which force is applied F • F along z direction A • A spanned by x and y axes (x-y plane) • surface normal corresponds to z direction z • this gives σzz y • corresponds to a force along z resulting in compression along z x • if forces extends material, this is referred to as tension

Cell Optimization cell optimization is also important for materials stress tensor: • consider a cubic cell with lattice vectors a, b, and c aligned with the x, y, and z axes, respectively i = direction of vector normal to area over which force is applied j = direction along which force is applied F • F along z direction • A spanned by x and z axes (x-z plane) A • surface normal corresponds to z direction z • this gives σyz y • corresponds to a force along z resulting in shear along y x • this is called a shear stress

Cell Optimization cell optimization is also important for materials stress tensor: • in general, the stress tensor is a 3x3 second rank tensor • diagonal elements are related to the pressure through: • these cause changes in the volume of the cell, but not the shape

Cell Optimization cell optimization is also important for materials stress tensor: • in general, the stress tensor is a 3x3 second rank tensor • off-diagonal elements correspond to shear • these cause changes in the shape of the cell, but not the volume • the shear stresses are symmetric: • otherwise, the cell would rotate

Cell Optimization cell optimization is also important for materials we also need to optimize the cell size and shape: • when the cell is at it’s optimal shape an size • to arrive at this state, we treat the components of the lattice vectors as variables that can be optimized just like atomic coordinates in geometry optimizations • however, calculating the stress depends on the derivative of the energy with respect to the cell vectors • once again, we use the Hellmann-Feynman theorem for this purpose

Cell Optimization cell optimization is also important for materials do not equal zero because planewave basis functions depend on lattice vectors • basically, if we change a, b, or c, then we alter the curvature of the plane wave • ability of planewave basis set to represent wave function depends on curvature of plane waves • so, changing cell size alters accuracy of calculation

Cell Optimization change cell size alters accuracy of basis set compressed cell original cell two ways to deal with this: use a really high kinetic energy cutoff so that has such high curvature that small changes don’t matter remove or add planewaves when cell size changes to ensure curvature is consistent

Cell Optimization convergence of pressure with Ecut • check convergence of stress tensor to see if Ecut is high enough to handle changes in cell size • pressure (and other stress components) should be converged to about 0.5 kbar pressure converged at Ecut = 550 eV • usually requires an Ecut that is ~30% higher than that needed when the cell size remains constant energy converged at Ecut = 450 eV • NOTE: this is only a problem if you keep the number of plane waves constant. That is, if you do two separate calculations on cells of different size, you should use the same Ecut for consistency convergence of pressure with kinetic energy cutoff

Cell Optimization we can also remove/add planewaves as needed to maintain curvature • curvature of the planewave is determined directly by the highest value of G contained in the plane wave expansion and has a kinetic energy of: • so define a value of Gmax that should be maintained throughout the calculation • then add or remove planewave using a step function to maintain to only use planewaves up to Gmax in the calculation • A and σ are parameters determining the width and height of the step function • this approach allows the use of low kinetic energy cutoffs, but requires testing to figure out appropriate values of A and σ M. Bernasconi, et al, J. Phys. Chem. Solids, 56, 501 (1995)

Electronic Structure we are often interested in band structure and density of states when studying materials • we can obtain the band structure by evaluating the wavefunction at several k-points • we can evaluate the density of states by integrating the band structure over the Brillouin zone

Band Structure evaluating the wavefunction at several k-points gives the band structure band structure of aluminum band structure of silicon • metal with no band gap • semiconductor with indirect band gap of 1.6 eV

Band Structure to properly model metals, you need to use dense k-point grids • a 1x1x1 grid (only Γ point used in calculation) would predict a large band gap • to capture metallic behavior, the L, X, W, and U points have to be included in the k-point grid • number of grid points needed can be an order of magnitude greater than for insulators band structure of aluminum

Band Structure to properly model metals, you probably also need to consider fractional occupation numbers • with no band gap, electrons can enter conduction band at thermal energies • EF = Fermi energy: the energy of the highest occupied state assuming that each state contains 1 electrons (or 2 if closed-shell) • often accounted for by using fractional occupations other smearing schemes exist, with various advantages and disadvantages

Band Structure k-points can also be a problem with semi-conductors • band gap is sufficiently large to prevent the need for smearing • band energies change rapidly in k-space, so you need a large k-point grid band structure of silicon • semiconductor with indirect band gap of 1.6 eV • correct structure only obtained with denser k-point grids

Density of States band structure plots are complicated to analyze – usually we condense them into a density of states plot • density of states tells us how many states are present at a particular energy • can determine the band gap • Fermi level is usually shifted to 0 eV • can project onto atomic orbitals to identify states on different atoms Egap DOS silicon chromia (Cr2O3)

Strongly Correlated Electron Systems some electronic structures are not described properly by DFT EFermi n(E) Energy • partially filled states  conductor • orbitals hybridize • electrons delocalize With strong correlations: EFermi n(E) • strong on-site Coulomb repulsion prevents delocalization X e- Energy e- X • Mott-Hubbard insulator these include materials ranging from superconductors to metal oxides

Electronic Structure of Chromia Experimentally: d states d states charge transfer/Mott-Hubbard insulator p states • valence band = Cr d + O p states • conduction band = Cr d states • band gap = 3.4 eV E With LDA: • valence band = Cr d states • conduction band = Cr d states • band gap = 0.8 eV inconsistent with experiment ultimately a self-interaction problem

delocalized state (r) • decreased self-interaction energy r Strongly-Correlated Electron Materials localized state (r) • large self-interaction energy DFT r For strongly-correlated electron materials: • ‘correct’ density  large self-interaction error in DFT • electrons delocalize to reduce artificial electron-electron repulsion • delocalization: • qualitatively alters electronic structure • shifts orbital energies • lowers band gap

average Coulomb repulsion between localized electrons average exchange repulsion between localized electrons DFT+U treat on-site interactions with a Hartree-Fock-like potential Dudarev et al., Phys. Rev. B, 57, 1505 (1998)

average Coulomb repulsion between localized electrons average exchange repulsion between localized electrons DFT+U treat on-site interactions with a Hartree-Fock-like potential Energy Functional: on-site density matrix of d electrons density of delocalized and localized states Dudarev et al., Phys. Rev. B, 57, 1505 (1998)

DFT+U Energy Functional: Effects on the electronic structure: • favors integer occupations  counters delocalization • opens band gap ni = 0 i ni = 0 ni = 1 ni = 1 DFT DFT+U Dudarev et al., Phys. Rev. B, 57, 1505 (1998)

Effect of DFT+U on Electronic Structure Does DFT+U improve the electronic structure? LDA LDA+U, (U-J) = 7.7 eV • no mixing of Cr d- and O p-states in valence bands • mixing of Cr d- and O p-states in valence bands • band gap = 0.8 eV • band gap = 3.8 eV inconsistent with experiment agrees with experiment

Ab initio Evaluation of U and J treat U and J as average HF Coulomb and exchange interactions between localized electrons ‘localized’ states in Cr2O912- Mosey, Liao & Carter, J. Chem. Phys. B, 129, 014103 (2008)

Ab initio Evaluation of U and J U and J are obtained through cluster calculations U-J = 3.2 eV U-J = 4.3 eV Mosey, Liao & Carter, J. Chem. Phys. B, 129, 014103 (2008)

Performance performance is pretty good with LDA functionals • LDA+U matches experiment quite well with ab initio values of U and J • GGA+U approach doesn’t seem to work as well: probably due to an inexact removal of double-counting terms • other types of DFT+U implementations exist, and may use slightly different interpretation of U and J Mosey, Liao & Carter, J. Chem. Phys. B, 129, 014103 (2008)

Recap of Scenarios here’s a general outline of requirements for different materials/scenarios metals: • dense k-point grids • electron smearing for fractional occupations semi-conductors: • dense k-point grids • no electron smearing insulators: • coarse k-point grids • no electron smearing strong correlations: • coarse k-point grids • no electron smearing • DFT+U, probably with LDA functional stress tensor: • high cutoffs • or addition/removal of planewaves

- (1120) b a c b c a Surfaces the properties of surfaces are of interest in materials science surfaces are made by cleaving bulk structures: - (0001) (0112) possible surfaces of chromia we use Miller notation to designate these surfaces

Miller Indices we define surfaces (planes) by the points they pass through along the lattice vectors • a surface plane with notation (l, m, n) will pass through the points a/l, b/m, c/n • in other words, figure out where the surface plane crosses each axis and take the reciprocal to get the Miller index • a Miller index of 0 means the surface plane never crosses the lattice vector • a bar over an index indicates that it is a negative number

- - (1120) (0112) Surfaces cleaving the bulk material leaves a slab surfaces are made by cleaving bulk structures: (0001)

Surfaces cleaving the bulk material leaves a slab surfaces are made by cleaving bulk structures: • system is periodic in all three dimensions (0001) • this is good for the two dimension defining the surface, since it yields an infinite slab vacuum • it’s bad for the direction perpendicular to the surface because the slab ends up repeated infinitely in that direction • to minimize interactions between slabs: slab use vacuum space (~10 Å, but test convergence of properties) ensure slab is symmetric to eliminate dipoles employ dipole corrections if dipoles are present vacuum

Surfaces cleaving the bulk material leaves a slab surfaces are made by cleaving bulk structures: (0001) • need to ensure that slab is sufficiently thick • basic idea is to use a slab that is thick enough for the middle to look bulk-like, while allowing the outer layers to fully relax vacuum • instead of using very thick slabs, you can also fix the positions of the middle atoms to their bulk values slab • ultimately, you should test the effects of changing the slab thickness vacuum

Surface Energies the surface energy is a key quantity obtained in calculations (0001) vacuum • M = number of stoichiometric units in slab model • N = number of stoichiometric units in bulk model • A = cross-sectional area of surface slab where a and b are the lattice vectors defining the plane of the surface vacuum

Surface Energies the surface energy is a key quantity obtained in calculations - Cr2O3 (0112) surface bulk Cr2O3 • 2 unit cells shown • 2 unit cells shown • 48 stoichiometric units / simluation cell • 6 stoichiometric units / simulation cell

Surface Energies the surface energy can be sensitive to the vacuum space and slab thickness • you should ensure convergence to less than 1 mJ/m2 • comparing surface energies allows you to determine which surface is most stable

Surface Energies the surface energy is also affected by relaxation of the slab • surface relaxation can have significant effects on surface energies • e.g. (0001) surface undergoes large stabilization with relaxation • due to large inward relaxation of terminal Cr3+ ions • reconstruction can also occur (e.g. dimerization on Si surfaces) relaxed (γr) and unrelaxed (γf) energies of different surfaces of Cr2O3

CHEM 938: Density Functional Theory