Transcript Computational Quantum chemistry

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, [email protected]
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
1
Last Time
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
2
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
2
3
© copyrightbetween
2011 William A.densities
Goddard III, allr
rights
reserved\
JLecture
the Coulomb repulsion
ab is 1Ch121a-Goddard-L02
a=[ψ
a(1)] and rb
The energy for an antisymmetrized product,
Aψaψb
The total energy is that of the product wavefunction plus the new
terms arising from exchange term which is negative with 4 parts
Eex=-< ψaψb|h(1)|ψb ψa >-< ψaψb|h(2)|ψb ψa >-< ψaψb|1/R|ψb ψa >
- < ψaψb|1/r12|ψb ψa >
The first 3 terms lead to < ψa|h(1)|ψb><ψbψa >+
<ψa|ψb><ψb|h(2)|ψa >+ <ψa|ψb><ψb|ψa>/R
But <ψb|ψa>=0
Thus all are zero
Thus the only nonzero term is the 4th term:
-Kab=- < ψaψb|1/r12|ψb ψa > which is called the exchange energy
(or the 2-electron exchange) since it arises from the exchange
term due to the antisymmetrizer.
Summarizing, the energy of the Aψaψb wavefunction for H2 is
E = haa + hbb + (Jab –Kab) + 1/R
Lecture 1Ch121a-Goddard-L02
One© new
term from the antisymmetrizer
copyright 2011 William A. Goddard III, all rights reserved\
4
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.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
5
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 Cm Xm, 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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
6
New stuff
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
7
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.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
8
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
9
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
10
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
11
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. B 76, 6062 (1964).
min
Showed that there exists some functional of the density FHK [ r ] 
r = V - rep
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.
Lecture 1Ch121a-Goddard-L02
Walter Kohn (1923-)
Nobel
Chemistry 199812
© copyright 2011 William A. Goddard
III, allPrize
rights reserved\
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,..N SA=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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
13
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.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
14
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 + ʃd3r r(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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
15
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
KS
Lecture 1Ch121a-Goddard-L02
Thus H
© copyright 2011 William
A. Goddard III, all rights reserved\
looks quite analogous
to HHF
16
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
17
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
E xc = E x  Ec
Lecture 1Ch121a-Goddard-L02
Ex =  ε x ρ(r),  ρ(r) ,...]ρ(r)dr
© copyright 2011 William A. Goddard III, all rights reserved\
18
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
ε
LDA
x
(ρ ) = A x ρ(r )
Lecture 1Ch121a-Goddard-L02
1
3
1
3
3 3
A x = -   .
4π
© copyright 2011 William A. Goddard III, all rights reserved\
19
Generalized gradient approximations
E xc = E x  Ec
3.5
3.0
ε
(ρ,  ρ ) = ε
LDA
x
 F(s)
Becke 88
2.5
F(S)
Ex =  ε x ρ(r),  ρ(r) ,...]ρ(r)dr
GGA
x
F(s) GGA functionals
X3LYP
2.0
1.5
PBE
PW91
1.0
0.5
s=
F
B88
0.0
ρ
(24 π ) ρ
1
2 3
0.0
5.0
s
10.0
S
4
3
2
9
1  sa1sinh1 (sa 2 )  a 3sBecke
b = 0.0042 a4 and a5 zero
(s) =
1
1
1  sa1sinh (sa 2 )
2
(
1  sa1sinh (sa 2 )  a 3  a 4e
1
1
d a2
2
(sa2 ,2a)3=a5s 2
1π sa
a 1 = 6 βa
)3 ,1sinh
Here a 2 = (48
FPW 91 (s) =
1
Lecture 1Ch121a-Goddard-L02
100 s 2
2
a
10
2
3 , a = 6 βa , a = 
a
=

(
)
a
=
48
π
Here

β
,
4
1
2
2
3
1/ 3
2
81
2 Ax
s
)
 a 42  106
10
 β , a 4 =  a 3 , a 5 = 1/ 3
