Transcript No Slide Title

Computational Modeling of
Macromolecular Systems
Dr. GuanHua CHEN
Department of Chemistry
University of Hong Kong
Computational Chemistry
• Quantum Chemistry
SchrÖdinger Equation
H = E
• Molecular Mechanics
F = Ma
F : Force Field
Computational Chemistry Industry
Company
Software
Gaussian Inc.
Gaussian 94, Gaussian 98
Schrödinger Inc.
Jaguar
Wavefunction
Spartan
Q-Chem
Q-Chem
Molecular Simulation Inc. (MSI) InsightII, Cerius2, modeler
HyperCube
HyperChem
Applications: material discovery, drug design & research
R&D in Chemical & Pharmaceutical industries in 2000: US$ 80 billion
Sales of Scientific Computing in 2000: > US$ 200 million
Cytochrome c (involved in the ATP synthesis)
heme
1997 Nobel Prize
in Biology:
ATP Synthase in
Mitochondria
Cytochrome c is a peripheral membrane protein
involved in the long distance electron transfers
Simulation of a pair of polypeptides
Duration: 100 ps. Time step: 1 ps (Ng, Yokojima & Chen, 2000)
Protein Dynamics
1. Atomic Fluctuations
10-15 to 10-11 s; 0.01 to 1 Ao
2. Collective Motions
10-12 to 10-3 s; 0.01 to >5 Ao
3. Conformational Changes
o
-9
3
10 to 10 s; 0.5 to >10 A
Theoretician leaded the way ! (Karplus at Harvard U.)
Scanning Tunneling Microscope
Manipulating Atoms by Hand
Large Gear Drives Small Gear
G. Hong et. al., 1999
Calculated Electron distribution at equator
Vitamin C
The electron density around the vitamin C molecule. The colors show the electrostatic
potential with the negative areas shaded in red and the positive in blue.
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 (Dl)2
where, kb : stretch force constant
Dl : difference between equilibrium
& actual bond length
Two-body interaction
Bond Angle Deformation Potential
Ea = 1/2 ka (D)2
where, ka : angle force constant
D : 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
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
Runge-Kutta methods:
x(t+Dt) = x(t) + (dx/dt) Dt
Fourth-order Runge-Kutta
x(t+Dt) = x(t) + (1/6) (s1+2s2+2s3+s4) Dt +O(Dt5)
s1 = dx/dt
s2 = dx/dt
[w/ t=t+Dt/2, x = x(t)+s1Dt/2]
s3 = dx/dt
[w/ t=t+Dt/2, x = x(t)+s2Dt/2]
s4 = dx/dt
[w/ t=t+Dt, x = x(t)+s3 Dt]
Very accurate but slow!
Algorithms for Molecular Dynamics
Verlet Algorithm:
x(t+Dt) = x(t) + (dx/dt) Dt + (1/2) d2x/dt2 Dt2 + ...
x(t -Dt) = x(t) - (dx/dt) Dt + (1/2) d2x/dt2 Dt2 - ...
x(t+Dt) = 2x(t) - x(t -Dt) + d2x/dt2 Dt2 + O(Dt4)
Efficient & Commonly Used!
Calculated Properties
• Structure, Geometry
• Energy & Stability
• Mechanic Properties: Young’s
Modulus
• Vibration Frequency & Mode
Crystal Structure of C60 solid
Crystal Structure of K3C60
Vibration Spectrum of K3C60
GH Chen, Ph.D. Thesis, Caltech (1992)
Quantum Chemistry Methods
• Ab initio Molecular Orbital Methods
Hartree-Fock, Configurationa Interaction (CI)
MP Perturbation, Coupled-Cluster, CASSCF
• Density Functional Theory
• Semiempirical Molecular Orbital
Methods
Huckel, PPP, CNDO, INDO, MNDO, AM1
PM3, CNDO/S, INDO/S
SchrÖdinger Equation
H=E
Wavefunction
Hamiltonian
H = (-h2/2m)2 - (h2/2me)ii2
+  ZZe2/r - i  Ze2/ri
+ i j e2/rij
Energy
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)
e-
+
e-
Fock operator for 1
Fock operator for 2
f(1)  -(h2/2me)12 - N ZN/r1N
one-electron term if no Coulomb interaction
J2(1)   dr2 f2*(2) e2/r12 f2(2)
Ave. Coulomb potential on electron 1 from 2
K2(1) f1(1)  f2(1)  dr2 f2*(2) e2/r12 f1(2)
Ave. exchange potential on electron 1 from 2
f(2)  -(h2/2me)22 - N ZN/r2N
J1(2)   dr1 f1*(1) e2/r12 f1(1)
K1(2) q(2)  f1(1)  dr1 f1*(1) e2/r12 q(1)
Average Hamiltonian for electron 1
F(1)  f(1)+ J2(1) - K2(1)
Average Hamiltonian for electron 2
F(2)  f(2)+ J1(2) - K1(2)
Hartree-Fock Method
1. Many-Body Wave Function is approximated
by Single 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 k LCAO-MO
{ k } is a set of atomic orbitals (or basis functions)
4. Hartree-Fock-Roothaan equation
j ( Fij - ei Sij ) cji = 0
Fij  < i|F | j >
Sij  < i| j >
5. Solve the Hartree-Fock-Roothaan equation
self-consistently (HFSCF)
Graphic Representation of Hartree-Fock Solution
0 eV
Ionization
Energy
Electron
Affinity
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.
Basis Set fi = p cip

