CHEM 938: Density Functional Theory

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

Examples of Condensed-Phase Systems many interesting processes occur in the condensed phase Crystalline: ordered network of atoms • examples are metals, ionic solids, covalent solids and molecular crystals bulk iron Amorphous: disordered network of atoms • examples are glass • could include liquids here amorphous SiO2

Examples of Condensed-Phase Systems many interesting processes occur in the condensed phase Combinations of solids and liquids: • films at interfaces, liquids moving through channels, etc. aldehydes within a sliding Al2O3 contact Surfaces: condensed system in 2-D (0001) - (0112) different surfaces of Cr2O3

Modeling Condensed Phases we represent condensed-phase systems using either of two models cluster models: • build up a small molecule (cluster) representing the system of interest • can be sufficient for surface calculations • not suited to bulk calculations (e.g. metallic systems rely on containing a large number of atoms and electronic states) • could contain point charges instead of atoms if a reactive site is known CrO33-: top Cr36O8766-: side [0001] CrO33-: side Cr36O8766-: top - [1010]

Modeling Condensed Phases we represent condensed-phase systems using either of two models cluster models: • build up a small molecule (cluster) representing the system of interest • can be sufficient for surface calculations • not suited to bulk calculations (e.g. metallic systems rely on containing a large number of atoms and electronic states) • could contain point charges instead of atoms if a reactive site is known quantum chemical part, remainder treated with point charges top view side view

Modeling Condensed Phases we represent condensed-phase systems using either of two models periodic models: • build up a small model representing the structure of interest, place it in a well-defined simulation cell, and repeat it infinitely in all three dimension • mirrors what occurs in solids • can be used to model liquids and amorphous solids, but some care needs to be taken because these systems are not periodic • surfaces and interfaces are often periodic in two dimension, so care needs to be taken in the treatment of the third • we are going to focus on periodic models when discussing condensed-phase systems

O Zn P Periodic Boundary Conditions we can represent a bulk system by replicating it infinitely in all three dimensions • define basic simulation cell • contains a certain number atoms and electrons • has cell vectors a, b, c • repeat cell infinitely in all three dimensions • done through mathematical means, we don’t actually treat a system with an infinite number of atoms • an atom at position r will also be at r + la + mb + nc, where l, m, and n are integers • wavefunction is also periodic • K is periodicity of crystal

Periodic Boundary Conditions (PBCs) care must be taken when setting up the simulation cell Crystalline materials: • should define cell on basis of known crystal structure • may want to consider using multiple unit cells because observable behavior is limited to size of cell • for example, unit cell vector lengths may only be a few Å, but relevant lengths of sound waves, phonons, stress fields, e.g. may be longer and thus don’t fit in the cell bcc unit cell bcc crystal with PBCs

O Zn P Periodic Boundary Conditions (PBCs) care must be taken when setting up the simulation cell Amorphous materials and liquids: • no basic repeating unit exists • usually use a large cell, so that the contents can be amorphous • repeating the ‘amorphous’ cell makes a crystal (no longer truly amorphous) • need to consider the relevant length scales for disorder amorphous simulation cell crystal of amorphous cells

Simulation Cells simulation cells can take on various shapes, but must be space-filling • cells are defined by three vectors: a, b, and c • can also be defined by lengths of vectors (a, b, c) and angles between vectors (α, β, γ) cubic crystal system a = b = c,  =  =  = 90° c    a b body-centered cubic face-centered cubic simple cubic

Simulation Cells simulation cells can take on various shapes, but must be space-filling • cells are defined by three vectors: a, b, and c • can also be defined by lengths of vectors (a, b, c) and angles between vectors (α, β, γ) hexagonal crystal system a = b  c,  =  = 90°,  = 120° a b c   c b  hexagonal a

