Ch121a Atomic Level Simulations of Materials and Molecules

Ch121a Atomic Level Simulations of Materials and Molecules BI 115 Hours: 2:30-3:30 Monday and Wednesday Lecture or Lab: Friday 2-3pm (+3-4pm) Lecture 2, March 30, 2011 QM-2: DFT William A. Goddard III, wag@wag.caltech.edu Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology Teaching Assistants Wei-Guang Liu, Fan Lu, Jose Mendozq

Last Time

Energy for 2-electron product wavefunction Consider the product wavefunction Ψ(1,2) = ψa(1) ψb(2) And the Hamiltonian for H2 molecule H(1,2) = h(1) + h(2) +1/r12 + 1/R In the details slides next, we derive E = < Ψ(1,2)| H(1,2)|Ψ(1,2)>/ <Ψ(1,2)|Ψ(1,2)> E = haa + hbb + Jab + 1/R where haa =<a|h|a>, hbb =<b|h|b> are just the 1 electron energies Jab ≡ <ψa(1)ψb(2) |1/r12 |ψa(1)ψb(2)>=ʃ d3r1[ψa(1)]2 ʃd3r2[ψb(2)]2/r12 = = ʃ [ψa(1)]2 Jb (1) = <ψa(1)| Jb (1)|ψa(1)> Where Jb (1) = ʃ [ψb(2)]2/r12 is the Coulomb potential at 1 due to the density distribution [ψb(2)]2 Jab is the Coulomb repulsion between densities ra=[ψa(1)]2 and rb

The general case of 2M electrons For the general case the HF closed shell wavefunction Ψ(1,2….2M) = A[(φ1a)(φ1b)(φ2a)(φ2b)… )(φMa)(φMb)] leads to HHFφk = lkφk where we solve for k=1,M occupied orbitals HHF = h + Σj=1,M [2Jj-Kj] This is the same as the Hamiltonian for a one electron system moving in the average electrostatic and exchange potential, 2Jj-Kj due to the other N-1 = 2M-1 electrons Problem: sum over 2Jj leads to 2M Coulomb terms, not 2M-1 This is because we added the self Coulomb and exchange terms But (2Jk-Kk) φk = (Jk) φk so that these self terms cancel. The HF equations describe each electron moving in the average potential due to all the other electrons.

The Matrix HF equations The HF equations are actually quite complicated because Kj is an integral operator, Kjφk(1) = φj(1)ʃ d3r2 [φj(2) φk(2)/r12] The practical solution involves expanding the orbitals in terms of a basis set consisting of atomic-like orbitals, φk(1) = Σm CmXm, where the basis functions, {Xm, m=1, MBF} are chosen as atomic like functions on the various centers As a result the HF equations HHFφk = lkφk Reduce to a set of Matrix equations ΣjmHjmCmk = ΣjmSjmCmklk This is still complicated since the Hjm operator includes exchange terms

New stuff

HF wavefunctions Good distances, geometries, vibrational levels But breaking bonds is described extremely poorly energies of virtual orbitals not good description of excitation energies cost scales as 4th power of the size of the system.

Electron correlation In fact when the electrons are close (rij small), the electrons correlate their motions to avoid a large electrostatic repulsion, 1/rij Thus the error in the HF equation is called electron correlation For He atom E = - 2.8477 h0 assuming a hydrogenic orbital exp(-zr)E = -2.86xx h0 exact HF (TA look up the energy) E = -2.9037 h0 exact Thus the elecgtron correlation energy for He atom is 0.04xx h0 = 1.x eV = 24.x kcal/mol. Thus HF accounts for 98.6% of the total energy

Configuration interaction Consider a set of N-electron wavefunctions: {i; i=1,2, ..M} where < i|j> = dij {=1 if i=j and 0 if i ≠j) Write approx = S (i=1 to M) Ci i Then E = < approx|H|approx>/< approx|approx> E= < Si Ci i |H| Si Cj j >/ < Si Ci i | Si Cj j > How choose optimum Ci? Require dE=0 for all dCi get Sj <i |H| Cj j > - Ei< i | Cj j > = 0 ,which we write as HCi = SCiEi in matrix notation, ie ΣjkHjkCki = ΣjkSjkCkiEi where Hjk = <j|H|k > and Sjk = < j|k > and Ci is a column vector for the ith eigenstate