Slater-type
orbitals (STO)
p
nlm = N rn-1exp(-r/a0) Ylm(q,f)
 the orbital exponent
Gaussian type functions (GTF)
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
(10s12p)  [3s6p]
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)
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
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
+
+ ...
+
Perturbation Theory
H = H0 + H’
H0n(0) = En(0) n(0)
n(0) is an eigenstate for unperturbed system
H’ is small compared with H0
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, MP3, MP4
Coupled-Cluster Method
 = eT (0)
(0): Hartree-Fock ground state wave function
: Ground state wave function
T = T1 + T2 + T3 + T4 + T5 + …
Tn : n electron excitation operator
T1
=
Coupled-Cluster Doubles (CCD) Method
CCD = eT (0)
2
(0): Hartree-Fock ground state wave function
CCD: 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:
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 <i |i2 |i >-   dr Ze2(r) / r1
+ (1/2)   dr1 dr2 e2/r12 + Exc[(r)]
Kohn-Sham Equation:
FKS i = ei i
FKS  - (h2/2me)ii2 -  Ze2 / r1 + j Jj + Vxc
Vxc  dExc[(r)] / d(r)
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 r
s ( Heffrs - ei Srs ) csi = 0
Heffrs  < r| Heff | s >
Srs  < r| s >
Parametrization:
Heffrr  < r| Heff | r >
= 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) = <r(1) t(2)| 1/r12 | s(1) u(2)>
CNDO: complete neglect of differential overlap
(rs|tu) = drs dtu (rr|tt)  drs dtu rt
INDO: intermediate neglect of differential overlap
(rs|tu) = 0 when r, s, t and u are not 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
*based on INDO & NDDO
*reproduce the binding energy
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).
Ground State: ab initio Hartree-Fock calculation
Computational Time: protein w/ 10,000 atoms
ab initio Hartree-Fock ground state calculation:
~20,000 years on CRAY YMP
In 2010:
~24 months on 100 processor machine
One Problem:
Transitor with a few atoms
Current Computer Technology will fail !
Quantum Chemist’s Solution
Linear-Scaling Method: O(N)
Computational time scales linearly with system size
Time
Size
Linear Scaling Calculation for Ground State
Divide-and-Conqure (DAC)
W. Yang, Phys. Rev. Lett. 1991
Density-Matrix Minimization (DMM)
Method
Minimize the Energy or the Grand Potential:
 = Tr [ (32 - 23) (H-I) ]
Li, Nunes and Vanderbilt, Phy. Rev. B. 1993
Orbital Minimization (OM) Method
Minimize the Energy or the Grand Potential:
 = 2 nij cni (H-I)ij cnj
