Computer Modeling Dr Guan Hua CHEN Department of

Computer Modeling Dr. Guan. Hua CHEN Department of Chemistry University of Hong Kong http: //yangtze. hku. hk/lecture/comput 06 -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 Hyper. Cube Informatix Celera Genomics Software Gaussian 94, Gaussian 98 Jaguar Spartan Q-Chem Insight. II, Cerius 2 Hyper. Chem 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

C 60 Vitamin C energy heme Cytochrome c OH + D 2 --> 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, AM 1 PM 3, CNDO/S, INDO/S

SchrÖdinger Equation H = E Wavefunction Hamiltonian H = (-h 2/2 m ) 2 - (h 2/2 me) i i 2 - i Z e 2/ri (+ Z Z e 2/r ) + i j e 2/rij Energy One-electron terms: (-h 2/2 m ) 2 - (h 2/2 me) i i 2 - i Z e 2/ri Two-electron term: i j e 2/rij

Hartree-Fock Method Orbitals 1. Hartree-Fock Equation F i = ei i F Fock operator i the i-th Hartree-Fock orbital ei the energy of the i-th Hartree-Fock orbital

2. Roothaan Method (introduction of Basis functions) i = 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 e. V Ionization Energy Electron Affinity

A Gaussian Input File for H 2 O # HF/6 -31 G(d) Route section water energy Title 0 1 Molecule Specification O -0. 464 0. 177 0. 0 (in Cartesian coordinates H -0. 464 1. 137 0. 0 H 0. 441 -0. 143 0. 0 Basis Set i = p cip p { k } is a set of atomic orbitals (or basis functions) STO-3 G, 3 -21 G, 4 -31 G, 6 -31 G*, 6 -31 G** ------------------------------------------- complexity & accuracy

Gaussian type functions gijk = N xi yj zk exp(-ar 2) (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 = H 0 + H’ H 0 n(0) = En(0) is an eigenstate for unperturbed system H’ is small compared with H 0

Moller-Plesset (MP) Perturbation Theory The MP unperturbed Hamiltonian H 0 = m F(m) where F(m) is the Fock operator for electron m. And thus, the perturbation H’ = H - H 0 Therefore, the unperturbed wave function is simply the Hartree-Fock wave function . Ab initio methods: MP 2, MP 3, MP 4

Coupled-Cluster Method = e. T (0): Hartree-Fock ground state wave function : Ground state wave function T = T 1 + T 2 + T 3 + T 4 + T 5 + … Tn : n electron excitation operator T 1 =

Coupled-Cluster Doubles (CCD) Method CCD = e. T (0) 2 (0): Hartree-Fock ground state wave function CCD: Ground state wave function T 2 : two electron excitation operator T 2 =

Complete Active Space SCF (CASSCF) Active space All possible configurations

Density-Functional Theory (DFT) Hohenberg-Kohn Theorem: Phys. Rev. 136, B 864 (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) E 0 = - (h 2/2 me) i < i | i 2 | i >- dr Z e 2 (r) / r 1 + (1/2) dr 1 dr 2 e 2/r 12 + Exc[ (r)] Kohn-Sham Equation Ground State: Phys. Rev. 140, A 1133 (1965) FKS i = ei i FKS - (h 2/2 me) i i 2 - Z e 2 / r 1 + j Jj + Vxc d. Exc[ (r)] / d (r) A popular exchange-correlation functional Exc[ (r)]: B 3 LYP

Hu, Wang, Wong & Chen, J. Chem. Phys. (Comm) (2003) B 3 LYP/6 -311+G(d, p) B 3 LYP/6 -311+G(3 df, 2 p) RMS=21. 4 kcal/mol RMS=12. 0 kcal/mol RMS=3. 1 kcal/mol RMS=3. 3 kcal/mol B 3 LYP/6 -311+G(d, p)-NEURON & B 3 LYP/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 -31 G*) HFSCF 1 0 OK DFT ~1 CIS <10 OK CISD 17 80 -90% (20 electrons) CISDTQ very large 98 -99% MP 2 1. 5 85 -95% (DZ+P) MP 4 5. 8 >90% CCD large >90% CCSDT very large ~100%

Search for Transition State: one negative frequency k e-DG/RT DG Reactant Product Reaction Coordinate

Gaussian Input File for Transition State Calculation #b 3 lyp/6 -31 G opt=qst 2 test the first is the reactant internal coordinate 0 1 O H 1 oh 1 2 ohh 1 oh 1 0. 90 ohh 1 104. 5 The second is the product internal coordinate 0 1 O H 1 oh 2 H 1 oh 3 oh 2 0. 9 oh 3 10. 0 ohh 2 160. 0 2 ohh 2

Semiempirical Molecular Orbital Calculation Extended Huckel MO Method (Wolfsberg, Helmholz, Hoffman) Independent electron approximation Hval = i Heff(i) = -(h 2/2 m) i 2 + Veff(i) Schrodinger equation for electron i Heff(i) i = i i