Configuration interaction upper bound theorm Consider the M solutions of the CI equations HCi = SCiEi ordered as i=1 lowest to i=M highest Then the exact ground state energy of the system Satisfies Eexact ≤ E1 Also the exact first excited state of the system satisfies E1st excited ≤ E2 etc This is called the Hylleraas-Unheim-McDonald Theorem

Alternative to Hartree-Fork, Density Functional Theory Walter Kohn’s dream: replace the 3N electronic degrees of freedom needed to define the N-electron wavefunction Ψ(1,2,…N) with just the 3 degrees of freedom for the electron density r(x,y,z). It is not obvious that this would be possible but P. Hohenberg and W. Kohn Phys. Rev. B76, 6062 (1964). Showed that there exists some functional of the density that gives the exact energy of the system Kohn did not specify the nature or form of this functional, but research over the last 46 years has provided increasingly accurate approximations to it. Walter Kohn (1923-) Nobel Prize Chemistry 1998

The Hohenberg-Kohn theorem The Hohenberg-Kohn theorem states that if N interacting electrons move in an external potential, Vext(1..N), the ground-state electron density r(xyz)=r(r) minimizes the functional E[r] = F[r] + ʃ r(r) Vext(r) d3r where F[r] is a universal functional of r and the minimum value of the functional, E, is E0, the exact ground-state electronic energy. Here we take Vext(1..N) = Si=1,..NSA=1..Z [-ZA/rAi], which is the electron-nuclear attraction part of our Hamiltonian. HK do NOT tell us what the form of this universal functional, only of its existence

Proof of the Hohenberg-Kohn theorem Mel Levy provided a particularly simple proof of Hohenberg-Kohn theorem {M. Levy, Proc. Nat. Acad. Sci.76, 6062 (1979)}. Define the functional O as O[r(r)] = min <Ψ|O|Ψ> |Ψ>r(r) where we consider all wavefunctions Ψ that lead to the same density, r(r), and select the one leading to the lowest expectation value for <Ψ|O|Ψ>. F[r] is defined as F[r(r)] = min <Ψ|F|Ψ> |Ψ>r(r) where F = Si [- ½ i2] + ½ Si≠k [1/rik]. Thus the usual Hamiltonian is H = F + Vext Now consider a trial function Ψapp that leads to the density r(r) and which minimizes <Ψ|F|Ψ> Then E[r] = F[r] + ʃ r(r) Vext(r) d3r = <Ψ|F +Vext|Ψ> = <Ψ|H|Ψ> Thus E[r] ≥ E0 the exact ground state energy.