, and d = 4.
81 Goddard III, all
© copyright
21 / 3 2011
A x William A.
2 rights
A x reserved\
20
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:
Exc  r ] =  U xc,l  r ]d l
1
0
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, 99829985.
Mori-Sanchez, P., Cohen, A.J. and Yang, W.T. J. Chem. Phys. (2006), 124,
091102-1-4.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
21
Becke half and half functional
assume a linear model
take
U xc,l =0 = Exexact
partition
the exact exchange of the KS orbitals
U xc,l =1  U
approximate
set
U xc,l = a  bl
LDA
xc,l =1
ExcLDA = ExLDA  EcLDA
=xexact
Exexact
=xcLDA
ExcLDA
xexact
Exexact
a =a E
; b;=b E
E
Get half-and-half functional
1 exact
1 LDA
LDA
Exc  r ] = ( Ex  Ex )  Ec
2
2
Becke, A.D. J. Chem. Phys.
(1993), 98, 1372-1377
© copyright 2011 William A. Goddard III, all rights reserved\
Lecture 1Ch121a-Goddard-L02
22
Becke 3 parameter functional
Empirically modify half-and-half
ExcB3  r ] = ExcLDA  c1 ( Exexact  ExLDA )  c2 ExGGA  c3EcGGA
where ExGGA is the gradient-containing correction terms to the LDA exchange
EcGGA is the gradient-containing correction to the LDA correlation,
c1, c2 , c3 are constants fitted against selected experimental thermochemical data.
The success of B3LYP in achieving high accuracy demonstrates that errors of ExcDFT for
covalent bonding arise principally from the λ  0 or exchange limit, making it important
to introduce some portion of exact exchange
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, 99829985.
Mori-Sanchez, P., Cohen, A.J. and Yang, W.T. J. Chem. Phys. (2006), 124,
091102-1-4.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
23
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
24
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
M06 = HF  tPBE + VSXC
Lecture 1Ch121a-Goddard-L02
Hybrid meta-GGA
© copyright 2011 William A. Goddard III, all rights reserved\
25
XYG3 approach to include London Dispersion in DFT
Görling-Levy coupling-constant perturbation expansion
Exc  r ] =  U xc,l  r ]d l
1
Take initial slope as the
2nd
0
order correlation energy:
Sum over virtual orbtials
i j ˆee ab
U xc ,l =0 =
U xc,l
l
= 2EcGL 2
l =0
i ˆx  fˆ a
2
2
1
E =  
 
4 ij ab e i  e j  ea  e b
e i  ea
i
a
where ˆee is the electron-electron repulsion operator, ˆx is the local exchange operator,
and fˆ is the Fock-like, non-local exchange operator.
GL 2 exact
Substitute into
Ex or ; b = ExcLDA  Exexact
U xc,l = a  bl with b = 2aEc=
where
GL 2
c
Combine both approaches (2 choices for b)
b = b1 EcGL 2  b2 ( ExcDFT  Exexact )
ExcR5  r ] = ExcLDA  c1 ( Exexact  ExLDA )  c2 ExGGA  c3 ( EcPT 2  EcLDA )  c4 EcGGA
DFT
a double hybrid DFT that mixes some exact exchange into Ex
DFT
PT 2
certain portion of Ec
into Ec
1
GL 2
PT 2 contains the double-excitation parts of E
=
E
c
c
while also introducing a 2