Simulation Cells simulation cells can take on various shapes, but must be space-filling • cells are defined by three vectors: a, b, and c • can also be defined by lengths of vectors (a, b, c) and angles between vectors (α, β, γ) • there are other crystal systems based on 3-D parallelipeds • a good description can be found at: http://en.wikipedia.org/wiki/Crystal_structure

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • recall, that the wavefunction is also periodic • K is periodicity of crystal • atom-centered basis sets are not periodic • instead we use plane wave basis sets • these are inherently periodic

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • recall, that the wavefunction is also periodic • K is periodicity of crystal • consider a single sin function Ψ = sin(2πr/L) • exhibits correct periodic behavior (in fact it has to, which can be problematic at times) • but, it doesn’t look much like a good wavefunction around an atom • so, we represent our molecular orbitals as linear combinations of plane-wave basis functions

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • to get a wavefunction we have to take a Fourier series • where f(x) is a periodic function • in practice, this sum is truncated at some point, N (more on this later) example: Fourier series representation of a periodic square function with N = 1, 2, 3, and 4

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • sin functions are just the solutions of a particle-in-a-box Schrodinger equation • then we define the wavenumber, kn: • then energy is then expressed as: • so, we can relate kn to the momentum:

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • with sin (or cosine) functions we have the opportunity of working on real or momentum space • recall, in quantum mechanics any wavefunction can be written as a linear combination of solutions for any operator • choosing solutions of the momentum operator, we can write a wavefunction as: • the eigenfunctions of the momentum operator are:

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • with sin (or cosine) functions we have the opportunity of working on real or momentum space • the real-space wavefunction is then written as these are called Fourier transforms • the inverse approach can be applied:

Planewave Basis Sets we need to use basis sets that are easy to repeat periodically • with sin (or cosine) functions we have the opportunity of working on real or momentum space why would we want to work in either space? • many terms are easier to evaluate in momentum space: vs. • while other terms are easier to evaluate in real space • fast-Fourier transform (FFT) routines permit an easy switch between spaces • these scales as Kln(K) ~ linear scaling for large K

Planewave Basis Sets of course, we are dealing with 3-D structures • particle in a 3-D box • the wavenumber, k, becomes a wavevector, k: • and the Fourier transforms become:

Planewave Basis Sets plane waves can be generalized to arbitrary box shapes • define cell with vectors a, b, and c • construct wavevector, G • planewave basis functions are: • periodic functions can then be expressed as:

Planewave Basis Sets plane waves can be generalized to arbitrary box shapes • define cell with vectors a, b, and c • construct wavevector, G • Bloch theorem states that Kohn-Sham orbitals can be expressed as: • using planewaves to represent the periodic function u yields: we’ll get back to Gandk in a few slides

Planewave Basis Sets plane waves can be generalized to arbitrary box shapes • define cell with vectors a, b, and c • construct wavevector, G • planewave basis functions: • are periodic with respect to the simulation cell • they are orthonormal • treat all regions of space equally

Planewave Basis Sets how many planewaves do you use in a calculation? • the plane wave basis set is defined by a kinetic energy cutoff, Ecut • kinetic energy of plane wave G • so, we define a value Ecut and all values of G below |G|2/2 are included in the basis set • Ecut is usually given in units of Rydbergs or eV (to make things less accessible to chemists) • values of 20 to 150 Ry or 300 to 800 eV are common in the literature • higher values are needed for elements with 2p or 3d electrons in the outermost valence shells

Planewave Basis Sets so, how many plane waves (basis functions) is this?? • the number of planes waves is given by: • this depends on both Ecut and the volume of the simulation cell Example: • in a recent study of 2 CF3CF2OCF3 molecules in a periodic simulation cell, I used Ecut = 120 Ry • this translates into ~ 30,000 basis functions • by way of comparison, a calculation on the same system with the 6-31G(d,p) basis set would require 360 contracted Gaussian basis functions (672 primitive Gaussians) • but plane waves make certain parts of the calculation much faster than with atom-centered basis sets