The Kohn-Sham equations Walter Kohn and Lou J. Sham. Phys. Rev.140, A1133 (1965). Provided a practical methodology to calculate DFT wavefunctions They partitioned the functional E[r] into parts E[r] = KE0 + ½ ʃʃd3r1 d3r2 [r(1) r(2)/r12 + ʃd3rr(r) Vext(r) + Exc[r(r)] Where KE0 = Si <φi| [- ½ i2 | φi> is the KE of a non-interacting electron gas having density r(r). This is NOT the KE of the real system. The 2nd term is the total electrostatic energy for the density r(r). Note that this includes the self interaction of an electron with itself. The 3rd term is the total electron-nuclear attraction term The 4th term contains all the unknown aspects of the Density Functional

Solving the Kohn-Sham equations Requiring that ʃ d3r r(r) = N the total number of electrons and applying the variational principle leads to [d/dr(r)] [E[r] – m ʃ d3r r(r) ] = 0 where the Lagrange multiplier m = dE[r]/dr = the chemical potential Here the notation [d/dr(r)] means a functional derivative inside the integral. To calculate the ground state wavefunction we solve HKSφi = [- ½ i2 + Veff(r)] φi = eiφi self consistently with r(r) = S i=1,N <φi|φi> where Veff (r) = Vext (r) + Jr(r) + Vxc(r) and Vxc(r) = dEXC[r]/dr Thus HKS looks quite analogous to HHF

The Local Density Approximation (LDA) We approximate Exc[r(r)] as ExcLDA[r(r)] = ʃ d3r eXC(r(r)) r(r) where eXC(r(r)) is derived from Quantum Monte Carlo calculations for the uniform electron gas {DM Ceperley and BJ Alder, Phys.Rev.Lett.45, 566 (1980)} It is argued that LDA is accurate for simple metals and simple semiconductors, where it generally gives good lattice parameters It is clearly very poor for molecular complexes (dominated by London attraction), and hydrogen bonding

Generalized gradient approximations The errors in LDA derive from the assumption that the density varies very slowly with distance. This is clearly very bad near the nuclei and the error will depend on the interatomic distances As the basis of improving over LDA a powerful approach has been to consider the scaled Hamiltonian

LDA exchange Here we say that in LDA each electron interacts with all N electrons but should be N-1. The exchange term cancels this extra term. If density is uniform then error is proportional to 1/N. since electron density is r = N/V

Generalized gradient approximations F(s) GGA functionals Becke 88 X3LYP PBE PW91 s

adiabatic connection formalism The adiabatic connection formalism provides a rigorous way to define Exc. It assumes an adiabatic path between the fictitious non-interacting KS system (λ = 0) and the physical system (λ = 1) while holding the electron density r fixed at its physical λ = 1 value for all λ of a family of partially interacting N-electron systems: is the exchange-correlation energy at intermediate coupling strength λ. The only problem is that the exact integrand is unknown. Becke, A.D. J. Chem. Phys. (1993), 98, 5648-5652. Langreth, D.C. and Perdew, J. P. Phys. Rev. (1977), B 15, 2884-2902. Gunnarsson, O. and Lundqvist, B. Phys. Rev. (1976), B 13, 4274-4298. Kurth, S. and Perdew, J. P. Phys. Rev. (1999), B 59, 10461-10468. Becke, A.D. J. Chem. Phys. (1993), 98, 1372-1377. Perdew, J.P. Ernzerhof, M. and Burke, K. J. Chem. Phys. (1996), 105, 9982-9985. Mori-Sanchez, P., Cohen, A.J. and Yang, W.T. J. Chem. Phys. (2006), 124, 091102-1-4.

Becke half and half functional assume a linear model take the exact exchange of the KS orbitals approximate partition set Get half-and-half functional Becke, A.D. J. Chem. Phys. (1993), 98, 1372-1377

Becke 3 parameter functional The success of B3LYP in achieving high accuracy demonstrates that errors of for covalent bonding arise principally from the λ  0 or exchange limit, making it important to introduce some portion of exact exchange Empirically modify half-and-half Becke, A.D. J. Chem. Phys. (1993), 98, 5648-5652. Becke, A.D. J. Chem. Phys. (1993), 98, 1372-1377. Perdew, J.P. Ernzerhof, M. and Burke, K. J. Chem. Phys. (1996), 105, 9982-9985. Mori-Sanchez, P., Cohen, A.J. and Yang, W.T. J. Chem. Phys. (2006), 124, 091102-1-4.

Some popular DFT functionals LDA: Slater exchange Vosko-Wilk-Nusair correlation, etc GGA: Exchange: B88, PW91, PBE, OPTX, HCTH, etc Correlations: LYP, P86, PW91, PBE, HCTH, etc Hybrid GGA:B3LYP, B3PW91, B3P86, PBE0, B97-1, B97-2, B98, O3LYP, etc Meta-GGA: VSXC, PKZB, TPSS, etc Hybrid meta-GGA:tHCTHh, TPSSh, BMK, etc

Truhlar’s DFT functionals MPW3LYP, X1B95, MPW1B95, PW6B95, TPSS1KCIS, PBE1KCIS, MPW1KCIS, BB1K, MPW1K, XB1K, MPWB1K, PWB6K, MPWKCIS1K MPWLYP1w,PBE1w,PBELYP1w, TPSSLYP1w G96HLYP, MPWLYP1M , MOHLYP M05, M05-2x M06, M06-2x, M06-l, M06-HF Hybrid meta-GGA M06 = HF + tPBE + VSXC

XYG3 approach to include London Dispersion in DFTGörling-Levy coupling-constant perturbation expansion where is the electron-electron repulsion operator, is the local exchange operator, and is the Fock-like, non-local exchange operator. where Combine both approaches (2 choices for b) a double hybrid DFT that mixes some exact exchange into while also introducing a certain portion of into contains the double-excitation parts of Take initial slope as the 2nd order correlation energy: Sum over virtual orbtials Substitute into with or This is a fifth-rung functional (R5) using information from both occupied and virtual KS orbitals. In principle can now describe dispersion

Final form of XYG3 DFT we adopt the LYP correlation functional but constrain c4 = (1 – c3) to exclude compensation from the LDA correlation term. This constraint is not necessary, but it eliminates one fitting parameter. Determine the final three parameters {c1, c2, c3} empirically by fitting only to the thermochemical experimental data in the G3/99 set of 223 molecules: Get {c1 = 0.8033, c2 = 0.2107, c3 = 0.3211} and c4 = (1 – c3) = 0.6789 Use 6-311+G(3df,2p) basis set XYG3 leads to mean absolute deviation (MAD) =1.81 kcal/mol, B3LYP: MAD = 4.74 kcal/mol. M06: MAD = 4.17 kcal/mol M06-2x: MAD = 2.93 kcal/mol M06-L: MAD = 5.82 kcal/mol . G3 ab initio (with one empirical parameter): MAD = 1.05 G2 ab initio (with one empirical parameter): MAD = 1.88 kcal/mol but G2 and G3 involve far higher computational cost.

Thermochemical accuracy with size G3/99 set has 223 molecules: G2-1: 56 molecules having up to 3 heavy atoms, G2-2: 92 additional molecules up to 6 heavy atoms G3-3: 75 additional molecules up to 10 heavy atoms. B3LYP: MAD = 2.12 kcal/mol (G2-1), 3.69 (G2-2), and 8.97 (G3-3) leads to errors that increase dramatically with size B2PLYP MAD = 1.85 kcal/mol (G2-1), 3.70 (G2-2) and 7.83 (G3-3) does not improve over B3LYP M06-L MAD = 3.76 kcal/mol (G2-1), 5.71 (G2-2) and 7.50 (G3-3). M06-2x MAD = 1.89 kcal/mol (G2-1), 3.22 (G2-2), and 3.36 (G3-3). XYG3, MAD = 1.52 kcal/mol (G2-1), 1.79 (G2-2), and 2.06 (G3-3), leading to the best description for larger molecules.

Accuracy (kcal/mol) of various QM methods for predicting standard enthalpies of formation

Comparison of QM methods for reaction surface ofH + CH4 H2 + CH3 H + CH4 H2 + CH3 HF CCSD(T) XYG3 Energy (kcal/mol) SVWN HF_PT2 B3LYP BLYP SVWN Reaction Coordinate: R(CH)-R(HH) (in Å)

Reaction barrier heights Note: no reaction barrier heights used in fitting the 3 parameters in XYG3) Zhao and Truhlarcompiled benchmarks of accurate barrier heights in 2004 includes forward and reverse barrier heights for 19 hydrogen transfer (HT) reactions, 6 heavy-atom transfer (HAT) reactions, 8 nucleophilic substitution (NS) reactions and 5 unimolecular and association (UM) reactions.

A. Total Energy (kcal/mol) B. Exchange Energy (kcal/mol) B B3LYP BLYP S HF XYG3 Distance (A) VWN B3LYP B3LYP LYP CCSD(T) LDA (SVWN) XYG3 HF_PT2 CCSD(T) PT2 XYG3 C. Correlation Energy (kcal/mol) Distance (A) Test for London Dispersion Conclusion: XYG3 provides excellent accuracy for London dispersion, as good as CCSD(T)

Accuracy of QM methods for noncovalent interactions. Note: no noncovalent complexes used in fitting the 3 parameters in XYG3) HB: 6 hydrogen bond complexes, CT 7 charge-transfer complexes DI: 6 dipole interaction complexes, WI:7 weak interaction complexes, PPS: 5 p-p stacking complexes. WI and PPS dominated by London dispersion.