LCAO-MO: i = r cri r eff - e S ) c = 0 ( H s rs i rs si Heffrs < r| Heff | s > Srs < r| s > Parametrization: Heffrr < r| Heff | r > = minus the valence-state ionization potential (VISP)

Atomic Orbital Energy -------e 5 -------e 4 -------e 3 -------e 2 -------e 1 Heffrs = ½ K (Heffrr + Heffss) Srs VISP -e 5 -e 4 -e 3 -e 2 -e 1 K: 1 3

CNDO, INDO, NDDO (Pople and co-workers) Hamiltonian with effective potentials Hval = i [ -(h 2/2 m) i 2 + Veff(i) ] + i j>i e 2 / rij two-electron integral: (rs|tu) = < r(1) t(2)| 1/r 12 | 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, 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, AM 1, PM 3 (Dewar and co-workers, University of Texas, Austin) MINDO: modified INDO MNDO: modified neglect of diatomic overlap AM 1: Austin Model 1 PM 3: MNDO parametric method 3 *based on INDO & NDDO *reproduce the binding energy

Relativistic Effects Speed of 1 s 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) = Eelec(QM/MM) + Evdw(MM) + EMM-bond(MM) EQM(QM) + Eelec(QM/MM) = < | Heff | > Heff = -1/2 i i 2 + ij 1/rij - i Z /ri - i q /ri + g 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) AM 1 Strain, Scuseria & Frisch, Science (1996): LSDA / 3 -21 G DFT calculation on 1026 atom RNA Fragment

Linear Scaling Calculation for Excited State Liang, Yokojima & Chen, JPC, 2000

Fast Multiple Method LDM-TDDFT: Cn. H 2 n+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-, PM 3 -, AM 1 -, INDO/S-, & TDDFT-LDM

Photo-excitations in Light Harvesting System II strong absorption: ~800 nm 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 ( 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

Non-bonding interaction van der Waals interaction for pairs of non-bonded atoms Coulomb potential for all pairs of charged atoms

Force Field Types • • • MM 2 AMBER CHAMM BIO OPLS Molecules Polymers 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) (s 1+2 s 2+2 s 3+s 4) t +O( t 5) s 1 = dx/dt s 2 = dx/dt [w/ t=t+ t/2, x = x(t)+s 1 t/2] s 3 = dx/dt [w/ t=t+ t/2, x = x(t)+s 2 t/2] s 4 = dx/dt [w/ t=t+ t, x = x(t)+s 3 t] Very accurate but slow!

Algorithms for Molecular Dynamics Verlet Algorithm: x(t+ t) = x(t) + (dx/dt) t + (1/2) d 2 x/dt 2 +. . . x(t - t) = x(t) - (dx/dt) t + (1/2) d 2 x/dt 2 -. . . x(t+ t) = 2 x(t) - x(t - t) + d 2 x/dt 2 + O( t 4) 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 Diabetes Glucose Sorbitol Diabetic Complications

Inhibitor Aldose Reductase Design of Aldose Reductase Inhibitors

Descriptors: Electron negativity Volume Database for Functional Groups

Possible drug leads: ~ 350 compounds

Aldose Reductase Active Site Structure CYS 298 LEU 300 TYP 219 TRP 111 Cerius 2 Ligand. Fit PHE 122 NADPH HIS 110 TRP 20 TYR 48 VAL 47 LYS 77 TRP 79

To further confirm the AR-ARI binding, We perform QM/MM calculations on drug leads. CHARMM 5'-OH, 6'-F, 7'-OH Binding energy is found to be – 45 kcal / mol

Docking of aldose reductase inhibitor Aldose reducatse Inhibitor (4 R)-6’-fluoro-7’-hydroxyl-8’-bromo-3’-methylspiro[imidazoli-dine-4, 4’(1’H)-quinazoline]-2, 2’, 5(3’H)-trione Cerius 2 Ligand. Fit Hu & Chen, 2003

Interaction energy between ligand protein Quantum Mechanics/Molecular Mechanics (QM / MM) Hu & Chen, 2003

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 (Cerius 2) QM / MM (binding energy) ?

Software in Department 1. Gaussian 2. Insight II CHARMm: molecular dynamics simulation, QM/MM Profiles-3 D: Predicting protein structure from sequences Seq. Fold: Functional Genomics, functional identification of protein w/ sequence and structure comparison NMR Refine: Structure determination w/ NMR data 3. Games 4. Hyper. Chem 5. Auto. Dock (docking) 6. Macro. Model 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 (CHEM 2503) (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) 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) 4. Repeat (2) and choose oxygen instead of carbon. Move the cursor to the carbon centre and drag the mouse from the carbon centre to an empty workspace. (A single bond is created between carbon atom and oxygen atom. ) 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 -31 G* basis set. 7. L-Click on Setup in the menu bar and L-click ab Initio; L-Click on 6 -31 G*; 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 0 K and 300 K with constant temperature. Simulation Time: 1 ps 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 300 K with constant temperature. Simulation Time: 1 ps 16. Record the total energy.