Ch121a Atomic Level Simulations of Materials and Molecules

Ch121a Atomic Level Simulations of Materials and Molecules Room BI 115 Lecture: Monday, Wednesday Friday 2-3pm Lab Session: Lecture 2, April 2, 2014 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 316 Beckman Institute TA’s Caitlin Scott and Andrea Kirkpatrick Special Advice and Help: Julius Su/SKIES

Homework and Research Project First 5 weeks: The homework each week uses generally available computer software implementing the basic methods on applications aimed at exposing the students to understanding how to use atomistic simulations to solve problems. Each calculation requires making decisions on the specific approaches and parameters relevant and how to analyze the results. Midterm: each student submits proposal for a project using the methods of Ch121a to solve a research problem that can be completed in the final 5 weeks. The homework for the last 5weeks is to turn in a one page report on progress with the project The final is a research report describing the calculations and conclusions

Last Time Overview of Quantum Mechanics, Hydrogen Atom, etc Please review again to make sure that you are comfortable with the concepts, which you should have seen before

The Hartree Fock Equations General concept: there are an infinite number of possible orbitals for the electrons. For a system with 2M electrons we will put the electrons into the M lowest orbitals, with two electrons in each orbital (one up or a spin, the other down or b spin) M occ orb 2M elect

z The Hartree Fock EquationsClosed shell M occ orb N=2M elect b The wavefunction is written as Ψ(1,2,3,4, ..N-1,N) = A[(φaa)(φab)(φba)(φbb)---------(φza)(φzb)] Where the A is the antisymmetrizer or determinant operator where the 1st column is φaa(1), φaa(2), φaa(3), etc The 2nd column is φab(1), φab(2), φab(3), etc Thus there are N! terms This guarantees that the wavefunction changes sign if any 2 electrons are interchanged (Pauli Principle) Properties of determinant: if two columns are identical get zero. Thus can never have 2 electrons in same orbital with same spin Can take every column to be orthogonal; thus <φa|φb>=0 Also can recombine any two orbitals and the wavefunction does a = (cos) φa + (sin) φb not change b= (-sin) φa + (cos) φb a N 1 2

The energy (closed shell) The electronic Hamiltonian is Hel(1,2,---N) = Si h(i) + Si<j 1/rij where h(i) = - ½ 2+ SiZA/Raiis the interaction of all nuclei A with electron 1 plus the kinetic energy, a total of N terms and the other term is the Coulomb interaction between each pair of electrons, a total of N(N-1)/2 terms If we ignore the antisymmetrizer, so that the wavefunction is a Hartree product Ψ(1,2,3,4, ..N-1,N) = [(φaa)(φab)(φba)(φbb)---------(φza)(φzb)] Then the energy is Eproduct = Sa 2<a|h|a> + SaJaa+ Sa<b 4Jab Thus N=2M 1e terms and M+2M(M-1)=M(2M-2+1)= N(N-1)/2 2e terms M=N/2 terms M(M-1)/2 terms M=N/2 terms

The Coulomb energy • Jab= <Φa(1)Φb(2) |1/r12 |Φa(1)Φb(2)> • = ʃ1,2[a(1)]2 [b(2)]2/r12 • =ʃ1[a(1)]2Jb(1) • where Jb(1) = ʃ [b(2)]2/r12is the coulomb potential evaluated at point 1 due to the charge density [b(2)]2integrated over all space • Thus Jabis the total Coulomb interaction between the electron density ra(1)=|a(1)|2 and rb(2)=|b(2)|2 • Since the integrand ra(1) rb(2)/r12 is positive for all positions of 1 and 2, the integral is positive, Jab> 0

Consider the effect of the Antisymmetrizer Two electrons with same spin Ψ(1,2)) = A[(φaa)(φba)]= (φaa)(φba) - (φba)(φaa) 1 1 2 2 New term in energy is the exchange term -<(φaa)(φba)|Hel(1,2)|(φba)(φaa)> is a sum of 3 terms <(φaa)(φba)|h(1)|(φba)(φ1a)> = <φaa|h(1)|φba><φba|φaa> <(φaa)(φba)|h(2)|(φba)(φ1a)>=<φaa|φba><φba|h(2)|φaa> <(φaa)(φba)|1/r12|(φba)(φ1a)>=Kab Thus the only new term is -Kabnote that it is negative because one side is exchanged but not the other Thus the total energy becomes E = <a|h|a> + <b|h|b> + Jab – Kab 0 0