Problem where is the electron-electron repulsion operator, is the local exchange operator, and is the Fock-like, non-local exchange operator. where XYG3 approach to include London Dispersion in DFTGörling-Levy coupling-constant perturbation expansion Take initial slope as the 2nd order correlation energy: Sum over virtual orbtials EGL2 involves double excitations to virtuals, scales as N5 with size MP2 has same critical step Yousung Jung (KAIST) has figured out how to get linear scaling for MP2 XYGJ-OS and XYGJ-OS

XYGJ-OS method most accurate DFT(including London Dispersion) at modest cost Yousung Jung a double hybrid DFT that mixes some exact exchange into while also introducing a certain portion of into contains the double-excitation parts of Ying Zhang , Xin Xu, Goddard; P. Natl. Acad. Sci. 106 (13) 4963-4968 (2009) Use Görling-Levy coupling-constant perturbation expansion c4 = (1 – c3)  3 parameters XYG3 most accurate DFT, but costs too high for large systems Yousung Jung figured out how to dramatically reduce the costs while retaining the accuracy XYGJ-OS

Timings XYGJ-OS and XYGJ-LOS for long alkanes XYGJ-LOS XYGJ-OS XYGJ-OS XYGJ-LOS

Accuracy of Methods (Mean absolute deviations MAD, in eV) HOF = heat of formation; IP = ionization potential, EA = electron affinity, PA = proton affinity, BDE = bond dissociation energy, NHTBH, HTBH = barrier heights for reactions, NCIE = the binding in molecular clusters

