Transcript Computational Chemistry G. H. CHEN Department of Chemistry University of Hong Kong

Computational Chemistry
G. H. CHEN
Department of Chemistry
University of Hong Kong
Beginning of Computational Chemistry
In 1929, Dirac declared, “The underlying physical
laws necessary for the mathematical theory of ...the
whole of chemistry are thus completely know, and
the difficulty is only that the exact application of
these laws leads to equations much too complicated
to be soluble.”
Dirac
Computational Chemistry
Quantum Chemistry
SchrÖdinger Equation
Molecular Mechanics
F=Ma
Nobel Prizes for Computational Chemsitry
Mulliken,1966
Pople, 1998
Fukui, 1981
Hoffmann, 1981
Kohn, 1998
Computational Chemistry Industry
Company
Software
Gaussian Inc.
Gaussian 94, Gaussian 98
Schrödinger Inc.
Jaguar
Wavefunction
Spartan
Q-Chem
Q-Chem
Accelrys
InsightII, Cerius2
HyperCube
HyperChem
Celera Genomics (Dr. Craig Venter, formal Prof., SUNY, Baffalo; 98-01)
Applications: material discovery, drug design & research
Carbon Nanotubes (Ijima, 1991)
Calculated STM Image
of a Carbon Nanotube
(Rubio, 1999)
STM Image of Carbon Nanotubes (Wildoer et. al., 1998)
Computer Simulations (Saito, Dresselhaus, Louie et. al.,
1992)
Carbon Nanotubes (n,m):
Conductor,
if n-m = 3I
Semiconductor, if n-m  3I
I=0,±1,±2,±3,…;or
Metallic Carbon Nanotubes:
Semiconducting Nanotubes:
Conducting Wires
Transistors
Molecular-scale circuits !
1 nm transistor!
30 nm transistor!
0.13 µm transistor!
Experimental Confirmations:
Lieber et. al. 1993; Dravid et. al., 1993; Iijima et. al. 1993;
Smalley et. al. 1998; Haddon et. al. 1998; Liu et. al. 1999
Wildoer, Venema, Rinzler, Smalley, Dekker, Nature 391, 59 (1998)
RL
L
Rc
C
7.39 kΩ
16.6 pH
6.45 kΩ (0.5g0-1)
0.073 aF
g0=2e2/h
Yam, Mo, Wang, Li, Chen, Zheng, Goddard (2008)
h
L ~   R L ~  d
 
h
2
e
≈ 18.8 pH
Microelectromechanical Systems (MEMS)
Micro-Electro-Mechanical Systems (MEMS) is the integration of mechanical elements,
sensors, actuators, and electronics on a common silicon substrate through
microfabrication technology. While the electronics are fabricated using integrated circuit
(IC) process sequences (e.g., CMOS, Bipolar, or BICMOS processes), the
micromechanical components are fabricated using compatible "micromachining"
processes that selectively etch away parts of the silicon wafer or add new structural layers
to form the mechanical and electromechanical devices.
Nanoelectromechanical Systems (NEMS)
K.E. Drexler, Nanosystems: Molecular Machinery,
anufacturing and Computation (Wiley, New York, 1992).
Large Gear Drives Small Gear
G. Hong et. al., 1999
Nano-oscillators
Nanoscopic Electromechanical Device
(NEMS)
Zhao, Ma, Chen & Jiang, Phys. Rev. Lett. 2003
Computer-Aided Drug Design
Human Genome Project
GENOMICS
ALDOSE REDUCTASE
HO
HO
OH HO
OH HO
Aldose Reductase
Diabetes
HO
O
HO
glucose
OH
Glucose
NADPH
NADP
HO
sorbitol
OH
Sorbitol
Diabetic
Complications
Design of Aldose Reductase Inhibitors
Inhibitor
Aldose Reductase
Hu, Chen & Chau, J. Mol. Graph. Mod. 24 (2006)
Prediction Results using AutoDock
LogIC50: 0.77,1.1
LogIC50: -1.87,4.05
LogIC50: -2.77,4.14
LogIC50: 0.68,0.88
Hu, Chen & Chau, J. Mol. Graph. Mod. 24 (2006)
Computer-aided drug design
Chemical Synthesis
Screening using in vitro assay
Animal Tests
Clinical Trials
Quantum Chemistry
G. H. Chen
Department of Chemistry
University of Hong Kong
Contributors:
Hartree, Fock, Slater, Hund, Mulliken,
Lennard-Jones, Heitler, London, Brillouin,
Koopmans, Pople, Kohn
Application:
Chemistry, Condensed Matter Physics,
Molecular Biology, Materials Science,
Drug Discovery
Emphasis
Hartree-Fock method
Concepts
Hands-on experience
Text Book
“Quantum Chemistry”, 4th Ed.
Ira N. Levine
http://yangtze.hku.hk/lecture/chem3504-3.ppt
Quantum Chemistry Methods
• Ab initio molecular orbital methods
• Semiempirical molecular orbital methods
• Density functional method
SchrÖdinger Equation
Hy=Ey
Wavefunction
Hamiltonian
H = (-h2/2m)2 - (h2/2me)ii2
+  ZZe2/r - i  Ze2/ri
+ i j e2/rij
Energy
Contents
1. Variation Method
2. Hartree-Fock Self-Consistent Field Method
3. Beyond Hartree-Fock
4. Perturbation Theory
5. Molecular Dynamics
The Variation Method
The variation theorem
Consider a system whose Hamiltonian operator
H is time independent and whose lowest-energy
eigenvalue is E1. If f is any normalized, wellbehaved function that satisfies the boundary
conditions of the problem, then
 f* H f d > E1
Proof:
Expand f in the basis set { yk}
f = k kyk
where
{k} are coefficients
Hyk = Ekyk
then
 f* H f d = k j k* j Ej dkj
= k | k|2 Ek > E 1 k | k|2 = E1
Since is normalized,
 f*f d = k | k|2 = 1
i. f : trial function is used to evaluate the upper limit
of ground state energy E1
ii. f = ground state wave function,  f* H f d = E1
iii. optimize paramemters in f by minimizing
 f* H f d /  f* f d
Application to a particle in a box of infinite depth
0
l
Requirements for the trial wave function:
i. zero at boundary;
ii. smoothness  a maximum in the center.
Trial wave function: f = x (l - x)
 f* H f dx = -(h2/82m)  (lx-x2) d2(lx-x2)/dx2 dx
= h2/(42m)  (x2 - lx) dx
= h2l3/(242m)
 f*f dx =  x2 (l-x)2 dx = l5/30
Ef = 5h2/(42l2m)  h2/(8ml2) = E1
Variational Method
(1) Construct a wave function f(c1,c2,,cm)
(2) Calculate the energy of f:
Ef  Ef(c1,c2,,cm)
(3) Choose {cj*} (i=1,2,,m) so that Ef is
minimum
Example: one-dimensional harmonic oscillator
Potential: V(x) = (1/2) kx2 = (1/2) m2x2 = 22m2x2
Trial wave function for the ground state:
f(x) = exp(-cx2)
 f* H f dx = -(h2/82m)  exp(-cx2) d2[exp(-cx2)]/dx2 dx
