CHEM 938: Density Functional Theory

1. CHEM 938: Density Functional Theory molecular properties February 2, 2010

2. Potential Energy Surfaces structure-energy relationship is fundamental to chemistry 1. geometry optimization • determine structures and energies of reactants, products and transitions states • these structures are as stationary points: • gives reaction energetics: 2. molecular dynamics • follow spatiotemporal evolution of a chemical system acceleration positions velocity

3. Potential Energy Surfaces structure-energy relationship is fundamental to chemistry 3. statistical mechanical treatments • predict thermodynamic properties by averaging over molecular systems with different geometries 4. normal mode analysis in the harmonic approximation • gives vibrational frequencies, IR/Raman spectra • depend on curvature of potential energy surface

4. Potential Energy Surfaces structure-energy relationship is accessible within Born-Oppenheimer approximation • constant energy for a given set of nuclear positions {RI} • depend parametrically on nuclear positions {RI} • changing {RI} alters nuclear-electron attraction term

5. Potential Energy Surfaces structure-energy relationship is accessible within Born-Oppenheimer approximation in density functional theory: nuclear-nuclear repulsion energy nuclear-electron attraction operator molecular geometry atomic positions, {RI} nuclear charges, {ZI} nuclear-electron attraction operator defines Kohn-Sham operator for specific geometry

6. Potential Energy Surfaces structure-energy relationship is accessible within Born-Oppenheimer approximation in density functional theory: dependent on positions and charges of nuclei

7. 1-D Example: N2 • calculate energy at several distinct values of dN-N • called a potential energy surface scan • interpolate between points to construct PES • gives relevant structural and energetic information

8. 2-D Example: Ozone Isomerization ‡ dO1-O2 O1 O1 O1 O2 O2 O2 O2 O2 O2 aO2-O1-O2 ozone iso-ozone build up the PES: • systematically vary aO2-O1-O2 and dO1-O2 on a grid and calculate the energy at each set of values • requires a large number of calculations due to 2-D nature of system 2-D potential energy surface

9. 2-D Example: Ozone Isomerization ‡ dO1-O2 O1 O1 O1 O2 O2 O2 O2 O2 O2 aO2-O1-O2 ozone iso-ozone can identify: • reactant and product structures iso-ozone ozone • local minima on PES 2-D potential energy surface

10. 2-D Example: Ozone Isomerization ‡ dO1-O2 O1 O1 O1 O2 O2 O2 O2 O2 O2 aO2-O1-O2 ozone iso-ozone can identify: • reactant and product structures transition state • local minima on PES • transition state structure • saddle point connecting local minima 2-D potential energy surface

11. 2-D Example: Ozone Isomerization ‡ dO1-O2 O1 O1 O1 O2 O2 O2 O2 O2 O2 aO2-O1-O2 ozone iso-ozone can identify: • reactant and product structures reaction path • local minima on PES • transition state structure • saddle point connecting local minima • reaction pathway • lowest energy pathway connecting local minima 2-D potential energy surface

12. 2-D Example: Ozone Isomerization ‡ dO1-O2 O1 O1 O1 O2 O2 O2 O2 O2 O2 aO2-O1-O2 ozone iso-ozone reaction coordinate: • lowest energy path on full PES connecting reactant, product, and transition state • often plotted to aid in interpreting data useful for identifying: • reaction energies • reaction barriers • reactant, product and transition state geometries

13. What About Higher-Dimensional Cases? virtually all chemical systems have more than 2 degrees of freedom • in general, a molecule has 3N-6 degrees of freedom (3N-5 for linear molecules)  PES is a hypersurface Problem: we cannot plot a surface with more than 2 degrees of freedom Solution: treat hypersurface as mathematical entity • employ mathematical definitions to identify relevant points on the PES • better than manually scanning the PES • discreteness of points in a PES scan means we may miss minima • mathematical search for minima and transition states usually involves fewer calculations than a scan (particularly for large systems) So, what are we actually interested in doing? • want to identify minima and saddle-points (transition states) on the PES • mathematically, these are stationary points

14. Stationary Points reactants, products and transition states are stationary points • mathematically, a stationary point is defined by: along all directions, qi There are different kinds of stationary points • reactants and products are local minima: • transition states are saddle points: • minimum in all directions • maximum in only one direction, minima in all others • lowest energy minimum on the entire PES is called the global minimum

15. Stationary Points Local minimum (reactant or product): • for all qi • for all qi • stationary point • positive curvature in all directions Saddle Point (transition state): • for all qi • for exactly one qi • stationary point • negative curvature in one direction • maximum along that direction • for all other qi • positive curvature in all other directions • minimum along all other directions

16. Geometry Optimization • process for locating a minimum or transition state on the PES • also called relaxation (at least for finding minima) • most common type of calculation performed in computational chemistry Why? • need geometries/energies of minima and transition states to: • determine reaction mechanisms and associated energetics • determine preferred molecular geometries for calculations of other properties Basic idea: • iterative process • start with a ‘guess’ structure and use information of the PES to change the coordinates such that the final structure is a stationary point • lets us find minima/transition states without any exact information regarding those points on the PES

