Introduction to Computational Chemistry

Introduction to Computational Chemistry NSF Computational Nanotechnology and Molecular Engineering Pan-American Advanced Studies Institutes (PASI) Workshop January 5-16, 2004 California Institute of Technology, Pasadena, CA Andrew S. Ichimura For the Beginner… There are three main problems: 1. Deciphering the language. 2. Technical implementation. 3. Quality assessment. Focus on… Calculating molecular structures and relative energies. 1. Hartree-Fock (Self-Consistent Field) 2. Electron Correlation 3. Basis sets and performance Molecular properties Transition States Reaction coords. Ab initio electronic structure theory Hartree-Fock (HF) Electron Correlation (MP2, CI, CC, etc.) Spectroscopic observables Geometry prediction Benchmarks for parameterization Prodding Experimentalists Goal: Insight into chemical phenomena. Setting up the problem… What is a molecule? A molecule is “composed” of atoms, or, more generally as a collection of charged particles, positive nuclei and negative electrons. The interaction between charged particles is described by; qi q j qi q j Vij  V (rij )   4 0 rij rij Coulomb Potential qi rij qj  Coulomb interaction between these charged particles is the only important physical force necessary to describe chemical phenomena. But, electrons and nuclei are in constant motion… In Classical Mechanics, the dynamics of a system (i.e. how the system evolves in time) is described by Newton’s 2nd Law: F = force F  ma a = acceleration 2 dV d r r = position vector  m 2 dr dt m = particle mass   In Quantum Mechanics, particle behavior is described in terms of a wavefunction, Y.Y  ˆ HY  i t  Time-dependent Schrödinger Equation i  1;  h 2    ˆ H   Hamiltonian Operator Time-Independent Schrödinger Equation ˆ ˆ H (r,t)  H (r) Y(r,t)  Y(r)eiEt /   If H is time-independent, the timedependence of Y may be separated out as a simple phase factor.  ˆ H (r)Y(r)  EY(r) Time-Independent Schrödinger Equation Describes the particle-wave duality of electrons.  Hamiltonian for a system with N-particles ˆ ˆ ˆ H  T V N N  N 2 2 2   2  2  2  ˆ ˆ T  Ti    i    2  2  2   y i zi  i1 2mi i1 2mi  x i  i1 Sum of kinetic (T) and potential (V) energy Kinetic energy    2  2  2     2  2  2  x i y i zi  2 i Laplacian operator  ˆ  V   qiq j V ij rij i1 j1 i1 j1 N N N N Potential energy  When these expressions are used in the time-independent Schrodinger Equation, the dynamics of all electrons and nuclei in a molecule or atom are taken into account. Born-Oppenheimer Approximation • So far, the Hamiltonian contains the following terms:  ˆ ˆ ˆ ˆ ˆ ˆ H  Tn  Te  Vne  Ve e  Vnn ˆ Tn ˆ Te ˆ Vne ˆ Vee ˆ Vnn Kinetic energy of nuclei, n Kinetic energy of electrons, e Electron-nuclear attraction Electron-electron repulsion Internuclear repulsion • Since nuclei are much heavier than electrons, their velocities are much smaller. To a good approximation, the Schrödinger equation can be separated into two parts: – One part describes the electronic wavefunction for a fixed nuclear geometry. – The second describes the nuclear wavefunction, where the electronic energy plays the role of a potential energy. Born-Oppenheimer Approx. cont. • In other words, the kinetic energy of the nuclei can be treated separately. This is the Born-Oppenheimer approximation. As a result, the electronic wavefunction depends only on the positions of the nuclei. Physically, this implies that the nuclei move on a potential energy surface (PES), which are solutions to the electronic Schrödinger equation. Under the BO approx., the PES is independent of the nuclear masses; that is, it is the same for isotopic molecules. E H + H • . . 0 H • H Solution of the nuclear wavefunction leads to physically meaningful quantities such as molecular vibrations and rotations. Limitations of the Born-Oppenheimer approximation • The total wavefunction is limited to one electronic surface, i.e. a particular electronic state. • The BO approx. is usually very good, but breaks down when two (or more) electronic states are close in energy at particular nuclear geometries. In such situations, a “ non-adiabatic” wavefunction - a product of nuclear and electronic wavefunctions - must be used. • In writing the Hamiltonian as a sum of electron kinetic and potential energy terms, relativistic effects have been ignored. These are normally negligible for lighter elements (Z<36), but not for the 4 th period or higher. • By neglecting relativistic effects, electron spin must be introduced in an ad hoc fashion. Spin-dependent terms, e.g., spin-orbit or spin-spin coupling may be calculated as corrections after the electronic Schrödinger equation has been solved. The electronic Hamiltonian becomes, ˆ ˆ ˆ ˆ ˆ H  Te  Vne   Ve e  Vnn B.O. approx.; fixed nuclear coord. Self-consistent Field (SCF) Theory GOAL: Solve the electronic Schrödinger equation, HeY=EY. PROBLEM: Exact solutions can only be found for one-electron systems, e.g., H2+. SOLUTION: Use the variational principle to generate approximate solutions. Variational principle - If an approximate wavefunction is used in HeY=EY, then the energy must be greater than or equal to the exact energy. The equality holds when Y is the exact wavefunction. In practice: Generate the “best” trial function that has a number of adjustable parameters. The energy is minimized as a function of these parameters. SCF cont. The energy is calculated as an expectation value of the Hamiltonian operator: ˆ  Y H Yd E  Y Yd  e  Introduce “bra-ket” notation, right Combined bracket denotes integration over all coordinates. ˆ  Y H Yd   Y Yd    e ˆ Y | He | Y bra ket n complex conjugate , left Y |Y m E ˆ Y |He |Y Y |Y If the wavefunctions are orthogonal and normalized (orthonormal), ij  1 Yi | Yj   ij (Kroenecker delta)  0 ij Then, ˆ E  Y | He | Y SCF cont. Antisymmetric wavefunctions can be written as Slater determinants. Since electrons are fermions, S=1/2, the total electronic wavefunction must be antisymmetric (change sign) with respect to the interchange of any two electron coordinates. (Pauli principle - no two electrons can have the same set of quantum numbers.) Consider a two electron system, e.g. He or H2. A suitable antisymmetric wavefunction to describe the ground state is: 1,2  1(1)2(2)  1 (2)2 (1) (He: 1 =2 = 1s) (H2: 1 = 2 = bonding MO) Each electron resides in a spin-orbital, a product of spatial and spin functions. (Spin functions are orthonormal:  |  =  |  =1;  |    |   0) Interchange the coordinates of the two electrons, 2,1  1 (2) 2 (1)  1 (1) 2 (2) 2,1   1,2 SCF cont. A more general way to represent antisymmetric electronic wavefunctions is in the form of a determinant. For the two-electron case, 1 (1)  2  (1) 1,2   1 (1) 2 (2)  1 (2) 2 (1) 1 (2)  2  (2) For an N-electron N-spinorbital wavefunction,  SD   1 1 1 2   2 (1)   N (1)  2 (2)   N (2)    ,  i |  j  ij 1 N   2 (N)   N (N) A Slater Determinant (SD) satisfies the antisymmetry requirement. are one-electron wavefunctions, molecular orbitals. Columns Rows contain the electron coordinates. One more approximation: The trial wavefunction will consist of a single SD. Now the variational principle is used to derive the Hartree-Fock equations... Hartree-Fock Equations (1) Reformulate the Slater Determinant as, ˆ ˆ   A1 (1) 2 (2) N (N) A    is the diagonal product ˆ A the antisymmetrizer 1 N 1 1  p ˆ ˆ ˆ  P   A 1  (1) P  N!   Pij  ˆijk  N! p 0  ij  ijk    ˆ ˆ P is the permutation operator. Pij permutes two electron coordinates . (2) ˆ ˆ ˆ ˆ ˆ He  Te  Vne   Ve e  Vnn 1 ˆ Te     2 i 2 i ˆ Vne     i a N N Za Ra  ri  ˆ Ve e    i j i N N One electron terms 1 ri  rj ZZ ˆ Vnn    a b a b a Ra  Rb Depends on two electrons 1 Za ˆ hi    2   i 2 a Ra  r i ˆ gij  1 ri  rj N N i N i 1 j i One-electron operator - describes electron i, moving in the field of the nuclei. Two-electron operator interelectron repulsion. Hamiltonian ˆ ˆ ˆ ˆ He   hi    gij  Vnn (3) Calculation of the energy. ˆ E e   | He |  ˆ ˆ ˆ ˆ ˆ E e  A | He | A    (1) p  | He | P p 0 N 1 Expectation value over Slater Determinant Examine specific integrals: ˆ  | Vnn |   Vnn Nuclear repulsion does not depend on electron coordinates. For coordinate 1, ˆ ˆ   | h1 |    1 (1) 2 (2) N (N)| h1 |  1 (1) 2 (2) N (N)   ˆ  1 (1) | h1 | 1 (1)  2 (2) |  2 (2)  N (N ) | N (N )  h1 The one-electron operator acts only on electron 1 and yields an energy, h1, that depends only on the kinetic energy and attraction to all nuclei. ˆ ˆ  | g12 |   1 (1) 2 (2) N (N)| g12 | 1(1) 2 (2) N (N ) ˆ  1 (1)2 (2) | g12 | 1 (1)2 (2)  3 (3) |  3 (3)  N (N ) |  N (N )  ˆ = 1 (1)2 (2) | g12 | 1 (1)2 (2)  J12 Coulomb integral, J12: represents the classical repulsion between two charge distributions 12(1) and 2 ˆ 2 (2). ˆ12 | P12  1 (1) 2 (2) N (N)| g12 |  2 (1)1 (2) N (N) ˆ  | g  ˆ  1 (1)2 (2) | g12 | 2 (1)1 (2)  3 (3) |  3 (3)  N (N ) |  N (N )  ˆ = 1 (1)2 (2) | g12 | 2 (1)1 (2)  K12 Exchange integral, K12: no classical analogue. Responsible for chemical bonds. The expression for the energy can now be written as: 1 N N E e   hi    (J ij  Kij )  Vnn 2 i j i 1 N Sum of one-electron, Coulomb, and exchange integrals, and Vnn. To apply the variational principle, the Coulomb and Exchange integrals are written as operators, Ee   i 1 N 1 N N ˆ ˆ ˆ  i | hi |  i     j | J i |  j   j | Ki |  j  Vnn 2 i j ˆ ˆ J i |  j (2)   i (1) | g12 |  i (1)  j (2) ˆ ˆ Ki |  j (2)   i (1) | g12 |  j (1)  i (2)   The objective now is to find the best orbitals ( i, MOs) that minimize the energy (or at least remain stationary with respect to further changes in i), while maintaining orthonormality of i. • Employ the method of Langrange Multipliers: f (x1 ,x 2 ,x N ) g(x1 ,x 2 ,x N )  0 L L Optimize L such that  0, 0 xi  i  L  E    ij  i |  j   ij ij N ij N Function to optimize. L(x1 ,x 2 ,x N , )  f (x1, x 2 , x N )  g(x1 ,x 2 ,x N ) Rewrite in terms of another function. Define Lagrange function. Constrained optimization of L. • In terms of molecular orbitals, the Langrange function is:    0 Change in L with respect to small changes in i should be zero. L  E    ij  i |  j  i |  j  • Change in the energy with respect changes in i. E   i 1 N  ˆ ˆ  i | hi |  i   i | hi | i    | Jˆ  Kˆ N i j ij j ˆ ˆ |  i  i | J j  K j |  i  Define the Fock Operator, Fi ˆ ˆ ˆ ˆ Fi  hi   J j  K j j N N   Effective one-electron operator, associated with the variation in the energy. ˆ ˆ E    i | Fi | i   i | Fi | i i 1 N   N ij Change in energy in terms of the Fock operator. ˆ ˆ L    i | Fi |  i   i | Fi |  i   ij  i |  j   i |  j i1     0 According to the variational principle, the best orbitals, i, will make L=0. After some algebra, the final expression N becomes: ˆ Fii   ij j Hartree-Fock Equations  j After a unitary transformation, ij0 and iii. ˆ Fii '  ii ' ˆ i   i '| Fi |  i ' HF equations in terms of Canonical MOs and diagonal Lagrange multipliers. Lagrange multipliers can be interpreted as MO energies. Note: 1. The HF equations cast in this way, form a set of pseudoeigenvalue equations. 2. A specific Fock orbital can only be determined once all the other occupied orbitals are known. 3. The HF equations are solved iteratively. Guess, calculate the energy, improve the guess, recalculate, etc. 4. A set of orbitals that is a solution to the HF equations are called Self-consistent Field (SCF) orbitals. 5. The Canonical MOs are a convenient set of functions to use in the variational procedure, but they are not unique from the standpoint of calculating the energy. Koopman’s Theorem The ionization energy is well approximated by the orbital energy, i. * Calculated according to Koopman’s theorem. Basis Set Approximation • • For atoms and diatomic molecules, numerical HF methods are available. In most molecular calculations, the unknown MOs are expressed in terms of a known set of functions - a basis set. Two criteria for selecting basis functions. I) They should be physically meaningful. ii) computation of the integrals should be tractable. • It is common practice to use a linear expansion of Gaussian functions in the MO basis because they are easy to handle computationally. • Each MO is expanded in a set of basis functions centered at the nuclei and are commonly called Atomic Orbitals. (Molecular orbital = Linear Combination of Atomic Orbitals - LCAO). MO Expansion i   ci   M M M LCAO - MO representation Coefficients are variational parameters ˆ Fi  c i    i  ci     HF equations in the AO basis FC  SC ˆ F   | F |   S    |   Matrix representation of HF eqns. Roothaan-Hall equations (closed shell) F - element of the Fock matrix S - overlap of two AOs Roothaan-Hall equations generate M molecular orbitals from M basis functions. • N-occupied MOs • M-N virtual or unoccupied MOs (no physical interpretation) Total Energy in MO basis E i 1 N 1 N N ˆ ˆ ˆ  i | hi |  i     i j | g |  i j   i j | g |  j i 2 i j  N  V nn Total Energy in AO basis E    c ic i i 1  N M 1 ˆ ˆ ˆ  | hi |      cic j c icj    | g |       | g |    2 ij  M   V nn One-electron integrals, M2 Two-electron integrals, M4 Computed at the start; do not change occ.MO Products of AO coeff form Density Matrix, D D  occ.MO  j c j cj ; D   i c ic i General SCF Procedure Obtain initial guess for coeff., ci,form the initial D Form the Fock matrix Two-electron integrals Iterate Diagonalize the Fock Matrix Form new Density Matrix Computational Effort • Formally, the SCF procedure scales as M4 (the number of basis functions to the 4th power). Accuracy • As the number of functions increases, the accuracy of the Molecular Orbitals improves. • As M, the complete basis set limit is reached  Hartree-Fock limit. • Result: The best single determinant wavefunction that can be obtained. (This is not the exact solution to the Schrodinger equation.) Practical Limitation • In practice, a finite basis set is used; the HF limit is never reached. • The term “Hartree-Fock” is often used to describe SCF calculations with incomplete basis sets. Restricted and Unrestricted Hartree-Fock Restricted Hartree-Fock (RHF) For even electron, closed-shell singlet states, electrons in a given MO with  and  spin are constrained to have the same spatial dependence. Restricted Open-shell Hartree-Fock (ROHF) The spatial part of the doubly occupied orbitals are restricted to be the same. Unrestricted Hartree-fock (UHF)  and  spinorbitals have different spatial parts.  5 4  Energy 3 2 1  Spinorbitals is(n) RHF Singlet ROHF Doublet UHF Doublet Comparison of RHF and UHF R(O)HF •  and  spins have same spatial part UHF •  and  spins have different spatial parts • Wavefunction, , is an eigenfunction of S2 operator. • For open-shell systems, the unpaired electron () interacts differently with  and  spins. The optimum spatial orbitals are different.  Restricted formalism is not suitable for spin dependent properties. • Starting point for more advanced calculations that include electron correlation. • Wavefunction is not an eigenfunction of S2.  may be contaminated with states of higher multiplicity (2S+1). • EUHF ≤ ER(O)HF • Yields qualitatively correct spin densities. • Starting point for more advanced calculations that include electron correlation. Ab Initio (latin, “from the beginning”) Quantum Chemistry Summary of approximations • • • • • • Born-Oppenheimer Approx. Non-relativistic Hamiltonian Use of trial functions, MOs, in the variational procedure Single Slater determinant Basis set, LCAO-MO approx. RHF, ROHF, UHF Consequence of using a single Slater determinant and the Self-consistent Field equations: Electron-electron repulsion is included as an average effect. The electron repulsion felt by one electron is an average potential field of all the others, assuming that their spatial distribution is represented by orbitals. This is sometimes referred to as the Mean Field Approximation. Electron correlation has been neglected!!!