+ 22m2  x2 exp(-2cx2) dx
= (h2/42m) (c/8)1/2 + 2m2 (/8c3)1/2
 f*f dx =  exp(-2cx2) dx = (/2)1/2 c-1/2
Ef = W = (h2/82m)c + (2/2)m2/c
To minimize W,
0 = dW/dc = h2/82m - (2/2)m2c-2
c = 22m/h
W = (1/2) h
Extension of Variation Method
.
.
.
E3
y3
E2
y2
E1
y1
For a wave function f which is orthogonal to
the ground state wave function y1, i.e.
d f*y1 = 0
Ef = d f*Hf / d f*f > E2
the first excited state energy
The trial wave function f:
d f*y1 = 0
f = k=1 ak yk
d f*y1 = |a1|2 = 0
Ef = d f*Hf / d f*f = k=2|ak|2Ek / k=2|ak|2
> k=2|ak|2E2 / k=2|ak|2 =
E2
Application to H2+
e
f=c y +c y
1 1
2 2
+
+
y
1
y
2
W =  f*H f d /  f*f d
= (c12 H11 + 2c1 c2 H12 + c22 H22 )
/ (c12 + 2c1 c2 S + c22 )
W (c12 + 2c1 c2 S + c22) = c12 H11 + 2c1 c2 H12 + c22 H22
Partial derivative with respect to c1 (W/c1 = 0) :
W (c1 + S c2) = c1H11 + c2H12
Partial derivative with respect to c2 (W/c2 = 0) :
W (S c1 + c2) = c1H12 + c2H22
(H11 - W) c1 + (H12 - S W) c2 = 0
(H12 - S W) c1 + (H22 - W) c2 = 0
To have nontrivial solution:
H11 - W
H12 - S W
H12 - S W
=
0
H22 - W
For H2+, H11 = H22; H12 < 0.
Ground State: Eg = W1 = (H11+H12) / (1+S)
f1 = (y1+y2) / 2(1+S)1/2
bonding orbital
Excited State: Ee = W2 = (H11-H12) / (1-S)
f2 = (y1-y2) / 2(1-S)1/2
Anti-bonding orbital
Results: De = 1.76 eV, Re = 1.32 A
Exact: De = 2.79 eV, Re = 1.06 A
1 eV = 23.0605 kcal / mol
Further Improvements
Optimization of 1s orbitals
H
-1/2 exp(-r)
He+
23/2 -1/2 exp(-2r)
Trial wave function: k3/2 -1/2 exp(-kr)
Eg = W1(k,R)
at each R, choose k so that W1/k = 0
Results:
De = 2.36 eV, Re = 1.06 A
Inclusion of other atomic orbitals
1s
2p
Resutls: De = 2.73 eV, Re = 1.06 A
Linear Equations
1. two linear equations for two unknown, x1 and x2
a11x1 + a12x2 = b1
a21x1 + a22x2 = b2
(a11a22-a12a21) x1 = b1a22-b2a12
(a11a22-a12a21) x2 = b2a11-b1a21
Introducing determinant:
a11 a12
= a11a22-a12a21
a21 a22
a11 a12
b1
a12
a21 a22
b2
a22
a11 a12
a11 b1
x1 =
x2 =
a21 a22
a21 b2
Our case: b1 = b2 = 0, homogeneous
1. trivial solution: x1 = x2 = 0
2. nontrivial solution:
a11 a12
=0
a21 a22
n linear equations for n unknown variables
a11x1 + a12x2 + ... + a1nxn=
b1
a21x1 + a22x2 + ... + a2nxn=
b2
a11 a12
a21 a22
det(aij) xk= .
.
an1 an2
...
...
...
...
a1,k-1
a2,k-1
.
an,k-1
...
...
...
...
a1n
a2n
.
ann
where,
det(aij) =
a11 a12
a21 a22
.
.
an1 an2
b1 a1,k+1
b2 a2,k+1
.
.
b2 an,k+1
... a1n
... a2n
... .
... ann
inhomogeneous case: bk = 0 for at least one k
a11 a12
a21 a22
.
.
an1 an2
...
...
...
...
a1,k-1
a2,k-1
.
an,k-1
xk =
det(aij)
b1
b2
.
b2
a1,k+1 ...
a2,k+1 ...
. ...
an,k+1 ...
a1n
a2n
.
ann
homogeneous case: bk = 0, k = 1, 2, ... , n
(a) travial case: xk = 0, k = 1, 2, ... , n
(b) nontravial case: det(aij) = 0
For a n-th order determinant,
n
det(aij) =  alk Clk
l=1
where, Clk is called cofactor
Trial wave function f is a variation function
which is a combination of n linear independent
functions { f1 , f2 , ... fn},
f = c1f1 + c2f2 + ... + cnfn
n
  [( Hik - SikW ) ck ] = 0
k=1
Sik  d fi fk
Hik   d fi H fk
W   d f H f /  d f f
i=1,2,...,n
Linear variational theorem
(i) W1  W2  ...  Wn are n roots of Eq.(1),
(ii) E1  E2  ...  En  En+1  ... are energies
of eigenstates;
then, W1  E1, W2  E2, ..., Wn  En
Molecular Orbital (MO):
f = c1y1 + c2y2
( H11 - W ) c1 + ( H12 - SW ) c2 = 0
S11=1
( H21 - SW ) c1 + ( H22 - W ) c2 = 0
S22=1
Generally : yi a set of atomic orbitals, basis set
LCAO-MO
f = c1y1 + c2y2 + ...... + cnfn
linear combination of atomic orbitals
n
 ( Hik - SikW ) ck = 0