Planewave Basis Sets Advantages of Plane Waves: 1. they are a natural choice for periodic systems 2. they make the mathematics very simple • using fast Fourier transforms, we can turn things like differentiation into multiplication, etc. 3. they treat all regions of space equally • focusing on regions around atoms introduces bias into calculations with atom-centered basis functions 4. improving basis set is straightforward • increase Ecut Disadvantages of Plane Waves: 1. we have to use a lot of basis functions • even with simplified mathematics, plane wave calculations are still computationally intensive 2. they treat all regions of space equally • we waste a lot of time optimizing the wavefunction in regions of space where there aren’t any electrons

Reciprocal Space using planewaves requires us to define a set of reciprocal space coordinates • we are representing our Kohn-Sham orbitals as: • we need to define coordinate systems for k • this will also be important when evaluating properties using Fourier transforms • the normal way of defining this space is by determining the set of axes that satisfy:

Reciprocal Space there are certain details that have to be considered to work in reciprocal (momentum) space • in real space, a point in the simulation cell is defined in terms of the lattice vectors a, b, and c • in reciprocal space, a point in the simulation cell is defined in terms of the reciprocal lattice vectors b1, b2, and b3 • vectors with large magnitude in real space become shorter in reciprocal space

Reciprocal Space there are certain details that have to be considered to work in reciprocal (momentum) space correct formulae if angle between a and b is 90° • consider a 2-D example reciprocal space real space

Reciprocal Space there are certain details that have to be considered to work in reciprocal (momentum) space correct formulae if angle between a and b is 90° • consider a 2-D example unit cell in reciprocal space: • draw lines midway along lattice vectors in all directions • area in enclosed by intersections is smallest repeating (primitive) unit cell in reciprocal space • called the Brillouin zone (BZ) • by performing calculations only in the BZ, one can obtain properties for a whole crystal reciprocal space

Reciprocal Space of course, we work in 3-D space fcc primitive cell Brillouin zone volume:

Properties in Reciprocal Space consider the evaluation of the electron density • ‘normally’: • with planewaves: all other real-space properties involve integration over the Brillouin zone

k-point Integration evaluating properties requires integration over the Brillouin zone • in practice, we do this numerically • to do this, we define a k-point grid • usually, this is done by specifying a number of divisions along each reciprocal lattice vector: k1 x k2 x k3 • each of these grid points are called k-points • to ensure integration is performed evenly, the spacing of the grid points should be consistent along all three directions (i.e. shorter reciprocal space lattice vectors get fewer divisions) • there are a few exceptions to the last point

k-point Integration evaluating properties requires integration over the Brillouin zone • in practice, we do this numerically • note that you have to evaluate the property (wavefunction, density, energy, etc.) at each k-point • thus obtaining the real-space wavefunction involves evaluating the wavefunction at k1xk2xk3 k-points • so, you want to ensure that you minimize the number of k-points used (ensuring that you retain good accuracy of course) • for large cells, a 1x1x1 grid is usually sufficient (metals are an exception)

k-point Integration some approaches take advantage of symmetry when setting k-point grids description symbol center of Brillouin zone Γ simple cube center of an edge M corner point R center of face X face-centered cubic middle of an edge joining two faces K center of a hexagonal face L middle of an edge joining a hexagonal and a square face U Brillouin zone of an fcc cell with high-symmetry points labeled corner point W center of square face X high symmetry k-points also play an important role in interpreting electronic structure

k-point Integration some approaches take advantage of symmetry when setting k-point grids description symbol center of Brillouin zone Γ body-centered cubic corner point joining four edges H center of a face N corner point joining three edges P hexagonal center of a hexagonal face A corner point H middle of an edge joining two rectangular faces K Brillouin zone of an fcc cell with high-symmetry points labeled middle of an edge joining a hexagonal and a rectangular face L high symmetry k-points also play an important role in interpreting electronic structure center of a rectangular face M

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 Effective Core Potentials: • 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 Effective Core Potentials: • 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

CHEM 938: Density Functional Theory