4 ab
ij
i j ˆee ab
e i  e j  ea  e b
This is a fifth-rung functional (R5) using information from both occupied and virtual KS
26
Lecture 1Ch121a-Goddard-L02
© copyrightdispersion
2011 William A. Goddard III, all rights reserved\
orbitals.
In principle can now describe
Final form of XYG3 DFT
ExcR5  r ] = ExcLDA  c1 ( Exexact  ExLDA )  c2 ExGGA  c3 ( EcPT 2  EcLDA )  c4 EcGGA
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.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
27
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.
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
28
Accuracy (kcal/mol) of various QM methods for
predicting standard enthalpies of formation
Functional
MAD
Max(+)
Max(-)
XYG3 a
1.81
16.67 (SF6)
-6.28 (BCl3)
M06-2x a
2.93
20.77 (O3)
-17.39 (P4)
M06 a
4.17
11.25 (O3)
-25.89 (C2F6)
B2PLYP a
4.63
20.37(n-octane)
-8.01(C2F4)
B3LYP a
4.74
19.22 (SF6)
-8.03 (BeH)
M06-L a
5.82
14.75 (PF5)
-27.13 (C2Cl4)
BLYP b
9.49
41.0 (C8H18)
-28.1 (NO2)
PBE b
22.22
10.8 (Si2H6)
-79.7 (azulene)
LDA b
121.85
0.4 (Li2)
-347.5 (azulene)
HFa
211.48
582.72(n-octane)
-0.46 (BeH)
MP2a
10.93
29.21(Si(CH3)4)
-48.34 (C2F6)
QCISD(T) c
15.22
42.78(n-octane)
-1.44 (Na2)
7.2 (SiF4)
-9.4 (C2F6)
DFT
Ab initio
G2(1 empirical parm) 1.88
Lecture
1Ch121a-Goddard-L02
G3(1
empirical
parm) 1.05
© copyright
2011) William A. Goddard III, all-4.9
rights(C
reserved\
7.1 (PF
F)
29
Comparison of QM methods for reaction surface of
H + CH4  H2 + CH3
30.00
H + CH4  H2 + CH3
HF
25.00
Energy (kcal/mol)
(kcal/mol)
Energy
HF
HF_PT2
XYG3
20.00
CCSD(T)
B3LYP
15.00
CCSD(T)
XYG3
BLYP
SVWN
HF_PT2 SVWN
B3LYP
10.00
BLYP
5.00
SVWN
0.00
-2.00
-1.50
-5.00
Lecture 1Ch121a-Goddard-L02
-1.00
-0.50
0.00
0.50
1.00
1.50
2.00
2.50
ReactionR(CH)-R(HH)
coordinate (in Å)
Reaction Coordinate:
© copyright 2011 William A. Goddard III, all rights reserved\
30
All (76)
HT38
HAT12
NS16
UM10
XYG3
1.02
0.75
1.38
1.42
0.98
M06-2x a
1.20
1.13
1.61
1.22
0.92
B2PLYP
1.94
1.81
3.06
2.16
0.73
M06 a
2.13
2.00
3.38
1.78
1.69
M06-La
3.88
4.16
5.93
3.58
1.86
B3LYP
4.28
4.23
8.49
3.25
2.02
BLYP a
8.23
7.52
14.66
8.40
3.51
PBEa
8.71
9.32
14.93
6.97
3.35
LDAb
14.88
17.72
23.38
8.50
Zhao and Truhlar
compiled benchmarks
Ab initio
of accurate barrier
HFb
11.28
13.66
16.87
6.67
heights in 2004
MP2 b
4.57
4.14
11.76
0.74
includes forward and
reverse barrier heights
QCISD(T) b
1.10
1.24
1.21
1.08
for
19 hydrogen transfer (HT) reactions,
6 heavy-atom transfer (HAT) reactions,
8 nucleophilic substitution (NS) reactions and
5Lecture
unimolecular
and association© (UM)
reactions.
1Ch121a-Goddard-L02
copyright
2011 William A. Goddard III, all rights reserved\
5.90
Reaction
barrier
heights
Note: no reaction
barrier heights used
in fitting the 3
parameters in
XYG3)
Functional
DFT
3.82
5.44
0.53
31
A. Total Energy (kcal/mol)
30.00
Test for
London
Dispersion
B. Exchange Energy (kcal/mol)
25.00
B
20.00
15.00
30.00
Ex_B
Ex_B3LYP
Ex_XYG3
Ex_HF
Ex_S
10.00 B3LYP
25.00
5.00
BLYP
Energy (kcal/mol)
20.00
15.00
10.00
5.00
B3LYP
0.00
BLYP
B3LYP
XYG3
CCSD(T)
SVWN
HF_PT2
-5.00 3.0
0.00
-5.00
-10.00
-15.00
4.0
HF_PT2
5.0
S 4.0
HF XYG3
5.0
6.0
Distance (A)
VWN
-3.00 B3LYP
0.00
3.0
(B)
6.0
LDA
(SVWN)
Intermolecular distance
CCSD(T)
(A)
XYG3
Distance (A)
Ec_VWN
Ec_B3LYP
LYP
CCSD(T)
Ec_LYP
-6.00
XYG3
Ec_XYG3
Ec_CCSD(T)
-9.00
PT2
Ec_PT2
(C)
C.
Correlation
Energy
(kcal/mol)
-12.00
Conclusion: XYG3 provides excellent accuracy for London dispersion, as good as
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
CCSD(T)
32
Accuracy of QM methods for noncovalent interactions.
Functional
Note: no
noncovalent
complexes used
in fitting the 3
parameters in
XYG3)
Total
HB6/04
CT7/04
DI6/04
WI7/05
PPS5/05
M06-2x b
0.30
0.45
0.36
0.25
0.17
0.26
XYG3 a
0.32
0.38
0.64
0.19
0.12
0.25
M06 b
0.43
0.26
1.11
0.26
0.20
0.21
M06-L b
0.58
0.21
1.80
0.32
0.19
0.17
B2PLYP
0.75
0.35
0.75
0.30
0.12
2.68
B3LYP
0.97
0.60
0.71
0.78
0.31
2.95
1.17
0.45
2.95
0.46
0.13
1.86
1.48
1.18
1.67
1.00
0.45
3.58
3.12
4.64
6.78
2.93
0.30
0.35
2.08
2.25
3.61
2.17
0.29
2.11
0.64
0.99
0.47
0.29
0.08
1.69
0.57
0.90
0.62
0.47
0.07
0.95
DFT
HB: 6 hydrogen bond PBE c
complexes,
BLYP c
CT 7 charge-transfer
LDA c
complexes
Ab initio
DI: 6 dipole
interaction complexes, HF
WI:7 weak interaction MP2c
complexes,
QCISD(T) c
PPS: 5 pp stacking
complexes.
Lecture 1Ch121a-Goddard-L02
WI and PPS dominated
by London
© copyright 2011 William A. Goddard III, all rights
reserved\dispersion.33
Problem
XYG3 approach to include London Dispersion in DFT
Görling-Levy coupling-constant perturbation expansion
Exc  r ] =  U xc,l  r ]d l
1
Take initial slope as the
2nd
0
order correlation energy:
Sum over virtual orbtials
i j ˆee ab
U xc ,l =0 =
2
U xc,l
l
= 2EcGL 2
l =0
i ˆx  fˆ a
2
1
E =  
 
