**introduction to Ab Initio Molecular Dynamics**

introduction to Ab Initio Molecular Dynamics

(Parte **1** de 3)

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 1

References Material for these lecture notes was compiled from the references below

9MIT’s 3.320 course notes (Prof. N. Marzari)

9R. M. Martin, Electronic Structure: Basic Theory and Methods, Cambridge University Press, Cambridge, U.K., 2004.

9W. Koch and M.C. Holthausen, A Chemist's Guide to Density Functional Theory, 2nd Edition

Introduction to Pseudopotentials, applying DFT to solids

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 2

Electronic structure calculations

¾In this lecture, we will review the basics of DFT and discuss practical aspects in the implementation of DFT. In particular, we will discuss:

9How do you set up a DFT calculation? (how do you obtain from the variational principle a linear algebra problem?) 9 Pseudopotentials 9Bloch theorem for an electron in an infinite solid 9Using plane waves as a basis set 9Discretization in electronic structure calculations –accuracy versus computational effort. ¾Examples of calculation of material properties using DFT.

¾Pitfalls of DFT calculations.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 3

Review of the fundamentals of DFT

Schrödinger equation to find the ground state charge density, the ground state energy and other ground state properties.

¾We can compute all these using a variational principle on the charge density. The functional that we need to minimize is the following:

¾We need to vary the charge density n’ and the minimum value that this functional will take is the ground state energy E 0 .

¾We also showed that eventually we have a Schrödinger like equation expressed in terms of the charge density alone (2nd Hohenberg-Kohn theorem).

¾However, note that while everything is well defined in principle it is not clear how it works in practice.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 4

The universal functional

¾The Hohenberg and Kohn 1st theoremstates that given any charge density n’, an external potential v’ext is well defined for which that charge density is going to be the ground state charge density:

n’ Æv’ext ÆΨ’ Æn’

¾ The solution to the Schrödinger equation corresponding to v’ ext is well defined and is denoted above as Ψ’.

¾One can construct a functional Fthat is just the expectation value of the many body kinetic energy operator plus the many body electron-electron interaction for the state Ψ’ .

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 5

¾In the variational principle, in addition to | we included an additional term |

Review of the fundamentals of DFT that is the integral of the external potential times the charge density.

¾This variational principle is a reformulation of the many body Schrödinger equation. To make this a practical method, Kohn and Sham proposed an approximation to the universal functional.

¾In this universal functional, | , the first term that we can extract is the |

Hartree electrostatic energythat is a simple and well defined functional of the charge density.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 6

The universal functional

¾It’s not easy to extract the kinetic energy (2nd derivative of the wave function) from the charge density.We already discussed this difficulty for a plane wave.

A plane wave of any wavelength gives you a constant charge density but gives you very different 2 nd derivative (kinetic energy).

¾Kohn and Sham introduced a related system of non-interacting electrons that they have the same ground state charge density as our original problem.Sincethese are non-interacting electrons, one can define exactly their quantum kinetic energy.

¾This kinetic energy term | together with the Hartree term |

makes our third term in the universal functional E xc [n](exchange correlation term) to be very small. Theexchange correlation term contains all the complexity of the problem. It is well defined in principle but we do not know how to compute.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 7

The interacting homogeneous gas approximation

¾Kohn and Sham then applied the Thomas and Fermi idea. They used the the exchange correlation energy density of the interacting homogeneous electron gas at a given density. Then they computed the exchange correlation energy term in the universal functional by integrating assuming in each infinitesimal volume that the exchange correlation energy density is that of the interacting homogeneous electron gas at that density.

¾Ceperley and Alder (1980) did the first quantum Monte Carlo calculations of the interacting homogeneous electron gas as a function of the density. This functional was for the first time parameterized over a whole set of densities (Perdew-Zunger parameterization ).

¾Our computational problem now becomes the calculation of the minimum of the energy functional:

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 8

Exchange correlation in Hartree-Fock vs DFT

¾Recall that inHartree-Fock we also had a variational principle but the wave function approximation (Slater determinant) led us to a well defined exchange term.

