Transcript No Slide Title

Computer Modeling
Dr. GuanHua CHEN
Department of Chemistry
University of Hong Kong
http://yangtze.hku.hk/lecture/comput06-07.ppt
Computational Chemistry
• Quantum Chemistry
SchrÖdinger Equation
H = E
• Molecular Mechanics
F = Ma
F : Force Field
• Bioinformatics
Computational Chemistry Industry
Company
Gaussian Inc.
Schrödinger Inc.
Wavefunction
Q-Chem
Accelrys
HyperCube
Informatix
Celera Genomics
Software
Gaussian 94, Gaussian 98
Jaguar
Spartan
Q-Chem
InsightII, Cerius2
HyperChem
Applications: material discovery, drug design & research
R&D in Chemical & Pharmaceutical industries in 2000: US$ 80 billion
Bioinformatics: Total Sales in 2001
Project Sales in 2006
US$ 225 million
US$ 1.7 billion
Vitamin C
C60
energy
heme
Cytochrome c
OH + D2 --> HOD + D
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
- i  Ze2/ri (+  ZZe2/r )
+ i j e2/rij
Energy
One-electron terms:
(-h2/2m)2 - (h2/2me)ii2 - i  Ze2/ri
Two-electron term:
i j e2/rij
Hartree-Fock Method
Orbitals
1. 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
2. Roothaan Method (introduction of Basis functions)
fi = k cki k LCAO-MO
{ k } is a set of atomic orbitals (or basis functions)
3. Hartree-Fock-Roothaan equation
j ( Fij - ei Sij ) cji = 0
Fij  < i|F | j >
Sij  < i| j >
4. Solve the Hartree-Fock-Roothaan equation
self-consistently (HFSCF)
Graphic Representation of Hartree-Fock Solution
0 eV
Ionization
Energy
Electron
Affinity
A Gaussian Input File for H2O
# 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
{ k } is a set of atomic orbitals (or basis functions)
STO-3G, 3-21G, 4-31G, 6-31G, 6-31G*, 6-31G**
-------------------------------------------------------------------------------------
complexity & accuracy
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}
Electron Correlation:
avoiding each other
The reason of the instantaneous correlation:
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 - H 0
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: 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 <i |i2 |i >-   dr Ze2(r) / r1
+ (1/2)   dr1 dr2 e2/r12 + Exc[(r)]
Kohn-Sham Equation Ground State: Phys. Rev. 140, A1133 (1965)
FKS i = ei i
FKS  - (h2/2me)ii2 -  Ze2 / r1 + j Jj + Vxc
Vxc  dExc[(r)] / d(r)
A popular exchange-correlation functional Exc[(r)]: B3LYP
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 ???
Further Extension: X. Zheng, F. Wang & G.H. Chen (2005)
Generalized TDDFT equation: exact for open systems
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%