4 ij ab e i  e j  ea  e b
e i  ea
i
a
where ˆee is the electron-electron repulsion operator, ˆx is the local exchange operator,
and fˆ is the Fock-like, non-local exchange operator.
where
GL 2
c
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
Lecture
1Ch121a-Goddard-L02
XYGJ-OS
and XYGJ-OS
© copyright 2011 William A. Goddard III, all rights reserved\
34
XYGJ-OS method most accurate DFT
(including London Dispersion) at modest cost
Ying Zhang , Xin Xu, Goddard; P. Natl. Acad. Sci. 106 (13) 4963-4968 (2009)
Use Görling-Levy coupling-constant perturbation expansion
ExcR5  r ] = ExcLDA  c1 ( Exexact  ExLDA )  c2 ExGGA  c3 ( EcPT 2  EcLDA )  c4 EcGGA
DFT
a double hybrid DFT that mixes some exact exchange into Ex
DFT
PT 2
certain portion of Ec
into Ec
1
GL 2
PT 2 contains the double-excitation parts of E
=
E
c
while also introducing a 2

4 ab
c
ij
c4 = (1 – c3)  3 parameters
i j ˆee ab
e i  e j  ea  e b
Yousung
Jung
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
(
) (
)
SVWN
HF
S
LYP
VWN
PT 2
EXYGJ-OS
r
=
E

c
E

E

c
E

c
E

c

ExcXYG4-OS
r
=
E

c
E

E

c
E

c
E

c

E
 ] ] xc xc 1 (1 x x x )x ( 2 2c c 3 3c c ) 4 4 osEos