¾Since the exchange term in Hartree-Fock was coming from a variational principle we had the possibility of systematically making the Hartree-Fock approach better and better by extending the class of wavefunctions from one single Slater determinant to e.g. combinations of Slater determinants.

¾Density functional theory does not give us any systematic way to improve our solution. It leads to a decomposition of terms with a well defined (in principle) exchange correlation functional. It is for this reason that took several years for DFT to become a practical electronic structure calculation method.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 9

The first applications of DFT

¾Only after the LDA approximation of Ceperley and Alder practical applications of

DFT started to appear. We have seen in the introductory lecture to DFT how the phase diagram of Siwas computed, also lattice parameters, the bulk modulus, etc. Here we plot the energy as a function of the lattice parameter. We also discussed phase transition to high pressure phases of Si.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 10

Total Energy in DFT

¾Once the exchange correlation potential is written in an explicit form, the total energy is a well defined non-linear functional of the independent single particle orbitals (the Kohn-Sham orbitals) of which the charge density is just the sum of the square moduli.

¾We have to address the (non-linear) minimization of the total energyin the DFT theory. There are two distinct routes that we can take.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 1

Minimization problem and the Kohn-Sham equations

¾We need to find the orbitals that minimize this expression or, as it often happens, we can write the associated Euler-Lagrange equation.

¾We can take the functional differential of the functional and that gives us the Kohn-Sham equations .

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 12

Kohn-Sham Equations

¾The Kohn-Sham equations are similar to the Hartree-Fock equations. We still have a quantum kinetic energy term, a Hartree term, an external potential term, and

Ψi are our single particle orbitals.

¾As discussed earlier,the main difference from the Hartree-Fock equations is in the way the exchange correlation term is calculated. In Hartree-Fock, the exchange correlation term was an explicit integral of the orbitals. In DFT is a complex function of the charge density:

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 13

Kohn-Sham Equations

¾The Hamiltonian acting on this single particle orbitals depend on the charge density (given by the sum of the square moduli) and the charge density is a function of the orbitals themselves.

¾As expected the Hamiltonian is self-consistent(i.e. the operator depends on its own solution).

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 14

LDA approximation

¾Asin Hartree-Fock you have found your ground state solution only once you have become self-consistent.You have a Hamiltonian with eigenstates that give you a charge density that gives you the same Hamiltonian you have calculated the eigenstates for!

¾The computational approaches to solve the DFT problem, as those that solve the Hartree-Fock problem, are iterative approaches.

¾The quality of the DFT calculations depends ultimately on the quality of your functional.For many years the only functional that was used was the local density approximation (LDA) functional.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 15

GGA –Generalized gradient approximations

¾Functionals of the charge density have also been constructed that didn’t only depend on the charge density but also on its gradient. As discussed in the earlier lecture, this was done not in the sense of a Taylor series expansionbecause a Taylor expansion would not satisfy symmetry properties that we know that the exact exchange correlation functional would do. These approximations are called generalized gradient approximations.

¾Many of the DFT approximations are denoted with the first letter of the people who invented them. For example, P (for John Perdew), or B (for Axel Becke), or B (for Tyrone Burke). Visit this web site for review of many exchange correlation functional approximations.

¾The most used method today beyond LDA is PBE (J.P. Perdew, K. Burke, and

M. Ernzerhof, Phys. Rev. Lett. 7, 3865, 1996 ).

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 16

Hybrid exchange correlation functionals

¾The quantum chemistry community has developed hybrid functionalswith a certain percentage of the exchange correlation functional coming from the DFT approximation (e.g. PBE) and some percentage of Hartree-Fock exchange mixed in.

¾B3L YP and PBE0work very well especially for molecules. Hartree-Fock comes from atoms and molecules and it tends to work better in that limit. DFT is based from the homogeneous electron gas (LDA approximation) and tends to work better for solids.

¾When you want to study a molecule on a solid surface (for adsorption of H on metallic surfaces), then none of these two approaches (Hartree-Fock or DFT) works really very well alone and mixing of the two methods is needed.