- nmij cni (H-I)ij cmj l cnl cml
Mauri (1993), Ordejon (1993), Galii (1994), Kim (1995)
Fermi Operator Expansion (FOE) Method
Expand Density Matrix in Chebyshev Polynomial:
(H) = c0I + c1H + c2H2 + …
= c0I / 2 + cjTj(H) + …
T0(H) = I
T1(H) = H
Tj+1 (H) = 2HTj(H) - Tj-1(H)
Goedecker & Colombo (1994)
York, Lee & Yang, JACS, 1996
Superoxide Dismutase (4380 atoms) AM1
Linear Scaling First Principle Method
Two-electron integrals :
Vabcd = <fa(1) fb(2) | e2 / r12 | fd(1) fc(2)>
Coulomb Integrals:
Fast Multiple Method (FMM)
Exchange-Correlation (XC):
Use of Locality
Strain, Scuseria & Frisch, Science (1996):
LSDA / 3-21G DFT calculation on 1026 atom
RNA Fragment
Linear Scaling Calculation for Ground State
Yang, Phys. Rev. Lett. 1991
Li, Nunes & Vanderbilt, Phy. Rev. B. 1993
Baroni & Giannozzi, Europhys. Lett. 1992.
Gibson, Haydock & LaFemina, Phys. Rev. B 1993.
Aoki, Phys. Rev. Lett. 1993.
Cortona, Phys. Rev. B 1991.
Galli & Parrinello, Phys. Rev. Lett. 1992.
Mauri, Galli & Car, Phys. Rev. B 1993.
Ordejón et. al., Phys. Rev. B 1993.
Drabold & Sankey, Phys. Rev. Lett. 1993.
Linear Scaling Calculation for
EXCITED STATE ?
A Much More Difficult Problem !
Linear-Scaling Calculation for excited states
Localized-Density-Matrix (LDM) Method
(0)

=
+ d
E(t)
(0)
ij = 0
rij > r0
dij = 0
rij > r1
Yokojima & Chen, Phys. Rev. B, 1999
Principle of the nearsightedness of
equilibrium systems (Kohn, 1996)
Heisenberg Equation of Motion

i  = [H ,  ]
Time-Dependent Hartree-Fock
Random Phase Approximation
CH
Polyacetylene
4
2
8
6
10
CH
N
2
N
N-2
12
...
1
3
5
7
9
11
N-3
PPP Semiempirical Hamitonian
Hˆ = Hˆ Huckel + Hˆ c + Hˆ ext
N-1
Liang, Yokojima & Chen, JPC, 2000
3000
CPU Tim e (s)
40,000
HF
CPU Time (s)
30,000
HF
LDM
 =80
2000
c =30
LDM
1000
 =50
0
0
20,000
200
400
600
Num ber of Atom s
LDM
800
LDM
 =80
10,000
c =20
 =50
c=30
0=20
0
0
5000
10000
15000
20000
Number of Atoms
Yokojima, Zhou & Chen, Chem. Phys. Lett., 1999
Liang, Yokojima & Chen, JPC, 2000
Flat Panel Display
Cambridge Display Technology
Weight: 15 gram
Resolution: 800x236
Size:
45x37 mm
Voltage:
DC, 10V
Intensity
electron
hole
Energy
Carbon Nanotube
Surprising!
DFT: no or very small gap
Liang, Wang, Yokojima & Chen, JACS (2000)
Absorption Spectra of (9,0) SWNTs
Smallest SWNT: 0.4 nm in diameter
Wang, Tang & etc., Nature (2000)
Three possibilities:
(4,2), (3,3) & (5,0) SWNTs
Tang et. al, 2000
Absorption of SWNTs (4,2), (3,3) & (5,0)
(4,2): C332H12
(3,3): C420H12
(5,0): C330
Liang, & Chen (2001)
Quantum Mechanics / Molecular
Mechanics (QM/MM) Method
Combining quantum mechanics and
molecular mechanics methods:
QM
MM
Human Genome Project
GENOMICS
Aldose Reductase
Design of Aldose Reductase Inhibitors
Goddard, Caltech