k=1
Hik   d yi* H yk
i = 1, 2, ......, n
Sik  d yi*yk
Skk = 1
The Born-Oppenheimer Approximation
Hamiltonian
H = (-h2/2m)2 - (h2/2me)ii2
+  ZZe2/r - i  Ze2/ri
+ i j e2/rij
H y(ri;r) = E y(ri;r)
The Born-Oppenheimer Approximation:
(1) y(ri;r) = yel(ri;r) yN(r)
(2) Hel(r )= - (h2/2me)ii2 - i Ze2/ri
+ ij e2/rij
VNN =  ZZe2/r
Hel(r) yel(ri;r) = Eel(r) yel(ri;r)
(3) HN = (-h2/2m)2 + U(r)
U(r) = Eel(r) + VNN
HN(r) yN(r) = E yN(r)
Assignment
Calculate the ground state energy and bond length of H2
using the HyperChem with the 6-31G
(Hint: Born-Oppenheimer Approximation)
Hydrogen Molecule H2
e
+
+
e
The Pauli principle
two electrons cannot be in the same state.
Wave function:
f(1,2) = ja(1)jb(2) + c1 ja(2)jb(1)
f(2,1) = ja(2)jb(1) + c1 ja(1)jb(2)
Since two wave functions that correspond to the same state
can differ at most by a constant factor
f(1,2) = c2 f(2,1)
ja(1)jb(2) + c1ja(2)jb(1) =c2ja(2)jb(1) +c2c1ja(1)jb(2)
c1 = c2
Therefore: c1 = c2 = 1
c2c1 = 1
According to the Pauli principle,
c1 = c2 =- 1
The Pauli principle (different version)
the wave function of a system of electrons must
be antisymmetric with respect to interchanging
of any two electrons.
Wave function f of H2 : Slater Determinant
y(1,2) = 1/2! [f(1)(1)f(2)(2) - f(2)(2)f(1)(1)]
=
1/2!
f(1)(1)
f(2)(2)
f(1)(1)
f(2)(2)
Energy: Ey
Ey=2 d1 f*(1) (Te+VeN) f(1) + VNN
+  d1 d2 |f2(1)| e2/r12 |f2(2)|
= i=1,2 fii + J12 + VNN
To minimize Ey under the constraint  d |f2| = 1,
use Lagrange’s method:
L = Ey - 2 e [ d1 |f2(1)| - 1]
dL = dEy - 4 e  d1 f*(1)df(1)
= 4 d1 df*(1)(Te+VeN)f(1)
+4 d1 d2 f*(1)f*(2) e2/r12 f(2)df(1)
- 4 e  d1 f*(1)df(1)
=0
[ Te+VeN + d2 f*(2) e2/r12 f(2) ] f(1) = e f(1)
Average Hamiltonian
Hartree-Fock equation
(f+J)f=ef
f(1) = Te(1)+VeN(1)
J(1) = d2 f*(2) e2/r12 f(2)
one electron operator
two electron Coulomb operator
f(1) is the Hamiltonian of electron 1 in the absence
of electron 2;
J(1) is the mean Coulomb repulsion exerted on
electron 1 by 2;
e is the energy of orbital f.
LCAO-MO:
f = c1y1 + c2y2
Multiple y1 from the left and then integrate :
c1F11 + c2F12 = e (c1 + S c2)
Multiple y2 from the left and then integrate :
c1F12 + c2F22 = e (S c1 + c2)
where,
Fij =  d yi* ( f + J ) yj = Hij +  d yi* J yj
S =  d y1 y2
(F11 - e) c1 + (F12 - S e) c2 = 0
(F12 - S e) c1 + (F22 - e) c2 = 0
Secular Equation:
F11 - e F12 - S e
F12 - Se F22 - e
bonding orbital:
= 0
e1 = (F11+F12) / (1+S)
f1 = (y1+y2) / 2(1+S)1/2
antibonding orbital: e2 = (F11-F12) / (1-S )
f2 = (y1-y2) / 2(1-S)1/2
Molecular Orbital Configurations of
Homo nuclear Diatomic Molecules H2, Li2, O, He2, etc
Moecule
H2+
H2
He2+
He2
Li2
Be2
C2
N2+
N2
O2+
O2
Bond order

1

0
1
0
2

3
2
2
De/eV
2.79
4.75
1.08
0.0009
1.07
0.10
6.3
8.85
9.91
6.78
5.21
The more the
Bond Order is,
the stronger
the chemical
bond is.
Bond Order:
one-half the difference
between the number of
bonding and antibonding
electrons
----------------
f1
----------------
f2
f1(1)(1)
f2(1)(1)
y(1,2) = 1 /2
f1(2)(2)
f2(2)(2)
= 1/2 [f1(1) f2(2) - f2(1) f1(2)] (1) (2)
Ey =  d1d2 y* H y
=  d1d2 y* (T1+V1N+T2+V2N+V12+VNN) y
= <f1(1)| T1+V1N|f1(1)>
+ <f2(2)| T2+V2N|f2(2)>
+ <f1(1) f2(2)| V12 | f1(1) f2(2)>
- <f1(2) f2(1)| V12 | f1(1) f2(2)> + VNN
= i <fi(1)| T1+V1N |fi(1)>
+ <f1(1) f2(2)| V12 | f1(1) f2(2)>
- <f1(2) f2(1)| V12 | f1(1) f2(2)> + VNN
= i=1,2 fii + J12 - K12 + VNN
Average Hamiltonian
Particle One:
Particle Two:
f(1) + J2(1) - K2(1)
f(2) + J1(2) - K1(2)
f(j)  -(h2/2me)j2 -  Z/rj
Jj(1) q(1)  q(1)  dr2 fj*(2) e2/r12 fj(2)
Kj(1) q(1)  fj(1) dr2 fj*(2) e2/r12 q(2)
Hartree-Fock Equation:
[ f(1)+ J2(1) - K2(1)] f1(1) = e1 f1(1)
[ f(2)+ J1(2) - K1(2)] f2(2) = e2 f2(2)
Fock Operator:
F(1)  f(1)+ J2(1) - K2(1)
F(2)  f(2)+ J1(2) - K1(2)
Fock operator for 1
Fock operator for 2
Summary
1. At the Hartree-Fock Level there are two possible
Coulomb integrals contributing the energy between
two electrons i and j: Coulomb integrals Jij and
exchange integral Kij;
2. For two electrons with different spins, there is only
Coulomb integral Jij;
3. For two electrons with the same spins, both
Coulomb and exchange integrals exist.
4. Total Hartree-Fock energy consists of the
contributions from one-electron integrals fii and
two-electron Coulomb integrals Jij and exchange
integrals Kij;
5. At the Hartree-Fock Level there are two possible
Coulomb potentials (or operators) between two
electrons i and j: Coulomb operator and exchange
operator; Jj(i) is the Coulomb potential (operator)
that i feels from j, and Kj(i) is the exchange
potential (operator) that that i feels from j.
6. Fock operator (or, average Hamiltonian) consists
of one-electron operators f(i) and Coulomb
operators Jj(i) and exchange operators Kj(i)