The Exchange energy Kab= <Φa(1)Φb(2) |1/r12 |Φb(1)Φa(2)> = ʃ1 [a(1)b(1)] ʃ2[b(2)a(2)]/r12 = ʃ1 [a(1) {P12 b(2)] ʃ2[b(2)]/r12 } a(1)] = ʃ1 [a(1) Kb(1) a(1)] No simple classical interpretation, but we have written it in terms of an integral operator Kb(1) so that is looks similar to the Coulomb case

Relationship between Jab and Kab The total electron-electron repulsion part of the energy for any wavefunction Ψ(1,2) = A[(φaa)(φba)] must be positive Eee =∫ (d3r1)((d3r2)|Ψ(1,2)|2/r12 > 0 This follows since the integrand is positive for all positions of r1 and r2then Thus Jab – Kab > 0 and hence Jab> Kab > 0 Thus the exchange energy is positive but smaller than the Coulomb energy Note that Kaa= <Φa(1)Φa(2) |1/r12 |Φa(1)Φa(2)> = Jaa

The final energy for closed shell wavefunction The total energy is E = Sa 2<a|h|a> + SaJaa + Sa<b(4Jab – 2Kab) One from aa and one from bb [2Jaa – Kaa] E= Sa2<a|h|a> + Sa [2Jaa – Kaa] + Sa<b (4Jab – 2Kab) E= Sa 2<a|h|a> + Sa,b (2Jab– Kab) There are M2 terms, so it appears that we have 2M2 = 2(N/2)(N/2) = N2/2 terms, but we should have N(N-1)/2 = N2 –N/2 This is because we added N/2 fake terms, Jaa that must be cancelled by the N/2 fake Kaa terms. Also note Sa,b2Jab = (½)ʃ1,2[r(1)][r(2)]2/r12 where r(1)=Sa [Φa(1)]2 is the total electron density, the classical electrostatic energy for this charge density

The Hartree Fock Equations • Variational principle: Require that each orbital be the best possible (leading to the lowest energy) leads to • HHF(1)φa(1)= eaφa(1) • where we solve for the occupied orbital, φa, to be occupied by both electron 1 and electron 2 • Here HHF(1)= h(1) + 2[Ja(1) - Ka(1)] • This looks like the Hamiltonian for a one-electron system in which the Hamiltonian has the form it would have neglecting electron-electron repulsion plus the average potential due the electron in the other orbital • Thus the two-electron problem is factored into M=N/2 one-electron problems, which we can solve to get φa, φb, etc

Self consistency • However to solve for φa(1)we need to know 2[Ja(1) - Ka(1)] • which depends on all M orbitals • Thus the HHFφa= eaφa equation must be solved iteratively until it is self consistent • But after the equations are solved self consistently, we can consider each orbital as the optimum orbital moving in the average field of all the other electrons • In fact the motions between these electrons would tend to be correlated so that the electrons remain farther apart than in this average field • Thus the error in the HF energy is called the correlation energy

He atom one slater orbital If one approximates each orbital as φ1s = N0exp(-zr) a Slater orbital then it is only necessary to optimize the scale parameter z In this case He atom: EHe = 2(½ z2) – 2Zz + (5/8)z Applying the variational principle, the optimum z must satisfy dE/dz = 0 leading to 2z- 2Z + (5/8) = 0 Thus z = (Z – 5/16) = 1.6875 KE = 2(½ z2) = z2 PE = - 2Zz + (5/8)z = -2 z2 E= - z2 = -2.8477 h0 Ignoring e-e interactions the energy would have been E = -4 h0 The exact energy is E = -2.9037 h0Thus this simple approximation of assuming that each electron is in a H1s orbital and optimizing the size accounts for 98.1% of the exact result.

The Koopmans orbital energy The next question is the meaning of the one-electron energy, eain the HF equation HHF(1)φa(1)= eaφa(1) Multiplying each side by φa(1) and integrating leads to ea <a|a> = <a|HHF|a> = <a|h|a> + 2Sb<a|Jb|a> - Sb<a|Kb|a> = <a|h|a> + Jaa + Sb≠a<a|2Jb-Kb|a> Thus in the approximation that the remaining electron does not change shape, eacorresponds to the energy to ionize an electron from the a orbital to obtain the N-1 electron system Sometimes this is referred to as the Koopmans theorem (pronounced with a long o). It is not really Koopmans theorem, which we will discuss later, but we will use the term anyway

The ionization potential There are two errors in using the ea to approximate the IP IPKT ~ -ea First the remaining N-1 electrons should be allowed to relax to the optimum orbital of the positive ion, which would make the Koopmans IP too large However the energy of the HF description is leads to a total energy less negative than the exact energy, Exact = EHF – Ecorr Where Ecorr is called the electron correlation energy (since HF does NOT allow correlation of the electron motions. Each electron sees the average potential of the other) which would make the Koopmans IP too small These effects tend to cancel so that the eafrom the HF wavefunction leads to a reasonable estimate of IP (N-1)e Ne HF from Ne exact Koopmans IP exact IP HF exact

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 We still refer to this as solving the HF equations

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.

Minimal Basis set – STO-3G For benzene the smallest possible basis set is to use a 1s-like single exponential function, exp(-zr) called a Slater function, centered on each the 6 H atoms and C1s, C2s, C2pz, C2py, C2pz functions on each of the 6 C atoms This leads to 42 basis functions to describe the 21 occupied MOs and is refered to as a minimal basis set. In practice the use of exponetial functions, such as exp(-zr), leads to huge computational costs for multicenter molecules and we replace these by an expansion in terms of Gaussian basis functions, such as exp(-ar2). The most popular MBS is the STO-3G set of Pople in which 3 gaussian functions are combined to describe each Slater function

Double zeta + polarization Basis sets – 6-31G** To allow the atomic orbitals to contract as atoms are brought together to form bonds, we introduce 2 basis functions of the same character as each of the atomic orbitals: Thus 2 each of 1s, 2s, 2px, 2py, and 2pz for C This is referred to as double zeta. If properly chosen this leads to a good description of the contraction as bonds form. Often only a single function is used for the C1s, called split valence In addition it is necessary to provide one level higher angular momentum atomic orbitals to describe the polarization involved in bonding Thus add a set of 2p basis functions to each H and a set of 3d functions to each C. The most popular such basis is referred to as 6-31G**

6-31G** and 6-311G** 6-31G** means that the 1s is described with 6 Gaussians, the two valence basis functions use 3 gaussians for the inner one and 1 Gaussian for the outer function The first *  use of a single d set on each heavy atom (C,O etc) The second *  use of a single set of p functions on each H The 6-311G** is similar but allows 3 valence-like functions on each atom. There are addition basis sets including diffuse functions (+) and additional polarization function (2d, f) (3d,2f,g), but these will not be relvent to EES810

Effective Core Potentials (ECP, psuedopotentials) • For very heavy atoms, say starting with Sc, it is computationally convenient and accurate to replace the inner core electrons with effective core potentials • For example one might describe: • Si with just the 4 valence orbitals, replacing the Ne core with an ECP or • Ge with just 4 electrons, replacing the Ni core • Alternatively, Ge might be described with 14 electrons with the ECP replacing the Ar core. This leads to increased accuracy because the • For transition metal atoms, Fe might be described with 8 electrons replacing the Ar core with the ECP. • But much more accurate is to use the small Ne core, explicitly treating the (3s)2(3p)6 along with the 3d and 4s electrons

Software packages • Jaguar: Good for organometallics • QChem: very fast for organics • Gaussian: many analysis tools • GAMESS • HyperChem • ADF • Spartan/Titan

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 + ½ ʃʃd3r1d3r2 [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 ʃ d3rr(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) General form of Energy for DFT (Kohn-Sham) formulation EKS = Si [<φi|- ½i2|φi >+Vext(ri)+Vxc(ri)]+½ʃʃd3r1d3r2[r(1)r(2)/r12] Coulomb repulsion Nuclear attraction KE Exchange correlation If the density is r =N/V then Coulomb repulsion leads to a total of ½(N/V)2 interactions, but it should be ½(N(N-1)/V2) Thus LDA include an extra self term that should not be present At the very minimum, Vxc needs to correct for this If density is uniform then error is proportional to 1/N. since electron density is r = N/V

The Local Density Approximation (LDA) ExcLDA[r(r)] = ʃ d3reXC(r(r)) r(r) where eXC(r(r)) is derived from Quantum Monte Carlo calculations for the uniform electron gas {DMCeperley 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 most serious 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

Generalized gradient approximations F(s) GGA functionals Becke 88 X3LYP PBE PW91 s S is big where the density gradient is large

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 the exact exchange of the KS orbitals take 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

Accuracy: DFT is basis for QM on catalysts Current flavors of DFT accurate for properties of many systems B3LYP and M06 useful for chemical reaction mechanisms Progress is being made on developing new systems

Accuracy: DFT is basis for QM on catalysts Current flavors of DFT accurate for properties of many systems B3LYP and M06 useful for chemical reaction mechanisms Example: Reductive elimination of CH4 from (PONOP)Ir(CH3)(H)+ Goldberg exper at 168Kbarrier DG‡ = 9.3 kcal/mol. These calculations use extended basis sets and PBF solvation DG(173K) B3LYP M06 M06L 0.0 0.0 0.0 10.8 5.8 11.4 (reductive elimination) • B3LYP and M06L perform well. • M06 underestimates the barrier.

Reductive Elimination Thermochemistry Mu-Jeng Cheng H/D exchange was measured from 153-173K by Girolami (J . Am. Chem. Soc., Vol. 120, 1998 6605) by NMR to have a barrier of DG‡ = 8.1 kcal/mol. DG(173K) B3LYP M06 0.0 0.0 8.7 9.5 (reductive elimination) 4.6 5.3 (s-bound complex) 6.4 5.2 (site-exchange) M06 and B3LYP functionals both consistent with experimental barrier site exchange. These calculations use extended basis sets and PBF solvation QM allows first principles predictions on new ligands, oxidation states, and solvents. But there are error bars in the QM having to do with details of the caculations (flavor of DFT, basis set). We use the best available methods and compare to any available experimental data on known systems to assess the accuracy for new systems. Some examples here and on the next slides

Reductive Elimination Thermochemistry B3LYP greatly underestimates the barrier since its repulsive non-bonding interactions underestimate the Pt-phosphine bond strength. M06L performs well and M06 underestimates the barrier. Reductive elimination of ethane from (dppe)Pt(CH3)4 was observed from 165-205˚C in benzene by Goldberg (J . Am. Chem. Soc., Vol. 125, 2003 9444) with a barrier of DG‡ = 36 kcal/mol (DS‡ = 15 e.u.). (As carbons are constrained to approach each other, the trans phosphine dissociates automatically.)

Metal-oxo Oxidations M06 performs well B3LYP overestimates bimolecular barriers involving bulky or polarizable species Phosphine oxidation by (Tp)Re(O)Cl2 and (Tpm)Re(O)Cl2+ was observed from 15-50˚C in 1,2-dichlorobenzene by Seymore and Brown (Inorg. Chem., Vol. 39, 2000, 325): Experiment: M06: B3LYP: DH‡(25C) 17.1 kcal/mol 16.6 24.1 Experiment: M06: B3LYP: DH‡(25C) 13.4 kcal/mol 11.8 17.1

Methods matter (must use the correct flavor DFT and the correct basis set) Commonly used methods (B3LYP, triple zeta basis set ) are insufficient for oxidation of main group elements. (Martin, J. Chem. Phys.1998,108(7), 2791.) B3LYP disfavors oxidation of main group elements by >10 kcal/mol Simple example Experimental DH (kcal/mol) -27 -80.1 S(CH3)2+ ½ O2→ O=S(CH3)2 S(CH3)2 + O2→ (CH3)2SO2 M06 6311G**++ -22.0 -70.7 M06 6311++G-3df(S) -29.2 -82.2 B3LYP 6311G**++ -17 -58.1 S(CH3)2+ ½ O2→ O=S(CH3)2 S(CH3)2 + O2→ (CH3)2SO2 OK Bad, but typical in publications

Methods matter (for reactions in polar media, must include solvation) Phosphine oxidation by (Tp)Re(O)Cl2 and (Tpm)Re(O)Cl2+ observed from 15-50˚C in 1,2-dichlorobenzene Seymore and Brown; Inorg. Chem., Vol. 39, 2000, 325) DH‡(25C) With solvation 17.1 16.6 24.1 Barrier with No solvation 16.9 Exper: M06: B3LYP: Barrier with No solvation 2.4 DH‡(25C) With solvation 13.4 11.8 17.1 Exper: M06: B3LYP: Much larger corrections in H2O Most QM publications ignore solvation or use unreliable methods

Essential issues: must include Solvation effects in the QM ThePoisson-Boltzmann Continuum Model in Jaguar/Schrödinger is extremely accurate Calculate Solvent Accessible Surface of the solute by rolling a sphere of radius Rsolv over the surface formed by the vdW radii of the atoms. Calculate electrostatic field of the solute based on electron density from the orbitals Calculate the polarization in the solvent due to the electrostatic field of the solute(need dielectric constant ) This leads to Reaction Field that acts back on solute atoms, which in turn changes the orbitals. Iterated until self-consistent. Calculate solvent forces on solute atoms Use these forces to determine optimum geometry of solute in solution. Can treat solvent stabilized zwitterions Difficult to describe weakly bound solvent molecules interacting with solute (low frequency, many local minima) Short cut: Optimize structure in the gas phase and do single point solvation calculation. Some calculations done this way Solvent:  = 99 Rsolv= 2.205 A PBF Implementation in Jaguar (Schrodinger Inc): pK organics to ~0.2 units, solvation to ~1 kcal/mol (pH from -20 to +20)

Comparison of PBF (Jaguar) pK with experiment 6.9(6.7) -3.89(-52.35) 5.8(5.8)-4.96(-49.64) 6.1(6.0) -3.98(-55.11) 5.3(5.3) -3.90(-57.94) 5.0(4.9) -4.80(-51.84) pKa: Jaguar(experiment) E_sol:zero(H+)

PBF (Jaguar) predictions of Metal-aquo pKa’s Protonated Complex (diethylenetriamine)Pt(OH2)2+ PtCl3(OH2)1- Pt(NH3)2(OH2)22+ Pt(NH3)2(OH)(OH2)1+ cis-(bpy)2Os(OH)(H2O)1+ Experimental pKa 6.3 7.1 5.5 7.4 11.0 Calculated (B3LYP) pKa(MAD: 1.1) 5.5 4.1 5.2 6.5 11.3 Calculated (M06//B3LYP) pKa (MAD: 1.6) 9.1 8.8 6.2 10.9 13.0 15.2 11.0 13.9 5.6 6.3 10.9 cis-(bpy)2Os(H2O)22+ cis-(bpy)2Os(OH)(H2O)1+ trans-(bpy)2Os(H2O)22+ trans-(bpy)2Os(OH)(H2O)1+ cis-(bpy)2Ru(H2O)22+ cis-(bpy)2Ru(OH)(H2O)1+ trans-(bpy)2Ru(H2O)22+ trans-(bpy)2Ru(OH)(H2O)1+ (tpy)Os(H2O)32+ (tpy)Os(OH)(H2O)21+ (tpy)Os(OH)2(H2O) Experimental pKa 7.9 11.0 8.2 10.2 8.9 >11.0 9.2 >11.5 6.0 8.0 11.0

Use theory to predict optimal pH for each catalyst 34.6 40.0 34.6 32.6 LnOsII(OH)3- is stable LnOsII(OH2)3+2 is stable LnOsII(OH2)(OH)2 is stable 37.9 Predict the relative free energies of possible catalyst resting states and transition states as a function of pH. Insertion transition states Resting states Optimum pH is 8

Predict PourbaixDiagrams to determine the oxidation states of transition metal complexes as function of pH and electrochemical potential Black experimental data from Dobson and Meyer, Inorg. Chem. Vol. 27, No.19, 1988. Red is from QM calculation (no fitting) using M06 functional, PBF implicit solvent Max errors: 200 meV, 2pH units Trans-(bpy)2Ru(OH)2 This is essential in using theory to predict new catalysts

Fundamental philosophy of First principles predictions QM calculations on small systems ~100 atoms get accurate energies, geometries, stiffness, mechanisms Fit QM to force field to describe big systems (104 -107atoms) Fit to obtain parameters for continuum systems macroscopic properties based on first principles (QM) Can predict novel materials where no empirical data available.

Ch121a Atomic Level Simulations of Materials and Molecules