Computational modeling of macromolecular systems
1 / 98

Computational Modeling of Macromolecular Systems - PowerPoint PPT Presentation

  • Uploaded on

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.

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about ' Computational Modeling of Macromolecular Systems' - pisces

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript
Computational modeling of macromolecular systems

Computational Modeling of Macromolecular Systems

Dr. GuanHua CHEN

Department of Chemistry

University of Hong Kong

Computational chemistry
Computational Chemistry

  • Quantum Chemistry

    SchrÖdinger Equation

    H = E

  • Molecular Mechanics

    F = Ma

    F : Force Field

Computational Chemistry Industry



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)


1997 Nobel Prize

in Biology:

ATP Synthase in


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

10-9 to 103 s; 0.5 to >10 Ao

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

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

t: 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 Small molecules

  • AMBER Polymers

  • CHAMM Polymers

  • BIO Polymers

  • OPLS Solvent Effects





Algorithms for Molecular Dynamics

Runge-Kutta methods:

x(t+t) = x(t) + (dx/dt) t

Fourth-order Runge-Kutta

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
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
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


SchrÖdinger Equation

Hy = Ey



H = (-h2/2ma)2 - (h2/2me)ii2

+  ZaZbe2/rab - i Zae2/ria

+ ije2/rij


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) Fock operator for 1

F(2) f(2)+ J1(2) -K1(2) 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

Ffi = ei fi

FFock 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 ckiyk LCAO-MO

{ yk }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





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.

Slater-type orbitals (STO)

nlm = Nrn-1exp(-r/a0) Ylm(,)

x the orbitalexponent

Basis Set i = p cip p

Gaussian type functions (GTF)

gijk = N xi yj zk exp(-ar2)

(primitive Gaussian function)

p = u dupgu

(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 and a set of p functions to hydrogen



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

6-31G for a carbon atom: (10s12p)  [3s6p]

1s 2s 2pi (i=x,y,z)


1CGTF 1CGTF 1CGTF 1CGTF 1CGTF (s) (s) (s) (p) (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

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’

H0yn(0) = En(0) yn(0)

yn(0) is an eigenstate for unperturbed system

H’ is small compared with H0

Moller-Plesset (MP) Perturbation Theory

The MP unperturbed Hamiltonian H0

H0 = mF(m)

whereF(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

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



Coupled-Cluster Doubles (CCD) Method

yCCD= eT2 y(0)

y(0): Hartree-Fock ground state wave function

yCCD: Ground state wave function

T2 : two electron excitation operator



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.


Density-Functional Theory (DFT)

E0 = - (h2/2me)i <i |i2|i >-  drZae2(r) /r1a

+ (1/2)   dr1 dr2e2/r12 + Exc[(r)]

Kohn-Sham Equation:

FKSyi = ei yi

FKS- (h2/2me)ii2-  Zae2 /r1a + jJj + Vxc

Vxc dExc[(r)] / d(r)

Extended Huckel MO Method

(Wolfsberg, Helmholz, Hoffman)

Independent electron approximation

Schrodinger equation for electron i

Hval = iHeff(i)

Heff(i) = -(h2/2m) i2 + Veff(i)

Heff(i) i = i i

Semiempirical Molecular Orbital Calculation


fi= r criyr

s (Heffrs- eiSrs ) csi = 0

Heffrs < r|Heff| s >Srs< r| s >

  • Parametrization:

  • Heffrr < r|Heff| r >

  • = minus the valence-state ionization

  • potential (VISP)

Atomic Orbital Energy VISP

--------------- e5 -e5

--------------- e4 -e4

--------------- e3 -e3

--------------- e2 -e2

--------------- e1 -e1

Heffrs = ½ K(Heffrr + Heffss) SrsK: 13


(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) = rs tu (rr|tt) rs tu 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, INDOare 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.


(Dewar and co-workers, University of Texas,


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
Relativistic Effects

Speed of 1s electron: Zc / 137

Heavy elements have large Z, thus relativistic effects are


Dirac Equation:

Relativistic Hartree-Fock w/ Dirac-Fock operator; or

Relativistic Kohn-Sham calculation; or

Relativistic effective core potential (ECP).

Quantum Mechanical Simulation of Nano-size Systems

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



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ijcni (H-I)ijcnj

- nmijcni (H-I)ijcmj lcnlcml

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)


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 EXCITEDSTATE ?

A Much More Difficult Problem !

Linear-Scaling Calculation for excited states

Localized-Density-Matrix (LDM) Method

r = r(0) + dr


rij(0) = 0 rij > r0

drij = 0 rij > r1

Yokojima & Chen, Phys. Rev. B, 1999

Principle of the nearsightedness of

equilibrium systems (Kohn, 1996)

Heisenberg Equation of Motion

Time-Dependent Hartree-Fock

Random Phase Approximation


PPP Semiempirical Hamitonian

Cambridge Display Technology

Weight: 15 gram

Resolution: 800x236

Size: 45x37 mm

Voltage: DC, 10V






DFT: no or very small gap

Liang, Wang, Yokojima & Chen, JACS (2000)

Smallest SWNT: 0.4 nm in diameter

Wang, Tang & etc., Nature (2000)

Three possibilities:

(4,2), (3,3) & (5,0) SWNTs

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
Quantum Mechanics / Molecular Mechanics (QM/MM) Method

Combining quantum mechanics and

molecular mechanics methods:



Human Genome Project


Drug Discovery

Aldose Reductase

Design of Aldose Reductase Inhibitors

Goddard, Caltech

Goddard, Caltech