17. finds minimum in 1 step!! guess structure (qi, Ei) step Newton-Raphson Method • basis for geometry optimization procedures used in most computational chemistry codes consider a 1-D harmonic potential energy, E bond length, q minimum (q0, E0)

18. Newton-Raphson Method • in general, PES is not harmonic  cannot find q0 in 1 step • instead, we iterate until forces on nuclei (derivatives, gradients) are small • with a reasonable guess structure, most programs will find the nearest local minimum in 10 – 20 iterations • transition states can be more difficult to locate • NOTE: doesn’t just follow reaction coordinate, but optimizes over whole PES initial guess structure (qi) energy, E . . . local minimum 1 local minimum 2 reaction coordinate, q

19. Newton-Raphson Method in general, optimization is done using 3N cartesian coordinates

20. Newton-Raphson Method to do an optimization, we need: qi • column vector of 3N cartesian coordinates defining the molecular geometry • initial guess structure should look like anticipated minimum or transition state • can build with molecular editors gi • column vector of 3N derivatives of the energy wrt cartesian coordinates • represent the forces on the nuclei: • calculated by the program  ~20% of effort in geometry optimization Hi • “Hessian” matrix containing 3Nx3N 2nd derivatives of the energy wrt cartesian coordinates • 2nd derivatives are expensive to calculate analytically • generally estimated during optimization:

21. Forces on Nuclei Hellmann-Feynman theorem: do not depend directly on {RI} depend directly on {RI}

22. Forces on Nuclei Hellmann-Feynman theorem: do not depend directly on {RI} depend directly on {RI} yields derivatives as: • XA = x coordinate of atom A, analogous expressions exist for other coordinates

23. Forces on Nuclei Hellmann-Feynman theorem doesn’t necessarily apply to DFT calculations using atom-centered basis functions not true when basis functions depend on atomic positions Gaussian basis function depends on nuclear positions molecular orbitals

24. Forces on Nuclei Hellmann-Feynman theorem doesn’t necessarily apply to DFT calculations using atom-centered basis functions • called Pulay forces • in principle require derivative of wavefunction with respect to coordinates (expensive) • in practice, are evaluate in more efficient ways • not unique to DFT, but rather to using atom-centered basis functions

25. Forces on Nuclei Hellmann-Feynman theorem doesn’t necessarily apply to DFT calculations using atom-centered basis functions do not depend directly on {RI} depend directly on {RI} • grids used to evaluate exchange-correlation potential often depend on nuclear positions • electron-electron repulsion term in Hamiltonian also depends on coordinates

26. Forces on Nuclei Hellmann-Feynman theorem doesn’t necessarily apply to DFT calculations using atom-centered basis functions 0 • require derivatives of grids with respect to nuclear positions • some programs neglect these, while others evaluate them • can be responsible for different geometries with different programs, but differences should be slight assuming a sufficiently dense grid

27. Geometry Optimization – Input • keyword ‘opt’ requests a geometry optimization to a minimum energy structure • ‘hf/3-21G’ specifies level of theory for the calculation • ‘nosymm’ tells the program to set the initial symmetry to C1 • geometry optimizations sometimes change the symmetry of the system, causing the calculation to fail • the minimum energy structure may not have the same symmetry as the initial structure

28. Transition State Optimization – Input • keyword ‘opt’ requests a geometry optimization to a minimum energy structure • option ‘ts’ requests a transition state optimization • option ‘calcfc’ requests that the Hessian is calculated analytically at the first optimization step • option ‘noeigen’ requests that the Hessian is not tested throughout the calculation. If testing is permitted, the calculation often fails because at some steps the Hessian may not have the correct number of imaginary frequencies. • ‘hf/3-21G’ specifies level of theory for the calculation • ‘nosymm’ tells the program to set the initial symmetry to C1 • geometry optimizations sometimes change the symmetry of the system, causing the calculation to fail • the minimum energy structure may not have the same symmetry as the initial structure

29. Structure Optimization in Webmo transition state optimization options are correct by default

30. Geometry Optimization – Output in a geometry optimization we may want to determine: • a stationary point corresponding to a minimum energy structure  need to monitor whether the convergence criteria are met first step: intermediate step: final step:

31. Geometry Optimization – Output in a geometry optimization we may want to determine: • a stationary point corresponding to a minimum energy structure  need to monitor whether the convergence criteria are met • energy of the optimized structure  energy statement in step where convergence criteria are met  last energy statement in the output file

32. Geometry Optimization – Output in a geometry optimization we may want to determine: • a stationary point corresponding to a minimum energy structure  need to monitor whether the convergence criteria are met • energy of the optimized structure  energy statement in step where convergence criteria are met  last energy statement in the output file • geometry of the optimized structure  follows statement where convergence criteria are met