xc
XYG4-OS
SVWN
Lecture 1Ch121a-Goddard-L02
HF
S
LYP
VWN
© copyright 2011 William A. Goddard III, all rights reserved\
PT 2
35
XYG4-OS and XYG4-LOS timings
CPU (hours)
Timings XYGJ-OS and XYGJ-LOS for long XYG4-LOS
alkanes
200.0
160.0
80.0
200.0
40.0
XYGJ-OS
XYG4-OS
120.0
B3LYP
200.0
B3LYP
160.0
200.0
0.0
0
160.0
XYG3
80.0
U (hours)
CPU (hours)
CPU (hours)
CPU (hours)
XYG4-LOS
XYGJ-LOS
120.0
160.0
40.0
120.0
80.0
120.0
0.0
0
20
Lecture 1Ch121a-Goddard-L02
XYG4-OS
120.0
XYG3
XYG4-OS and XYG4-LO
XYG4-OS and XYG4-
XYG4-OS and XYG4-LO
20
40
60
XYG4-LOS
alkane chain l
XYGJ-OS
XYG4-OS
XYG4-LOS
B3LYP
XYG4-OS
XYG4-LOS
XYGJ-LOS
XYG3
B3LYP
XYG4-OS
80
100
120
XYG3
80.0
B3LYP
© alkane
copyright 2011
William
A. Goddard III, all rights reserved\
chain
length
40
60
36
Accuracy of Methods (Mean absolute deviations MAD, in eV)
Methods
HOF
(223)
IP
(38)
DFT methods
SPL (LDA)
5.484 0.255
BLYP
0.412 0.200
PBE
0.987 0.161
TPSS
0.276 0.173
B3LYP
0.206 0.162
PBE0
0.300 0.165
M06-2X
0.127 0.130
XYG3
0.078 0.057
0.072 0.055
XYGJ-lOS
MC3BB
0.165 0.120
B2PLYP
0.201 0.109
Wavefunction based methods
HF
9.171 1.005
MP2
0.474 0.163
G2
0.082 0.042
G3
0.046 0.055
EA
(25)
PA
(8)
BDE
(92)
NHTBH
(38)
HTBH
(38)
NCIE
(31)
All
(493)
0.311
0.105
0.102
0.104
0.106
0.128
0.103
0.080
0.084
0.175
0.090
0.276
0.080
0.072
0.071
0.061
0.057
0.092
0.070
0.067
0.046
0.067
0.754
0.292
0.177
0.245
0.226
0.155
0.069
0.068
0.033
0.111
0.124
0.542
0.376
0.371
0.391
0.202
0.154
0.056
0.056
0.049
0.062
0.090
0.775
0.337
0.413
0.344
0.192
0.193
0.055
0.033
0.038
0.036
0.078
0.140
0.063
0.052
0.049
0.041
0.031
0.013
0.014
0.015
0.023
0.023
2.771
0.322
0.562
0.250
0.187
0.213
0.096
0.065
0.056
0.123
0.143
1.148
0.166
0.057
0.049
0.133
0.084
0.058
0.046
0.104
0.363
0.078
0.047
0.397
0.249
0.042
0.042
0.582
0.166
0.054
0.054
0.098
0.028
0.025
0.025
4.387
0.338
0.068
0.046
HOF = heat of formation; IP = ionization potential,
EA = electron affinity, PA = proton affinity,
BDE = bond dissociation energy,
NHTBH, HTBH = barrier heights for reactions,
Lecture =
1Ch121a-Goddard-L02
© copyright
2011 William A. Goddard III, all rights reserved\
NCIE
the binding in molecular
clusters
37
Comparison of speeds
HOF
IP TimeEA
P
HTBH
NCIE
All
Methods
(223) C100
(38)
(38)
(31)
(493)
H202 (25)
C100H100 (8
DFT methods
SPL (LDA)
5.484 0.255
0.311
0.2
0.542
0.775
0.140
2.771
0.376
0.337
BLYP 0.063 0.322
0.412 0.200 0.105 0.0
0.371
0.413
0.052
0.562
PBE
0.987 0.161 0.102 0.0
0.391
0.344
TPSS 0.049 0.250
0.276 0.173 0.104 0.0
0.202
0.192
2.8
12.3 0.0
B3LYP 0.041 0.187
0.206 0.162
0.106
0.154
0.193
PBE0 0.031 0.213
0.300 0.165 0.128 0.0
0.056
0.055
M06-2X0.013 0.096
0.127 0.130 0.103 0.0
0.056
0.033
81.4 0.0
XYG3 0.014 0.065
0.078 200.0
0.057 0.080
7.8
46.4 0.0
0.049
0.038
0.015
0.056
0.072 0.055
0.084
XYGJ-lOS
0.062
0.036
MC3BB 0.023 0.123
0.165 0.120 0.175 0.0
0.090
0.078
0.143
Lecture 1Ch121a-Goddard-L02
copyright 2011 William
A. Goddard III,
all rights reserved\
B2PLYP©0.023
0.201
0.109
0.090 38
0.0
NHTBH
(38)
Lecture 1Ch121a-Goddard-L02
G3
G2
XYG4-OS
XYG3
B2PLYP
M06-L
M06-2x
G2-1 small molecules
G2-2
G3-3 Large molecules
M06
10.0
9.0
8.0
7.0
6.0
5.0
4.0
3.0
2.0
1.0
0.0
B3LYP
MAD (kcal/mol)
Heats of formation (kcal/mol)
© copyright 2011 William A. Goddard III, all rights reserved\
39
Truhlar NHTBH38/04 set and HTBH38/04 set
HAT12
NS16
UM10
HT38
20.0
15.0
10.0
Reaction barrier heights (kcal/mol)
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
XYG4-OS
XYG3
QCISD(T)
MP2
HF
LDA
PBE
0.0
BLYP
5.0
B3LYP
MAD (kcal/mol)
25.0
40
Truhlar NCIE31/05 set
Nonbonded interaction (kcal/mol)
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
XYG4-OS
XYG3
QCISD(T)
MP2
HF
LDA
PBE
HB6
CT7
DI6
WI7
PPS5
BLYP
7.0
6.0
5.0
4.0
3.0
2.0
1.0
0.0
B3LYP
MAD (kcal/mol)
8.0
41
Comparison of QM methods for reaction
surface of H + CH4  H2 + CH3
30.00
H + CH4  H2 + CH3
HF
25.00
Energy (kcal/mol)
(kcal/mol)
Energy
HF
HF_PT2
XYG3
20.00
CCSD(T)
B3LYP
15.00
CCSD(T)
XYG3
BLYP
SVWN
HF_PT2 SVWN
B3LYP
10.00
BLYP
5.00
SVWN
0.00
-2.00
-1.50
-5.00
Lecture 1Ch121a-Goddard-L02
-1.00
-0.50
0.00
0.50
1.00
1.50
2.00
2.50
ReactionR(CH)-R(HH)
coordinate (in Å)
Reaction Coordinate:
© copyright 2011 William A. Goddard III, all rights reserved\
42
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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
43
Current challenge in DFT calculation for energetic
materials
• 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 longrange1/R6 attraction of the London dispersion but at same cost as
standard DFT
EDFT D = EDFT  Edisp
Clg,ij
N
Elg = - 
r  dReij6
6
ij , i  j ij
C6 single parameter from QM-CC
d =1
R
Rei + Rek (UFF vdW
radii)
Lecture
© copyright
2011 William A. Goddard III, all rights reserved\
eik =1Ch121a-Goddard-L02
44
PBE-lg for benzene dimer
T-shaped
Sandwich
PBE-lg parameters
Elg = - 
Clg-CC=586.8, Clg-HH=31.14, Clg-HH=8.691
Clg,ij
N
r  dR
6
ij , i  j ij
Parallel-displaced
6
eij
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
45
2550–2555
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
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
sublimation at 0 K: Exp:11.295
kcal/mol;
PBE: A.
0.913;
PBE-lg:
6.762reserved\
Lecture of
1Ch121a-Goddard-L02
© copyright
2011 William
Goddard
III, all rights
46
Binding energy (kcal/mol)
DFT-lg description for graphite
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)
© copyright 2011 William A. Goddard III, all rights reserved\
Lecture 1Ch121a-Goddard-L02
47
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 (HLr) 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
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
48
blg Parameter Modifies Short-range Interactions
12-6 LJ potential (UFF parameter)
blg =1.0
lg potential
lg potential
blg =0.7
When blg =0.6966,
ELJ(r=1.1R0) = Elg(r=1.1R0)
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
49
Benzene Dimer
T-shaped
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
50
Benzene Dimer
Sandwich
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
51
Benzene Dimer
Paralleldisplaced
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
52
Parameter Optimization
Implemented in VASP 5.2.11