Search for Transition State
Transition State:
one negative frequency
k   e-DG/RT
DG
Reactant
Product
Reaction Coordinate
Gaussian Input File for Transition State Calculation
#b3lyp/6-31G opt=qst2
test
the first is the reactant internal coordinate
0 1
O
H 1 oh1
H 1 oh1
2 ohh1
oh1
0.90
ohh1 104.5
The second is the product internal coordinate
0 1
O
H 1 oh2
H 1 oh3
oh2
0.9
oh3
10.0
ohh2 160.0
2 ohh2
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).
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
Quantum Chemistry for
Complex Systems
Quantum Mechanics / Molecular
Mechanics (QM/MM) Method
Combining quantum mechanics and
molecular mechanics methods:
QM
MM
Hamiltonian of entire system:
H = HQM +HMM +HQM/MM
Energy of entire system:
E = EQM(QM) + EMM(MM) + EQM/MM(QM/MM)
EQM/MM(QM/MM) = Eelec(QM/MM) + Evdw(MM) + EMM-bond(MM)
EQM(QM) + Eelec(QM/MM) = <| Heff |>
Heff = -1/2 ii2 + ij 1/rij - i Z/ri - i q/ri
+ i Vv-b(ri) + d Z Zd/rd +  Zq/r
QM
MM
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
York, Lee & Yang, JACS, 1996
Superoxide Dismutase (4380 atoms) AM1
Strain, Scuseria & Frisch, Science (1996):
LSDA / 3-21G DFT calculation on 1026 atom
RNA Fragment
Linear Scaling Calculation for Excited State
Liang, Yokojima & Chen, JPC, 2000
Fast Multiple Method
LDM-TDDFT: CnH2n+2
LODESTAR: Software Package for Complex Systems
Characteristics：
O(N) Divide-and-Conquer
O(N) TDHF (ab initio & semiemptical)
O(N) TDDFT
Nonlinear Optical Light Harvesting System
CNDO/S-, PM3-, AM1-, INDO/S-, & TDDFT-LDM
Photo-excitations in Light Harvesting System II
strong absorption: ~800 nm
generated by VMD
generated by VMD
Carbon Nanotube
Quantum mechanical investigation of the field emission from
the tips of carbon nanotubes
Zheng, Chen, Li, Deng & Xu, Phys. Rev. Lett. 2004
Zettl, PRL 2001
Molecular Mechanics Force Field
•
•
•
•
Bond Stretching Term
Bond Angle Term
Torsional Term
Electrostatic Term
• van der Waals interaction
Molecular Mechanics
F = Ma
F : Force Field
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
Force Field Types
•
•
•
•
•
MM2
AMBER
CHAMM
BIO
OPLS
Molecules
Polymers
Polymers
Polymers
Solvent Effects
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!
Multiple Scale Simulation
Goddard, Caltech
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
Computer-aided drug design
Chemical Synthesis
Screening using in vitro assay
Animal Tests
Clinical Trials
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
Inhibitor
Aldose Reductase
Design of Aldose Reductase Inhibitors
Descriptors:
Electron negativity
Volume
Database for Functional Groups
5.0
O
NH
5' HN
Experimental values
4.5
O
6'
NMe
X
7'
8'
N
H
O
4.0
3.5
3.0
2.5
2.5
3.0
3.5
4.0
4.5
Fig 3 QSAR OF INHIBITOR CONCENTRATION OF INHIBITING AR Log(IC50)
0.8
0.7
O
NH
Experimetal values
0.6
5' HN
0.5
O
6'
NMe
X
7'
8'
0.4
N
H
O
0.3
0.2
0.1
0.0
-0.1
-0.1
0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
Fig 2 QSAR OF LOWER THE SCIATIC NERVE SORBITOL LEVEL(%)
Possible drug leads: ~ 350 compounds
Aldose Reductase Active Site Structure
CYS298
LEU300
TYP219
TRP111
Cerius2 LigandFit
PHE122
NADPH
HIS110
TRP20
TRP79
TYR48
VAL47
LYS77
To further confirm the AR-ARI binding,
We perform QM/MM calculations on
drug leads.
CHARMM
O
NH
5'-OH, 6'-F, 7'-OH
5' HN
O
6'
NMe
X
7'
8'
N
H
O
Binding energy is found to be –45 kcal / mol
Docking of aldose reductase inhibitor
Aldose reducatse
Inhibitor
(4R)-6’-fluoro-7’-hydroxyl-8’-bromo-3’-methylspiro[imidazoli-dine-4,4’(1’H)-quinazoline]-2,2’,5(3’H)-trione
Cerius2 LigandFit
Hu & Chen, 2003
Interaction energy between ligand and protein
Quantum Mechanics/Molecular Mechanics
(QM / MM)
Hu & Chen, 2003
O
NH
5' HN
O
6'
NMe
X
7'
8'
N
H
O
a:Inhibitor concentration of inhibit Aldose Reductase;
b: the percents of lower sciatic nerve sorbitol levels
c: interaction with AR in Fig. 4
Our Design Strategy
QSAR determination & prediction (Neural Network)
Docking (Cerius2)
QM / MM (binding energy)
?
Software in Department
1. Gaussian
2. Insight II
CHARMm: molecular dynamics simulation, QM/MM
Profiles-3D: Predicting protein structure from sequences
SeqFold:
Functional Genomics, functional identification
of protein w/ sequence and structure comparison
NMR Refine: Structure determination w/ NMR data
3. Games
4. HyperChem
5. AutoDock (docking)
6. MacroModel
6. In-House Developed Software
LODERSTAR
Neural Network for QSAR
Monte Carlo & Molecular Dynamics
Lecture Notes for Physical Chemistry
Year 2
1.Intermediate Physical Chemistry (CHEM2503) (Powerpoint
format .ppt)
Year 3
1.Advanced Physical Chemistry (Powerpoint format .ppt)
2.Electronic Spectroscopy (Powerpoint format .ppt)
3.Electronic Spectroscopy (assignment) (rar file)
Postgraduate Course
1.Research Techniques in Chemistry (Powerpoint format .ppt)
Course Work Download Molecule
M.Sc Course
1.Computational Modeling of Macromolecular Systems
(Powerpoint format .ppt) Download Molecule
HYPERCHEM Exercise
Part A: Study the electronic structure and vibrational spectrum of formaldehyde
Step 1: Build up the structure of the formaldehyde.
1. Run HYPERCHEM software in the start menu.
2. Double click the drawing tool to open the elements table dialogue box and select carbon atom.
Formaldehyde
Procedures:
Close the element table.
(Drawing tool)
O
3. L-click the cursor on the workspace. A carbon atom will appear.
(Make sure drawing tool is selected. R-click on the atom if you want to delete it)
C
4. Repeat (2) and choose oxygen instead of carbon. Move the cursor to the carbon centre and drag the
H atom.) H
mouse from the carbon centre to an empty workspace. (A single bond is created between carbon atom and oxygen
5. L-click the bond between carbon and oxygen to create a double bond.
6. L-click on Build in the menu bar and switch on ‘add H & model build’ (i.e. make sure a tick
appeared on the left of this function.).
Step 2: Optimize the structure using RHF and 6-31G* basis set.
7. L-Click on Setup in the menu bar and L-click ab Initio;
L-Click on 6-31G*; then, L-Click on Options button;
Select RHF, set Charge to 0 and Multiplicity to 1 (default for charge 0);
L-Click OK buttons after modifications were done.
8. L-Click on Compute in the menu bar and select Geometry Optimization;
Select Polak-Ribiere and set RMS gradient to 0.05 and max cycles to 60;
L-Click OK button (The calculation will be started. Repeat the step till “Conv=YES” appears in the status line.).
Record the energy appeared in the status line
9. L-Click on Compute in the menu bar and select Orbitals.
Record energy levels and point groups of required molecular orbitals (MO)
(Optional: You can draw the contour plot of the selected orbital and visualize the shape of the orbital.)
10. L-Click on Compute in the menu bar and select Vibrations.
11. L-Click on Compute in the menu bar and select Vibrational Spectrum.
Record the frequencies of different vibrational modes and their corresponding oscillator strengths.
(Optional: You can turn on animate vibrations, select any vibrational modes, and L-Click on OK button. The molecule begins to vibrate. To
suspend the animation, L-Click on Cancel button.)
Part B: Molecular Dynamics of Tetrapeptide
1. L-click Databases on the menu bar. Choose Amino Acids.
2. Select Beta sheet.
3. L-click Ala, Tyr, Asp and Gly to create tetrapeptide Ala-Tyr-Asp-Gly.
4. L-click on rotate-out-of-plane tool and use it to rotate the molecule to a proper angle for observation and measurements.
(Rotate-out-of-plane tool)
5. L-Click on Setup in the menu bar and L-click Molecular Mechanics;
L-Click on MM+;
L-Click OK buttons after modifications were done.
6. L-Click on Compute in the menu bar and select Geometry Optimization;
7. Record the total energies.
8. L-Click on Compute in the menu bar and L-click Molecular Dynamics;
Run molecular dynamics at 0K and 300K with constant temperature.
Simulation Time: 1ps
9. Record the total energies.
Part C: Molecular Dynamics of Ribosomal Protein
Procedures:
10. Use a web-browser and Go to http://yangtze.hku.hk/lecture_notes.htm.
11. R-click the title labeled “Download molecule” and save it in a folder in your local disk (C:\).
12. L-click on File in the menu bar and select open to load in the molecule.
(You should notify that this file has extension filename .ENT and is in PDB format.)
13. L-click on rotate-out-of-plane tool and use it to rotate the molecule to a proper angle for
observation and measurements.
(Rotate-out-of-plane tool)
14. L-Click on Setup in the menu bar and L-click Molecular Mechanics;
L-Click on MM+;
L-Click OK buttons after modifications were done.
15. L-Click on Compute in the menu bar and L-click Molecular Dynamics;
Run molecular dynamics at 300K with constant temperature.
Simulation Time: 1ps
16. Record the total energy.