33. Geometry Optimization – Webmo Output

34. Performance DFT usually yields accurate geometries main group compounds: SVWN = local functional, others are GGA • all data within ~0.02 – 0.03 Å of experiment LDA and GGA functionals both provide relatively good accuracy for bond lengths

35. Performance DFT usually yields accurate geometries main group compounds: SVWN = local functional, others are GGA • all data within ~0.02 – 0.03 Å of experiment • usually, LDA bond lengths are too short and GGA bond lengths are too long LDA and GGA functionals both provide relatively good accuracy for bond lengths

36. Performance hybrid functionals often outperform GGAs and ab initio methods main group compounds:

37. Performance hybrid functionals often outperform GGAs and ab initio methods main group compounds: • Hartree-Fock bond lengths are usually too short • GGA bond lengths are usually too long • hybrid functionals combine Hartree-Fock and GGA • fortuitous cancellation of errors • yields better results

38. Performance DFT usually yields accurate geometries transition metal compounds: “Transition metal chemistry (…) is a graveyard for UHF-based MP methods” P. Taylor, 1992, “Accurate Calculations and Calibration” in Lecture Notes in Quantum Chemistry, B.O. Roos (ed). Springer, Heidelberg.

39. Performance DFT usually yields accurate geometries transition metal compounds: • BP86 often provides good results for TM complexes • recently supplanted by B3LYP, which seems to perform better • hybrid functionals combine Hartree-Fock and GGA

40. Performance DFT usually yields accurate geometries transition metal compounds: • BP86 often provides good results for TM complexes • recently supplanted by B3LYP, which seems to perform better note that other factors matter, too: • TM complexes require large basis sets • relativistic effects can be important (particularly further down the periodic table) • capturing orbital (quasi-)degeneracy, open-shell character, and static correlation effects so, there is probably some cancellation of error going on: “Always make two approximations in quantum chemistry – never one.” paraphrased from R. Car

41. Relative Energies once you identify the reactant, product, and transition state • energies can be calculated relative to reactant energy (kcal/mol) • sets reactant energy to 0 • all energies must be calculated at the same level of theory!!! reactants TS product For reactions with multiple reactants or products: e.g. A + B TS†  C + D • energies are calculated separately for each molecule • e.g. do a calculation on A and a calculation on B, not a calculation on A and B together

42. Performance reaction energies are usually pretty good, but dependent on the XC functional atomization energies are a major testing ground for energetics: CH4 C + 4H • very accurate energetics are available from high level ab initio methods • represents a significant challenge (almost worst case scenario) because reactant and products are very different • reactants: (usually) closed-shell with bonds between atoms → large correlation effects • products: open-shell systems with no bonds → lower correlation effects • because of difficulty associated with these calculations, the results can probably be taken as upper bounds on typical errors for reaction energies • the difficulty in these calculations due to large differences in the magnitude of correlation, once again reflects the inherent need for cancellation of errors in DFT (and other quantum chemical) methods

43. Performance reaction energies are usually pretty good, but dependent on the XC functional atomization energies are a major testing ground for energetics: • JGP test set contains 32 first and second row species • DFT clearly outperforms HF, MP2 and QCISD in this case • poor performance of MP2 and QCISD can be attributed partially to use of a small basis set • functionals with GGA exchange (BVWN and BLYP) provide best results • ab initio methods underbind, LDA overbinds, and GGA is just about right

44. Performance reaction energies are usually pretty good, but dependent on the XC functional atomization energies are a major testing ground for energetics: • SVWN = LDA exchange, LDA correlation, SLYP = LDA exchange, GGA correlation, BVWN = GGA exchange, LDA correlation, BLYP = GGA exchange, GGA correlation • using (at least) GGA exchange is necessary to reliably obtain good energetics • exchange is the largest term in the XC energy (approx. 10x greater than correlation)

45. Performance hybrid functionals often outperform their GGA counterparts atomization energies are a major testing ground for energetics:

46. Performance hybrid functionals often outperform their GGA counterparts atomization energies are a major testing ground for energetics: mean absolute errors = no brackets maximum absolute errors = brackets

47. Performance hybrid functionals often outperform their GGA counterparts atomization energies are a major testing ground for energetics: • SVWN = LDA; PBE = GGA, non-empirical; BLYP = GGA, empirical; B3LYP = hybrid, empirical; PBE1PBE = hybrid, empirical; VSXC = meta-GGA, empirical in general, expected accuracy: SVWN << BP86 < BLYP ≈ BLYP < B3P86 < B3LYP ≈ B3PW91

48. Performance GGA and hybrid functionals generally perform well for transition metals TM complexes:

49. Performance GGA and hybrid functionals generally perform well for transition metals TM complexes:

50. Performance GGA and hybrid functionals generally perform well for transition metals TM complexes: M. Quintal et al., J. Phys. Chem. A, 110, 709 (2007)