¾There is a lot of work in developing more complex exchange correlation functionals that tend to work better than any of the ones mentioned earlier. These functionals in addition to the density and its gradient may also depend on the Laplacian or even the orbital themselves. There are not used in practical applications due to their computational complexity.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 17

Results from DFT: Lattice parameter

¾Using LDA or GGA we can accurately describe a number of properties for very different materials. The table above from (C.J. Pickard, 2002) compares for metals, semiconductors, oxides, alloys, the experimental lattice parameter with the DFT prediction.

Chris J. Pickard

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 18

Results from DFT: Lattice parameter

¾Note that the error is ~1%. The error on a lattice parameter of a material that does not have exotic electronic properties (e.g. not a strongly correlated electronic material) can be expected to be in the range of an error ~1-2%.

¾Electronic structure calculations tend to be extremely powerful because they are not fitted (they do not fit any potential).

Chris J. Pickard

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 19

¾We will discuss here the total energy pseudopotential approach (read paper by

Joannopoulos et al .) that was developed to study solids.

¾Describing accurately the core electrons in an atom or in a solid is exceedingly complex .

¾The core electrons are unaffected by the chemical environment. For example, if we take an Fe atom and you put it in Fe as a metal, or you put it as Fe in a transition metal oxide, the valence electrons of Fe will redistribute themselves.

However, the core electrons of Fe are so tightly bound to the nucleus by order of magnitude in energy with respect to the typical energy of valence electrons, that they are basically unaffected.

¾We do not want to carry the computational expense of describing the core electronswhen we know that their rigid contribution does not change depending on the chemical environment.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 20

¾The core electrons are bound to the nucleus and they screen the nucleus charge. In Fe, the two 1s electrons will be so tightly bound to the nucleus that the nucleus does not look to the 2s or 2p or 3s, 3p electrons as having 26 protons but it really looks like having 24 protons because the two 1s electrons screen completely.Similarly, the 2s and the 2p will also screen almost completely the nucleus from the point of view of the 3s and the 3p electrons, etc.

¾We need to account for the presence of the core electrons but we don’t want to carry them in our electronic structure calculation as that will be too expensive. Recall from earlier lectures that the spatial variation of core electrons is very sharp and we will thus need a lot of computational information to describe all the sharp wiggles of the core electrons around the nucleus.

¾This problem was addressed using the so called pseudopotential approaches.

¾Once we have removed the core electrons from our problem, we need to find out what are the Kohn-Sham orbitals that minimize the DFT functional and we need to find an appropriate computational representation for these orbitals.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 21

DFT in practice

¾We will expand the Kohn-Sham orbitals in a basis. This basis set needs to be flexible enough to describe all the possible wiggles of our orbitals but also needs to be computationally efficient.

¾Once we have our external potential being represented by the nuclei and by the pseudopotential approach accounting for the core electrons, and once we have selected a basis set, one can proceed in solving self-consistently either the Kohn- Sham equations or minimize the non-linear energy functional.

¾Since the Hamiltonian depends on the charge density, we need to have a guess for the initial charge density to construct our operator. We need to start with a trial charge density and trial orbitals.

¾We can calculate the quantum kinetic energy, the Hartree energy, the exchange correlation terms and solve our problem.

¾From the Hamiltonian we’l find new orbitals. With these new orbitals we’l calculate the ground state charge density and we’l obtain a new Hamiltonian that then will keep iterating finding new orbitals, new charge density, new Hamiltonians until we reach self-consistency.

¾Alternatively, we can minimize the total energy functional to self-consistency.

MAE 715 –Atomistic Modeling of Materials

N. Zabaras (2/25/2009) 2

Summary of the DFT approach

¾Use the pseudopotential method to remove the core electrons ¾Represent the Kohn-Sham orbitalsas plane waves

¾Calculate the total energy for trial orbitals

¾Compute the kinetic energy and Hartree energy in reciprocal space ¾Compute exchange-correlation and external potential in real space

(Parte **1** de 3)