0.7012
0.6966
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
53
Hydrocarbon Crystals
• Sublimation energy (kcal/mol/molecule)
Molecules
PBE
PBE-ℓg
Exp.
Benzene
1.051
12.808
11.295
Naphthalene
2.723
20.755
20.095
28.356
27.042
4.308 3/cell)
• Anthracene
Cell volume (angstrom
Molecules
PBE
PBE-ℓg
Exp.
Benzene
511.81
452.09
461.11
Naphthalene
380.23
344.41
338.79
Anthracene
515.49
451.55
451.59
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
55
Simple Molecular Crystals
• Sublimation energy (kcal/mol/molecule)
Molecules
PBE
PBE-ℓg
Exp.
F2
0.27
1.38
2.19
Cl2
2.05
5.76
7.17
Br2
5.91
10.39
11.07
I2
8.56
14.47
15.66
O2
0.13
1.50
2.07
 Average
N2 error: 3.86 (PBE)
0.02and 0.96 (PBE-ℓg)1.22
 Maximal
CO error: 7.10 (PBE)
0.11and 1.90 (PBE-ℓg)1.54
1.78
CO2
Lecture 1Ch121a-Goddard-L02
1.99
4.37
2.08
6.27
© copyright 2011 William A. Goddard III, all rights reserved\
56
Simple Molecular Crystals
• Cell volume (angstrom3/cell)
Molecules
PBE
PBE-ℓg
Exp.
F2
126.47
126.32
128.24
Cl2
282.48
236.23
231.06
Br2
317.30
270.06
260.74
I2
409.03
345.13
325.03
O2
69.38
69.35
69.47
N2
180.04
179.89
179.91
CO
178.96
178.99
179.53
CO2
218.17
179.93
177.88
Lecture 1Ch121a-Goddard-L02
© copyright 2011 William A. Goddard III, all rights reserved\
57
Chem 121 - Applied Quantum Chemistry
Method:
·
Semi-Empirical, used for very big systems, or for rough approximations of
geometry (extended Huckel theory, CNDO/INDO, AM1, MNDO)
·
HF (Hartree Fock). Simplest Ab Initio method. Very cheap, fairly inaccurate
·
MP2 (Moeller-Plasset 2). Advanced version of HF. Usually not as cheap or as
accurate as B3LYP, but can function as a complement.
·
CASSCF (Complete Active Space, Self Consisting Field). Advanced version
of HF, incorporating excited states. Mainly used for jobs where photochemistry is
important. Medium cost, Medium Accuracy. Quite complicated to run…
·
QCISD (Quadratic Configuration Interaction Singles Doubles). Very
advanced version of HF. Very Expensive, Very accurate. Can only be used on
systems smaller than 10 heavy atoms.
·
CCSD (Coupled Cluster Singles Doubles). Very much like QCISD.
Density Functional Theory
LDA (local density approximation)
PW91, PBE
·
B3LYP (density functional theory). Cheap, Accurate.
Generally, B3LYP is the method of choice. If the system allows it, QCISD or CCSD
can be used. HF and/or MP2 can be used to verify the B3LYP results.
Lecture 1Lecture 2
60
Chem 121 - Applied Quantum Chemistry
Basis Set: What mathematical expressions are used to describe orbitals. In
general, the more advanced the mathematical expression, the more accurate
the wavefunction, but also more expensive calculation.
·
STO-3G - The ‘minimal basis set’. Not particularly accurate, but cheap and
robust.
·
3-21G - Smallest practical Basis Set.
·
6-31G - More advanced, i.e. more functions for both core and valence.
·
6-31G** - As above, but with ‘polarized functions’ added. Essentially
makes the orbitals look more like ‘real’ ones. This is the standard basis set
used, as it gives fairly good results with low cost.
·
6-31++G - As above, but with ‘diffuse functions’ added. Makes the orbitals
stretch out in space. Important to add if there is hydrogen bonding, pi-pi
interactions, anions etc present.
·
6-311++G** - As above, with even more functions added on… The more
stuff, the more accurate… But also more expensive. Seldom used, as the
increase in accuracy usually is very small, while the cost increases drastically.
·
Frozen Core: Basis sets used for higher row elements, where all the core
electrons are treated as one big frozen chunk. Only the valence electrons are
treated explicitly
Lecture 1Lecture 2
61
Chem 121 - Applied Quantum Chemistry
• Software packages
–
–
–
–
–
–
–
Jaguar
GAMESS
TurboMol
Gaussian
Spartan/Titan
HyperChem
ADF
Lecture 1Lecture 2
62
Chem 121 - Applied Quantum Chemistry
Running an actual calculation
– Determine the starting geometry of the
molecule you wish to study
– Determine what you’d like to find out
– Determine what methods are suitable and/or
affordable for the above calculation
– Prepare input file
– Run job
– Evaluate result
Lecture 1Lecture 2
63
Chem 121 - Applied Quantum Chemistry
Example: Good ol’ water
Starting geometry: water is bent, (~104º), a normal
O-H bond is ~0.96 Å. For illustration, however, we’ll
start with a pretty bad guess.
Simple Z-matrix:
O1
1.00 Å
H2 O1 1.00
H3 O1 1.00 H2 110.00
Lecture 1Lecture 2
1.00 Å
110º
64
Chem 121 - Applied Quantum Chemistry
What do we wish to find out?
How about the IR spectra?
What is a suitable method for this calculation?
Well, any, really, since it is so small. But 99% of the
time the answer to this question is “B3LYP/631G**” – a variant of density functional theory that
is the main workhorse of applied quantum
chemistry, with a standard basis set. Let’s go with
that.
Lecture 1Lecture 2
65
Chem 121 - Applied Quantum Chemistry
Actual jaguar input:
&gen
igeopt=1
ifreq=1
dftname=b3lyp
basis=6-31g**
&
&zmat
O1
H2 O1 0.95
H3 O1 0.95 H2 120.00
&
Lecture 1Lecture 2
66
Chem 121 - Applied Quantum Chemistry
Running time!
Jaguar calculates the wave function for the
atomic coordinates we provided
From the wave function it determines the
energy and the forces on the current geometry
Based on this, it determines in what direction it
should move the atoms to reach a better
geometry, i.e. a geometry with a lower energy
Lecture 1Lecture 2
67
Chem 121 - Applied Quantum Chemistry
Forces
Our horrible guess
1.00 Å
1.00 Å
Target geometry
0.96 Å
110º
0.96 Å
104º
Think elastic springs:
The bonds are too long,
so there will be a force
towards shorter bonds
Lecture 1Lecture 2
68
Chem 121 - Applied Quantum Chemistry
Optimization –
minimization of the
forces. When all forces
are zero the energy will
not change and we
have the resting
geometry
O1
H2 O1 0.9500000000
H3 O1 0.9500000000 H2 120.0000000000
SCF energy: -76.41367730925
-O1
H2 O1 0.9566666804
H3 O1 0.9566666820 H2 106.8986301461
SCF energy: -76.41937497895
-O1
H2 O1 0.9653619358
H3 O1 0.9653619375 H2 103.0739287925
SCF energy: -76.41969584939
-O1
H2 O1 0.9653155294
H3 O1 0.9653155310 H2 103.6688074046
SCF energy: -76.41970381840
--
Lecture 1Lecture 2
69
Chem 121 - Applied Quantum Chemistry
Accuracy
“actual” accuracy
Computer accuracy
0.9653155294 Å
0.96 Å
0.96 Å
0.9653155294 Å
103.6688074046º
103.7º
Accuracy is a relative concept
Lecture 1Lecture 2
70
Chem 121 - Applied Quantum Chemistry
frequencies
1666.01 3801.19 3912.97
No negative frequencies!
(Compare IR spectra for gas-phase water)
Lecture 1Lecture 2
71
Chem 121 - Applied Quantum Chemistry
Zero Point Energies
Vibrational levels
Zero Point Energy (ZPE)
“zero” level
Optimized energy is at the zero level, but in reality the molecule has a higher
energy due to populated vibrational levels.
At 0 K, all molecules populate the lowest vibrational level, and so the
difference between the “zero” level and the first vibrational level is the Zero
Point Energy (ZPE)
From our calculation:
The zero point energy (ZPE):
13.410 kcal/mol
Lecture 1Lecture 2
72
Chem 121 - Applied Quantum Chemistry
Thermodynamic data at higher temperatures
T = 298.15 K
trans.
rot.
vib.
elec.
total
U
--------0.889
0.889
0.002
0.000
1.779
Cv
--------2.981
2.981
0.041
0.000
6.003
S
--------34.609
10.503
0.006
0.000
45.117
H
--------1.481
0.889
0.002
0.000
2.371
G
---------8.837
-2.243
0.000
0.000
-11.080
Most thermodynamic data can be computed with very good
accuracy in the gas phase. Temperature dependant
Lecture 1Lecture 2
73
Chem 121 - Applied Quantum Chemistry
Transition states
Transition State (TS)
Stationary points:
points on the surface
where the derivative
of the energy = 0
Line represents the
reacting coordinate, in this
case the forming C-Cl and
breaking C-Br bonds
Product
Reactant
Reaction coordinate
CH3Br +
Cl-
TS
CH3Cl + Br-
Lecture 1Lecture 2
74
Chem 121 - Applied Quantum Chemistry
Transition state =
stationary point where all forces
except one is at a minimum.
Not a hill, but a
mountain pass
The exception is at its maximum
Reaction coordinate
CH3Br +
Cl-
TS
CH3Cl + Br-
Lecture 1Lecture 2
75
Chem 121 - Applied Quantum Chemistry
Derivative of the energy = 0
TS
Second derivative:
For a minimum > 0
For a maximum < 0
So a TS should have a
negative second derivative of
the energy
Second derivative of the
energy = force
Reactant
Lecture 1Lecture 2
Product
76
Chem 121 - Applied Quantum Chemistry
A transition state should have one
negative (imaginary) frequency!!!
(and ONLY one)
Lecture 1Lecture 2
77
Chem 121 - Applied Quantum Chemistry
Inflection points
Optimizing transition states:
Simultaneously optimize all
modes (forces) towards their
minimum, except the reacting
mode
But for the computer to know
which mode is the reacting
mode, you must have one
imaginary frequency in your
starting point
TS
Reactant
Product
Region with
imaginary frequency
Must start with a good guess!!!
Lecture 1Lecture 2
78
Chem 121 - Applied Quantum Chemistry
Example:
CH3Br + Cl-
CH3Cl + Br-
What do we know about this reaction? It’s an SN2
reaction, so the Cl- must come in from the backside of
the CH3Br. The C-Cl forms at the same time as the CBr forms. The transition state should be five
coordinate
Lecture 1Lecture 2
79
Chem 121 - Applied Quantum Chemistry
H
H
Br
C
Cl
2.2
2.0
H
Initial guess: C-Cl = 2.0 Å, C-Br = 2.2 Å
Single point frequency on the above geometry:
frequencies
frequencies
98.64 99.58 109.11 310.66 1339.10 1348.64
1349.46 1428.45 1428.73 2838.52 3017.70 3017.93
No negative frequencies! Bad initial guess
Lecture 1Lecture 2
80
Chem 121 - Applied Quantum Chemistry
Refinement :
Initial guess most likely wrong because of erronous CBr and C-Cl bond lengths
Let the computer optimize the five-coordinate structure
Frozen optimizations:
Just like a normal optimization, but with one or more
geometry parameters frozen
In this case, we optimize the structure with all the H-CCl angles frozen at 90º
Lecture 1Lecture 2
81
Chem 121 - Applied Quantum Chemistry
Result:
Cl
Br
2.32
2.62
C-Cl and C-Br bonds quite a bit longer in the new structure
Frequency calculation:
frequencies
frequencies
-286.26 168.54 173.32 173.43 874.16 874.76
976.23 1413.99 1414.65 3220.91 3420.84 3421.80
One negative frequency! Good initial guess
Lecture 1Lecture 2
82
Chem 121 - Applied Quantum Chemistry
Time for the actual optimization:
Jaguar follows the negative frequency towards the maximum
Geometry optimization 1: SCF Energy = -513.35042353681
Geometry optimization 2: SCF Energy = -513.34995058422
Geometry optimization 3: SCF Energy = -513.35001640704
Geometry optimization 4: SCF Energy = -513.34970196448
Geometry optimization 5: SCF Energy = -513.34968682825
Geometry optimization 6: SCF Energy = -513.34968118535
Final energy higher than starting energy (although only 0.5 kcal/mol)
Frequency calculation
frequencies
-268.67 162.64 174.22 174.31 848.15 848.24
frequencies
960.97 1415.75 1415.96 3220.77 3420.80 3421.15
One negative frequency! We found a true transition state
Lecture 1Lecture 2
83
Chem 121 - Applied Quantum Chemistry
Cl
Br
2.46
Final geometry:
2.51
C-Cl = 2.46 Å
C-Br = 2.51 Å
Cl-C-H = 88.7º
Br-C-H = 91.3º
Structure not quite symmetric, the
hydrogens are bending a little bit away
from the Br.
Lecture 1Lecture 2
84
Chem 121 - Applied Quantum Chemistry
Solvation calculations
Explicit solvents:
Calculations where solvent molecules
are added as part of the calculation
Implicit solvents:
Calculations where solvation effects
are added as electrostatic interactions
between the molecule and a virtual
continuum of “solvent”.
Lecture 1Lecture 2
85
Chem 121 - Applied Quantum Chemistry
Reaction energetics and barrier heights
Collect the absolute energies of the reactants, products and transition states
CH3Br +
Cl-
-53.078938 + -460.248741
CH3Cl + Br-
TS
-513.349681
-500.108371 + -13.237607
Sum each term
CH3Br +
Cl-
-513.327679
TS
CH3Cl + Br-
-513.349681
-513.345978
Define reactants as “0”, and deduct the reactant energy from all terms
CH3Br +
0
Cl-
TS
-.022002
CH3Cl + Br-.018299
Convert to kcal/mol (1 hartree = 627.51 kcal/mol)
Lecture 1Lecture 2
86
Chem 121 - Applied Quantum Chemistry
Reaction energetics and barrier heights
Convert to kcal/mol (1 hartree = 627.51 kcal/mol)
Cl-
CH3Br +
0
TS
-13.8
CH3Cl + Br-11.5
But this doesn’t make sense 
Lecture 1Lecture 2
87
Chem 121 - Applied Quantum Chemistry
Reaction energetics and barrier heights
Cl-
CH3Br +
0
TS
-13.8
CH3Cl + Br-11.5
Solvation not included!
Include solvation corrections!
Cl-
CH3Br +
0
TS
9.2
Lecture 1Lecture 2
CH3Cl + Br-6.4
88