N electrons spin up and N electrons spin down.
Fock matrix for an electron 1 with spin up:
F(1) = f (1) + j [ Jj(1) - Kj(1) ] + j Jj(1)
j=1,N
j=1,N
Fock matrix for an electron 1 with spin down:
F(1) = f (1) + j [ Jj(1) - Kj(1) ] + j Jj(1)
j=1,N
j=1,N
f(1)  -(h2/2me)12 - N ZN/r1N
Jj(1)   dr2 fj*(2) e2/r12 fj(2)
Kj(1) q(1)  fj(1)  dr2 fj*(2) e2/r12 q(2)
Energy = j fjj+j fjj+(1/2) i j ( Jij - Kij )
+ (1/2) i j ( Jij - Kij ) + i j Jij
+ VNN
i=1,N j=1,N
fjj  fjj  <fj| f |fj>
Jij  Jij  <fj(2)| Ji(1) |fj(2)>
Kij  Kij  <fj(2)| Ki(1) |fj(2)>
Jij  Jij  <fj(2)| Ji(1) |fj(2)>
Close subshell case: ( N= N= n/2 )
F(1) = f (1) + j=1,n/2 [ 2Jj(1) - Kj(1) ]
Energy = 2 j=1,n/2 fjj + i=1,n/2 j=1,n/2 ( 2Jij - Kij ) +VNN
Hartree-Fock Method
1. Many-Body Wave Function is approximated
by Slater Determinant
2. Hartree-Fock Equation
F fi = ei fi
F Fock operator
fi the i-th Hartree-Fock orbital
ei the energy of the i-th Hartree-Fock orbital
3. Roothaan Method (introduction of Basis functions)
fi = k cki yk LCAO-MO
{ yk } is a set of atomic orbitals (or basis functions)
4. Hartree-Fock-Roothaan equation
j ( Fij - ei Sij ) cji = 0
Fij  < yi| F | yj >
Sij  < yi| yj >
5. Solve the Hartree-Fock-Roothaan equation
self-consistently
The Condon-Slater Rules
<fa(1)fb(2)fc(3)...fd(n) | f(1) | fe(1)ff(2)fg(3)...fh(n)>
= <fa(1) | f(1) | fe(1)> < fb(2)fc(3)...fd(n) | ff(2)fg(3)...fh(n)>
= <fa(1) | f(1) | fe(1)>
if b=f, c=g, ..., d=h;
0, otherwise
<fa(1)fb(2)fc(3)...fd(n) | V12 | fe(1)ff(2)fg(3)...fh(n)>
= <fa(1) fb(2) | V12 | fe(1) ff(2)> < fc(3)...fd(n) | fg(3)...fh(n)>
= <fa(1) fb(2) | V12 | fe(1) ff(2)>
if c=g, ..., d=h;
0, otherwise
LUMO
------the lowest unoccupied molecular orbital  -------
HOMO
the highest occupied molecular orbital  ------------Koopman’s Theorem
The energy required to remove an electron from a
closed-shell atom or molecules is well approximated
by minus the orbital energy e of the AO or MO from
which the electron is removed.
# HF/6-31G(d)
Route section
water energy
Title
0
O
H
H
Molecule Specification
(in Cartesian coordinates
1
-0.464 0.177 0.0
-0.464 1.137 0.0
0.441 -0.143 0.0
Basis Set fi = p cip

p
Slater-type
orbitals (STO)
nlm = N rn-1exp(-r/a0) Ylm(q,f)
 the orbital exponent
* is used instead of y in the textbook
Gaussian type functions
gijk = N xi yj zk exp(-r2)
(primitive Gaussian function)
p = u dup gu
(contracted Gaussian-type function, CGTF)
u = {ijk}
p = {nlm}
Basis set of GTFs
STO-3G, 3-21G, 4-31G, 6-31G, 6-31G*, 6-31G**
-------------------------------------------------------------------------------------
complexity & accuracy
Minimal basis set: one STO for each atomic orbital (AO)
STO-3G: 3 GTFs for each atomic orbital
3-21G: 3 GTFs for each inner shell AO
2 CGTFs (w/ 2 & 1 GTFs) for each valence AO
6-31G: 6 GTFs for each inner shell AO
2 CGTFs (w/ 3 & 1 GTFs) for each valence AO
6-31G*: adds a set of d orbitals to atoms in 2nd & 3rd rows
6-31G**: adds a set of d orbitals to atoms in 2nd & 3rd rows
Polarization
and a set of p functions to hydrogen
Function
Diffuse Basis Sets:
For excited states and in anions where electronic density
is more spread out, additional basis functions are needed.
Diffuse functions to 6-31G basis set as follows:
6-31G* - adds a set of diffuse s & p orbitals to atoms
in 1st & 2nd rows (Li - Cl).
6-31G** - adds a set of diffuse s and p orbitals to atoms
in 1st & 2nd rows (Li- Cl) and a set of diffuse
s functions to H
Diffuse functions + polarisation functions:
6-31+G*, 6-31++G*, 6-31+G** and 6-31++G** basis sets.
Double-zeta (DZ) basis set:
two STO for each AO
(10s4p)  [3s2p]
6-31G for a carbon atom:
1s
2s
2pi (i=x,y,z)
6GTFs
3GTFs
1GTF
3GTFs
1GTF
1CGTF
(s)
1CGTF
(s)
1CGTF
(s)
1CGTF
(p)
1CGTF
(p)
Minimal basis set:
One STO for each inner-shell and
valence-shell AO of each atom
example: C2H2 (2S1P/1S)
C: 1S, 2S, 2Px,2Py,2Pz
H: 1S
total 12 STOs as Basis set
Double-Zeta (DZ) basis set:
two STOs for each and
valence-shell AO of each atom
example: C2H2 (4S2P/2S)
C: two 1S, two 2S,
two 2Px, two 2Py,two 2Pz
H: two 1S (STOs)
total 24 STOs as Basis set
Split -Valence (SV) basis set
Two STOs for each inner-shell and valence-shell AO
One STO for each inner-shell AO
Double-zeta plus polarization set(DZ+P, or DZP)
Additional STO w/l quantum number larger
than the lmax of the valence - shell
 ( 2Px, 2Py ,2Pz ) to H
 Five 3d Aos to Li - Ne , Na -Ar
 C2H5 O Si H3 :
(6s4p1d/4s2p1d/2s1p)
Si
C,O
H
Assignment: Calculate the structure, ground
state energy, molecular orbital energies, and
vibrational modes and frequencies of a water
molecule using Hartree-Fock method with
3-21G basis set. (due 30/10)
Ab Initio Molecular Orbital Calculation: H2O
(using HyperChem)
1. L-Click on (click on left button of Mouse) “Startup”, and select and
L-Click on “Program/Hyperchem”.
2. Select “Build’’ and turn on “Explicit Hydrogens”.
3. Select “Display” and make sure that “Show Hydrogens” is on; L-Click
on “Rendering” and double L-Click “Spheres”.
4. Double L-Click on “Draw” tool box and double L-Click on “O”.
5. Move the cursor to the workspace, and L-Click & release.
6. L-Click on “Magnify/Shrink” tool box, move the cursor to the
workspace; L-press and move the cursor inward to reduce the size of
oxygen atom.
7. Double L-Click on “Draw” tool box, and double L-Click on “H”; Move
the cursor close to oxygen atom and L-Click & release. A hydrogen
atom appears. Draw second hydrogen atom using the same procedure.
8. L-Click on “Setup” & select “Ab Initio”; double L-Click on 3-21G;
then L-Click on “Option”, select “UHF”, and set “Charge” to 0 and
“Multiplicity” to 1.
9. L-Click “Compute”, and select “Geometry Optimization”, and L-Click
on “OK”; repeat the step till “Conv=YES” appears in the bottom bar.
Record the energy.
10.L-Click “Compute” and L-Click “Orbitals”; select a energy level,
record the energy of each molecular orbitals (MO), and L-Click “OK”
to observe the contour plots of the orbitals.
11.L-Click “Compute” and select “Vibrations”.
12.Make sure that “Rendering/Sphere” is on; L-Click “Compute” and
select “Vibrational Spectrum”. Note that frequencies of different
vibrational modes.
13.Turn on “Animate vibrations”, select one of the three modes, and
L-Click “OK”. Water molecule begins to vibrate. To suspend the
animation, L-Click on “Cancel”.
The Hartree-Fock treatment of H2
e+
e+
The Valence-Bond Treatment of H2
f1 = y1(1) y2(2)
f2 = y1(2) y2(1)
f = c1 f1 + c2 f2
H11 - W
H21 - S W
H12 - S W
H22 - W
=0
H11 = H22 = <y1(1) y2(2)|H|y1(1) y2(2)>
H12 = H21 = <y1(1) y2(2)|H|y1(2) y2(1)>
S = <y1(1) y2(2)|y1(2) y2(1)> [ = S2 ]
The Heitler-London ground-state wave function
{[y1(1) y2(2) + y1(2) y2(1)]/2(1+S)1/2} [(1)(2)-(2)(1)]/2
Comparison of the HF and VB Treatments
HF LCAO-MO wave function for H2
[y1(1) + y2(1)] [y1(2) + y2(2)]
= y1(1) y1(2) + y1(1) y2(2) + y2(1) y1(2) + y2(1) y2(2)
H- H+
H
H
VB wave function for H2
y1(1) y2(2) + y2(1) y1(2)
H
H
H
H
H
H
H+ H-
At large distance, the system becomes
H ............ H
MO: 50% H ............ H
50% H+............ HVB: 100% H ............ H
The VB is computationally expensive and requires
chemical intuition in implementation.
The Generalized valence-bond (GVB) method is a
variational method, and thus computationally feasible.
(William A. Goddard III)
f1 = y (1)y (2)
f 2 = y (2)y (1)
 = c1 f1 + c2 f 2
R
R
H 11 = H 22 = y 1 (1)y 2 (2) H y 1 (1)y 2 (2)
H 11 - W H 12 - SW
H 21 - SW H 22 - W
=0
H 12 = H 21 = y 1 (1)y 2 (2) H y 1 (2)y 2 (1)
[
S = y 1 (1)y 2 (2) y 1 (2)y 2 (1) = S 2
]
The Heitler-London ground-state wave function
[y (1)y
1
2

(2) + y 1 (2)y 2 (1)]/ 2(1 + S ) [ (1)  (2) -  (2)  (1)]/ 2
Assignment
8.4, 8.10, 8.12b, 8.40, 10.5, 10.6, 10.7, 10.8,
11.37, 13.37
Electron Correlation
Human Repulsive Correlation
Electron Correlation:
avoiding each other
Two reasons of the instantaneous correlation:
(1) Pauli Exclusion Principle (HF includes the effect)
(2) Coulomb repulsion (not included in the HF)
Beyond the Hartree-Fock
Configuration Interaction (CI)*
Perturbation theory*
Coupled Cluster Method
Density functional theory
-e
-e
r12
r2
r1
+2e
H = - (h2/2me)12 - 2e2/r1 - (h2/2me)22 - 2e2/r2 + e2/r12
H10
H20
H’
H 0 = H 10 + H 2 0
y(0)(1,2) = F1(1) F2(2)
H10 F1(1) = E1 F1(1)
H20 F2(1) = E2 F2(1)
E1 = -2e2/n12a0 n1 = 1, 2, 3, ...
E2 = -2e2/n22a0 n2 = 1, 2, 3, ...
Ground state wave function
y(0)(1,2) = (1/1/2)(2/a0)3/2exp(-2r1/a0) * (1/1/2)(2/a0)3/2exp(-2r1/a0)
E(0) = - 4e2/a0
E(1) = <y(0)(1,2)| H’ |y(0)(1,2)> = 5e2/4a0
E  E(0) + E(1) = -108.8 + 34.0 = -74.8 (eV)
[compared with exp. -79.0 eV]
Nondegenerate Perturbation Theory
(for Non-Degenerate Energy Levels)
H = H0 + H’
H0yn(0) = En(0) yn(0)
yn(0) is an eigenstate for unperturbed
system
H’ is small compared with H0
Introducing a parameter l
H(l) = H0 + lH’
H(l) yn(l) = En(l) yn(l)
yn(l) = yn(0) + l yn(1) + l2 yn(2) + ... + lk yn(k) + ...
En(l) = En(0) + l En(1) + l2 En(2) + ... + lk En(k) + ...
l = 1, the original Hamiltonian
yn = yn(0) + yn(1) + yn(2) + ... + yn(k) + ...
En = En(0) + En(1) + En(2) + ... + En(k) + ...
Where, < yn(0) | yn(j) > = 0, j=1,2,...,k,...
H0yn(0) = En(0) yn(0)
solving for En(0), yn(0)
H0yn(1) + H’ yn(0) = En(0) yn(1) + En(1)yn(0)
solving for En(1), yn(1)
H0yn(2) + H’ yn(1) = En(0) yn(2) + En(1)yn(1) + En(2)yn(0)
 solving for En(2),yn(2)
The first order:
Multiplied ym(0) from the left and integrate,
<ym(0) | H0 | yn(1) > + < ym(0) | H' | yn(0) > = < ym(0)|yn(1) >En(0) + En(1) dmn
<ym(0)|yn(1) > [Em(0)- En(0)] + < ym(0) | H' | yn(0) > = En(1) dmn
For m = n,
En(1) = < yn(0) | H' | yn(0) > Eq.(1)
For m  n, <ym(0)|yn(1) > = < ym(0) | H' | yn(0) > / [En(0)- Em(0)]
If we expand yn(1) =  cnm ym(0),
cnm = < ym(0) | H' | yn(0) > / [En(0)- Em(0)] for m  n;
(1)
(0) | H' | y (0) > / [E (0)- E (0)] y (0)
cynn
=
0.
=

<
y
n
m
m
n
n
m
m
Eq.(2)
The second order:
<ym(0)|H0|yn(2) > + < ym(0)|H'|yn(1) > = < ym(0)|yn(2)
>En(0) + < ym(0)|yn(1) >En(1) + En(1) dmn
Set m = n, we have
En(2) = m  n |<ym(0) | H' | yn(0) >|2 / [En(0)- Em(0)]
Eq.(3)
Discussion: (Text Book: page 522-527)
a. Eq.(2) shows that the effect of the perturbation
on the wave function yn(0) is to mix in
contributions from the other zero-th order states
ym(0) mn. Because of the factor 1/(En(0)-Em(0)),
the most important contributions to the yn(1)
come from the states nearest in energy to state n.
b. To evaluate the first-order correction in energy,
we need only to evaluate a single integral H’nn;
to evaluate the second-order energy correction,
we must evalute the matrix elements H’ between
the n-th and all other states m.
c. The summation in Eq.(2), (3) is over all the states,
not the energy levels.
Moller-Plesset (MP) Perturbation Theory
The MP unperturbed Hamiltonian H0
H0 = m F(m)
where F(m) is the Fock operator for electron m.
And thus, the perturbation H’
H’ = H - H0
Therefore, the unperturbed wave function is
simply the Hartree-Fock wave function .
Ab initio methods: MP2, MP4
Example One:
Consider the one-particle, one-dimensional system
with potential-energy function
V=b
for L/4 < x < 3L/4,
V=0
for 0 < x  L/4 & 3L/4  x < L
and V =  elsewhere. Assume that the magnitude
of b is small, and can be treated as a perturbation.
Find the first-order energy correction for the ground
and first excited states. The unperturbed wave
functions of the ground and first excited states are
y1 = (2/L)1/2 sin(x/L) and y2 = (2/L)1/2 sin(2x/L),
respectively.
Example Two:
As the first step of the Moller-Plesset perturbation
theory, Hartree-Fock method gives the zeroth-order
energy. Is the above statement correct?
Example Three:
Show that, for any perturbation H’, E1(0) + E1(1)  E1
where E1(0) and E1(1) are the zero-th order energy
and the first order energy correction, and E1 is the
ground state energy of the full Hamiltonian H0 + H’.
Example Four:
Calculate the bond orders of Li2 and Li2+.
Perturbation Theory for a Degenerate Energy Level
Hydrogen Atom
n=3
n=2
n=1
3s, 3px , 3py , 3pz , 3d1-5
2s, 2px , 2py , 2pz
1s
B/e
H = H0 + H’
H0yn(0) = Ed(0) yn(0)
H’ is small compared with H0
n=1,2,...,d
(1)Apply the results of nondegenerate perturbation
theory
cnm = < ym(0) | H' | yn(0) > / [En(0)- Em(0)]   for 1  m, n  d
WRONG ! something very different !
(2) What happened ?
c1 y1(0) + c2 y2(0) + ... + cd yd(0) is an eigenstate for H0
There are infinite number of such states that are
degenerate.
When H’is switched on, these states are no longer
degenerate, and nondegenerate eigenstates of
H0 + H’ appear !
Therefore, even for zero-th order of eigenstates,
there are sudden changes !
(3) Introducing a parameter l
H(l) = H0 + lH’
H(l) yn(l) = En(l) yn(l)
l = 1, the original Hamiltonian
yn(l) = fn(0) + l yn(1) + l2 yn(2) + ... + lk yn(k) + ...
En(l) = Ed(0) + l En(1) + l2 En(2) + ... + lk En(k) + ...
fn(0) = k ck yk(0)
H0yn(1) + H’ fn(0) = Ed(0) yn(1) + En(1)fn(0)
solving for En(1), fn(0) , yn(1)
Multiplied ym(0) from the left and integrate,
<ym(0) | H0 | yn(1) > + < ym(0) | H' | fn(0) >
= < ym(0)|yn(1) >Ed(0) + En(1) <ym(0)| fn(0) >
<ym(0)|yn(1) > [Em(0)- Ed(0)] + < ym(0) | H' | fn(0) >
= En(1) < ym(0)| fn(0) >
For 1  m  d,
n [< ym(0) | H' | yn(0) > - Em(1)dmn] cn = 0
Em(1) = < fm(0) | H' | fm(0) >
Assignment 2: 9.2, 9.4a, 9.9, 9.18, 9.24
Configuration Interaction (CI)
+
+ …
Single Electron Excitation or Singly Excited
Double Electrons Excitation or Doubly Excited
Singly Excited Configuration Interaction (CIS):
Changes only the excited states
+
Doubly Excited CI (CID):
Changes ground & excited states
+
Singly & Doubly Excited CI (CISD):
Most Used CI Method
Full CI (FCI):
Changes ground & excited states
+
+ ...
+
Coupled-Cluster Method
y = eT y(0)
y(0): Hartree-Fock ground state wave function
y: Ground state wave function
T = T1 + T2 + T3 + T4 + T5 + …
Tn : n electron excitation operator
T1
=
Coupled-Cluster Doubles (CCD) Method
yCCD = eT y(0)
2
y(0): Hartree-Fock ground state wave function
yCCD: Ground state wave function
T2 : two electron excitation operator
T2
=
Complete Active Space SCF (CASSCF)
Active space
All possible configurations
Density-Functional Theory (DFT)
Hohenberg-Kohn Theorem: Phys. Rev. 136, B864 (1964)
The ground state electronic density (r) determines
uniquely all possible properties of an electronic system
(r)  Properties P (e.g. conductance), i.e. P  P[(r)]
Density-Functional Theory (DFT)
E0 = - (h2/2me)i <yi |i2 |yi >-   dr Ze2(r) / r1
+ (1/2)   dr1 dr2 e2/r12 + Exc[(r)]
Kohn-Sham Equation Ground State: Phys. Rev. 140, A1133 (1965)
FKS yi = ei yi
FKS  - (h2/2me)ii2 -  Ze2 / r1 + j Jj + Vxc
Vxc  dExc[(r)] / d(r)
A popular exchange-correlation functional Exc[(r)]: B3LYP
180 small- or medium-size organic molecules:
1. C.L. Yaws, Chemical Properties Handbook,
(McGraw-Hill, New York, 1999)
2. D.R. Lide, CRC Handbook of Chemistry and Physics,
3rd ed. (CRC Press, Boca Raton, FL, 2000)
3. J.B . Pedley, R.D. Naylor, S.P. Kirby,
Thermochemical data of organic compunds,
2nd ed. (Chapman and Hall, New York, 1986)
Differences of heat of formation in three references
for same compound are less than 1 kcal/mol; and
error bars are all less than 1kcal/mol
Hu, Wang, Wong & Chen, J. Chem. Phys. (Comm) (2003)
B3LYP/6-311+G(d,p)
B3LYP/6-311+G(3df,2p)
RMS=21.4 kcal/mol
RMS=12.0 kcal/mol
RMS=3.1 kcal/mol
RMS=3.3 kcal/mol
B3LYP/6-311+G(d,p)-NEURON & B3LYP/6-311+G(d,p)-NEURON: same
accuracy
Time-Dependent Density-Functional Theory (TDDFT)
Runge-Gross Extension: Phys. Rev. Lett. 52, 997 (1984)
Time-dependent system
(r,t)  Properties P (e.g. absorption)
TDDFT equation: exact for excited states
Isolated system
Open system
Density-Functional Theory for Open System ?
 Holographic electron density theorem for
time-independent systems
•
Riess and Munch (1981)
•
Mezey (1999)
•
Fournais (2002)
(r)
D
Analytical
continuation
D(r)
(r)
HK
system properties
 Holographic electron density theorem for
time-dependent systems
It is difficult to prove the analyticity for (r,t) rigorously!
D(r,t)
v(r,t)
system properties
Holographic electron
density theorem
(r,t)
D
Zheng, Wang, Yam, Mo & Chen, PRB (2007)
The electron density distribution of the
reduced system determines all physical
properties or processes of the entire system!
Existence of a rigorous TDDFT for Open System
A Grain of Sand
William Blake
To see a world in a grain of sand,
And a heaven in a wild flower,
Hold infinity in the palm of your hand,
And eternity in an hour.
一粒沙子
威廉.布莱克
从一粒沙子看到一个世界，
从一朵野花看到一个天堂，
把握在你手心里的就是无限，
永恒也就消融于一个时辰。
Transient current (red lines) & applied
bias voltage (green lines) for the AlCNT-Al system. (a) Bias voltage is
turned on exponentially, Vb = V0 (1-et/a) with V = 0.1 mV & a = 1 fs. Blue
0
line in (a) is a fit to transient current,
I0(1-e-t/τ) with τ = 2.8 fs & I0 =13.9 nA.
(b) Bias voltage is sinusoidal with a
period of T = 5 fs. The red line is for
the current from the right electrode &
squares are the current from the left
electrode at different times.
Vb = V0 (1-e-t/a)
V0 = 0.1 mV & a = 1 fs
Switch-on time: ~ 10 fs
V (t) = V (1- e
(-t / a )

)
Ground State Excited State CPU Time Correlation Geometry Size Consistent
(CHNH,6-31G*)
HFSCF


1
0
OK

DFT


~1


CIS


<10
OK

CISD


17


CISDTQ




MP2


1.5
MP4


CCD

CCSDT

80-90%
(20 electrons)
very large 98-99%


5.8
85-95%
(DZ+P)
>90%



large
>90%



very large
~100%


Relativistic Effects
Speed of 1s electron: Zc / 137
Heavy elements have large Z, thus relativistic effects are
important.
Dirac Equation:
Relativistic Hartree-Fock w/ Dirac-Fock operator; or
Relativistic Kohn-Sham calculation; or
Relativistic effective core potential (ECP).
Four Sources of error in ab initio Calculation
(1) Neglect or incomplete treatment of electron correlation
(2) Incompleteness of the Basis set
(3) Relativistic effects
(4) Deviation from the Born-Oppenheimer approximation
Semiempirical Molecular Orbital Calculation
Extended Huckel MO Method
(Wolfsberg, Helmholz, Hoffman)
Independent electron approximation
Hval = i Heff(i)
Heff(i) = -(h2/2m) i2 + Veff(i)
Schrodinger equation for electron i
Heff(i) fi = ei fi
LCAO-MO:
fi = r cri yr
s ( Heffrs - ei Srs ) csi = 0
Heffrs  < yr| Heff | ys >
Srs  < yr| ys >
Parametrization:
Heffrr  < yr| Heff | yr >
= minus the valence-state ionization
potential (VISP)
-----------------------------------------------------------------------
Atomic Orbital Energy
e5
e4
e3
e2
e1
Heffrs = ½ K (Heffrr + Heffss) Srs
VISP
-e5
-e4
-e3
-e2
-e1
K:
13
CNDO, INDO, NDDO
(Pople and co-workers)
Hamiltonian with effective potentials
Hval = i [ -(h2/2m) i2 + Veff(i) ] + ij>i e2 / rij
two-electron integral:
(rs|tu) = <yr(1) yt(2)| 1/r12 | ys(1) yu(2)>
CNDO: complete neglect of differential overlap
(rs|tu) = drs dtu (rr|tt)  drs dtu rt
INDO: intermediate neglect of differential overlap
(rs|tu) = drs dtu (rr|tt) when r, s, t & u not on same atom;
(rs|tu)  0 when r, s, t and u are on the same atom.
NDDO: neglect of diatomic differential overlap
(rs|tu) = 0 if r and s (or t and u) are not on the
same atom.
CNDO, INDO are parametrized so that the overall
results fit well with the results of minimal basis ab
initio Hartree-Fock calculation.
CNDO/S, INDO/S are parametrized to predict
optical spectra.
MINDO, MNDO, AM1, PM3
(Dewar and co-workers, University of Texas,
Austin)
MINDO: modified INDO
MNDO: modified neglect of diatomic overlap
AM1: Austin Model 1
PM3: MNDO parametric method 3
MINDO, MNDO, AM1 & PM3:
*based on INDO & NDDO
*reproduce the binding energy
Molecular Mechanics (MM) Method
F = Ma
F : Force Field
Molecular Mechanics Force Field
•
•
•
•
Bond Stretching Term
Bond Angle Term
Torsional Term
Non-Bonding Terms: Electrostatic
Interaction & van der Waals Interaction
Bond Stretching Potential
Eb = 1/2 kb (l)2
where, kb : stretch force constant
l : difference between equilibrium
& actual bond length
Two-body interaction
Bond Angle Deformation Potential
Ea = 1/2 ka ()2
where, ka : angle force constant
 : difference between equilibrium
& actual bond angle

Three-body interaction
Periodic Torsional Barrier Potential
Et = (V/2) (1+ cosn )
where, V : rotational barrier
 : torsion angle
n : rotational degeneracy
Four-body interaction
C2H3Cl
Non-bonding interaction
van der Waals interaction
for pairs of non-bonded atoms
Coulomb potential
for all pairs of charged atoms
MM Force Field Types
•
•
•
•
•
MM2
AMBER
CHAMM
BIO
OPLS
Small molecules
Polymers
Polymers
Polymers
Solvent Effects
CHAMM FORCE FIELD FILE
########################################################
##
##
## TINKER Atom Class Numbers to CHARMM22 Atom Names ##
##
##
##
1 HA
11 CA
21 CY
31 NR3
##
##
2 HP
12 CC
22 CPT
32 NY
##
##
3 H
13 CT1
23 CT
33 NC2
##
##
4 HB
14 CT2
24 NH1
34 O
##
##
5 HC
15 CT3
25 NH2
35 OH1
##
##
6 HR1
16 CP1
26 NH3
36 OC
##
##
7 HR2
17 CP2
27 N
37 S
##
##
8 HR3
18 CP3
28 NP
38 SM
##
##
9 HS
19 CH1
29 NR1
##
##
10 C
20 CH2
30 NR2
##
##
##
########################################################
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
atom
1
1
1
HA
"Nonpolar Hydrogen"
1
1.008
2
2
HP
"Aromatic Hydrogen"
1
1.008
3
3
H
"Peptide Amide HN"
1
1.008
4
4
HB
"Peptide HCA"
1
1.008
5
4
HB
"N-Terminal HCA"
1
1.008
6
5
HC
"N-Terminal Hydrogen"
1
1.008
7
5
HC
"N-Terminal PRO HN"
1
1.008
8
3
H
"Hydroxyl Hydrogen"
1
1.008
9
3
H
"TRP Indole HE1"
1
1.008
10
3
H
"HIS+ Ring NH"
1
1.008
11
3
H
"HISDE Ring NH"
1
1.008
12
6
HR1
"HIS+ HD2/HISDE HE1"
1
1.008
vdw
vdw
vdw
vdw
vdw
vdw
vdw
vdw
vdw
vdw
################################
##
##
## Van der Waals Parameters ##
##
##
################################
/Ao
1
2
3
4
5
6
7
8
9
10
1.3200
1.3582
0.2245
1.3200
0.2245
0.9000
0.7000
1.4680
0.4500
2.0000
/(kcal/mol)
-0.0220
-0.0300
-0.0460
-0.0220
-0.0460
-0.0460
-0.0460
-0.0078
-0.1000
-0.1100
bond
bond
bond
bond
bond
bond
bond
bond
bond
##################################
##
##
## Bond Stretching Parameters ##
##
##
##################################
/(kcal/mol/Ao2)
1
1
1
1
1
1
1
1
1
10
11
12
13
14
15
17
18
21
330.00
340.00
317.13
309.00
309.00
322.00
309.00
309.00
330.00
/Ao
1.1000
1.0830
1.1000
1.1110
1.1110
1.1110
1.1110
1.1110
1.0800
angle
angle
angle
angle
angle
angle
angle
angle
angle
angle
angle
angle
################################
##
##
## Angle Bending Parameters ##
##
##
################################
/(kcal/mol/rad2)
3
13
13
13
14
14
14
15
15
15
16
16
10
10
10
10
10
10
10
10
10
10
10
10
34
24
27
34
24
27
34
24
27
34
24
27
50.00
80.00
20.00
80.00
80.00
20.00
80.00
80.00
20.00
80.00
80.00
20.00
/deg
121.70
116.50
112.50
121.00
116.50
112.50
121.00
116.50
112.50
121.00
116.50
112.50
############################
##
##
## Torsional Parameters ##
##
##
############################
torsion
1
11
11
1
torsion
1
11
11
11
torsion
1
11
11
22
torsion
2
11
11
2
torsion
2
11
11
11
torsion
2
11
11
14
torsion
2
11
11
15
torsion
2
11
11
22
torsion
2
11
11
35
torsion
2
11
11
36
torsion
11
11
11
11
torsion
11
11
11
14
torsion
11
11
11
15
torsion
11
11
11
22
torsion
11
11
11
35
torsion
11
11
11
36
/(kcal/mol)
2.500
3.500
3.500
2.400
4.200
4.200
4.200
3.000
4.200
4.200
3.100
3.100
3.100
3.100
3.100
3.100
/deg
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
180.0
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
2
Algorithms for Molecular Dynamics
Euler method:
x(t+t) = x(t) + (dx/dt) t
Fourth-order Runge-Kutta method:
x(t+t) = x(t) + (1/6) (s1+2s2+2s3+s4) t +O(t5)
s1 = dx/dt
s2 = dx/dt
[w/ t=t+t/2, x = x(t)+s1t/2]
s3 = dx/dt
[w/ t=t+t/2, x = x(t)+s2t/2]
s4 = dx/dt
[w/ t=t+t, x = x(t)+s3 t]
Very accurate but slow!
Algorithms for Molecular Dynamics
Verlet Algorithm:
x(t+t) = x(t) + (dx/dt) t + (1/2) d2x/dt2 t2 + ...
x(t -t) = x(t) - (dx/dt) t + (1/2) d2x/dt2 t2 - ...
x(t+t) = 2x(t) - x(t -t) + d2x/dt2 t2 + O(t4)
Efficient & Commonly Used!
Calculated Properties
•
•
•
•
Structure, Geometry
Energy & Stability
Vibration Frequency & Mode
Real Time Dynamics
Dynamics simulation for Protein folding
Simulation Time: 100ps
Temperature: 300K
RIBOSOMAL PROTEIN (C-TERMINAL DOMAIN)
PDB code: 1CTF (68 amino acid)
Summary
Hamiltonian
H = (-h2/2m)2 - (h2/2me)ii2 +  ZZe2/r
- i  Ze2/ri + i j e2/rij
The variation theorem
Consider a system whose Hamiltonian operator H is time independent
and whose lowest-energy eigenvalue is E1. If f is any well-behaved
function that satisfies the boundary conditions of the problem, then
 f* H f d /  f* f d > E1
Variational Method
(1) Construct a wave function f(c1,c2,,cm)
(2) Calculate the energy of f:
Ef  Ef(c1,c2,,cm)
(3) Choose {cj*} (i=1,2,,m) so that Ef is minimum
Extension of Variation Method
For a wave function f which is orthogonal to the ground state
wave function y1, i.e.
d f*y1 = 0
Ef = d f*Hf / d f*f > E2
the first excited state energy
The Pauli principle
two electrons cannot be in the same state
the wave function of a system of electrons must be antisymmetric with
respect to interchanging of any two electrons.
Slater determinant of H2 :
y(1,2) = 1/2! [f(1)(1)f(2)(2) - f(2)(2)f(1)(1)]
f(1)(1) f(2)(2)
= 1/2!
f(1)(1) f(2)(2)
Hartree-Fock Equation:
[ f(1)+ J2(1) - K2(1)] f1(1) = e1 f1(1)
[ f(2)+ J1(2) - K1(2)] f2(2) = e2 f2(2)
Fock Operator:
F(1)  f(1)+ J2(1) - K2(1)
F(2)  f(2)+ J1(2) - K1(2)
LCAO-MO:
Fock operator for 1
Fock operator for 2
f = c1y1 + c2y2
Express Hartree-Fock energy
in terms of fi, Jij & Kij
Molecule
H2+
H2
He2+
He2
Li2
Be2
C2
N2+
N2
O2
Bond order
1/2
1
1/2
0
1
0
2
1/2
3
2
De/eV
2.79
4.75
1.08
0.0009
1.07
0.10
6.3
8.85
9.91
5.21
Basis set of GTFs
STO-3G, 3-21G, 4-31G, 6-31G, 6-31G*, 6-31G**
-------------------------------------------------------------------------------------
complexity & accuracy
Gaussian 98 Input file
# HF/6-31G(d)
Route section
water energy
Title
0
O
H
H
Molecule Specification
(in Cartesian coordinates
1
-0.464 0.177 0.0
-0.464 1.137 0.0
0.441 -0.143 0.0
Comparison of the HF and VB Treatments
Electron Correlation
Beyond the Hartree-Fock
Configuration Interaction (CI)*
Perturbation theory*
En(1) = < yn(0) | H' | yn(0) >
Coupled Cluster Method
Density functional theory
Moller-Plesset (MP) Perturbation Theory
The MP unperturbed Hamiltonian H0
H0 = m F(m)
where F(m) is the Fock operator for electron m.
And thus, the perturbation H’
H ’ = H - H0
Ground State Excited State CPU Time Correlation Geometry Size Consistent
(CH3NH2,6-31G*)
HFSCF


1
0
OK

DFT


~1


CIS


<10
OK

CISD


17


CISDTQ




MP2


1.5
MP4


CCD

CCSDT

80-90%
(20 electrons)
very large 98-99%


5.8
85-95%
(DZ+P)
>90%



large
>90%



very large
~100%


Molecular Mechanics Force Field
•
•
•
•
Bond Stretching Term
Bond Angle Term
Torsional Term
Non-Bonding Terms: Electrostatic
Interaction & van der Waals Interaction
F = Ma
F : Force Field