Comparison of speeds

Heats of formation (kcal/mol) small molecules Large molecules

Truhlar NHTBH38/04 set and HTBH38/04 set Reaction barrier heights (kcal/mol)

Truhlar NCIE31/05 set Nonbonded interaction (kcal/mol)

Comparison of QM methods for reaction surface of H + CH4 H2 + CH3 H + CH4 H2 + CH3 HF CCSD(T) XYG3 Energy (kcal/mol) SVWN HF_PT2 B3LYP BLYP SVWN Reaction Coordinate: R(CH)-R(HH) (in Å)

DFT-ℓg for accurate Dispersive Interactions for Full Periodic Table Hyungjun Kim, Jeong-Mo Choi, William A. Goddard, III 1Materials and Process Simulation Center, Caltech 2Center for Materials Simulations and Design, KAIST

Current challenge in DFT calculation for energetic materials C6 single parameter from QM-CC d =1 Reik = Rei + Rek (UFF vdW radii) Current implementations of DFT describe well strongly bound geometries and energies, but fail to describe the long range van der Waals (vdW) interactions. Get volumes ~ 10% too large XYGJ-lOS solves this problem but much slower than standard methods DFT-low gradient (DFT-lg) model accurate description of the long-range1/R6 attraction of the London dispersion but at same cost as standard DFT

PBE-lg for benzene dimer T-shaped Sandwich Parallel-displaced PBE-lg parameters Clg-CC=586.8, Clg-HH=31.14, Clg-HH=8.691 RC = 1.925 (UFF), RH = 1.44 (UFF) First-Principles-Based Dispersion Augmented Density Functional Theory: From Molecules to Crystals’ Yi Liu and wag; J. Phys. Chem. Lett., 2010, 1 (17), pp 2550–2555

DFT-lg description for benzene • PBE-lg predicted the EOS of benzene crystal (orthorhombic phase I) in good agreement with corrected experimental EOS at 0 K (dashed line). • Pressure at zero K geometry: PBE: 1.43 Gpa; PBE-lg: 0.11 Gpa • Zero pressure volume change: PBE: 35.0%; PBE-lg: 2.8% • Heat of sublimation at 0 K: Exp:11.295 kcal/mol; PBE: 0.913; PBE-lg: 6.762

DFT-lg description for graphite Binding energy (kcal/mol) PBE PBE-lg  Exper E 0.8, 1.0, 1.2 c lattice constant (A) Exper c 6.556 graphite has AB stacking (also show AA eclipsed graphite)

Universal PBE-ℓg Method UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations; A. K. Rappé, C. J. Casewit, K. S. Colwell, W. A. Goddard III, and W. M. Skiff; J. Am. Chem. Soc. 114, 10024 (1992) Derived C6/R6 parameters from scaled atomic polarizabilities for Z=1-103 (H-Lr) and derived Dvdw from combining atomic IP and C6 Universal PBE-lg: use same Re, C6, and De as UFF, add a single new parameter slg

blg Parameter Modifies Short-range Interactions 12-6 LJ potential (UFF parameter) blg =1.0 blg =0.7 lg potential lg potential When blg =0.6966, ELJ(r=1.1R0) = Elg(r=1.1R0)

Benzene Dimer T-shaped

Ch121a Atomic Level Simulations of Materials and Molecules