UNIT II - Computational Chemistry

UNIT II - Computational Chemistry

I. Introduction: Definition - taking known theory, developing and using computer software to solve chemical problems.

A. Categories.

1. Information retrieval:

a. Chemical databases.

b. Internet, www.

2. Spreadsheets for data organization and manipulation.

*3. Molecular graphics and visualization.

4. Numerical methods:

a. Data fitting.

b. Statistical analysis.

5. Computer-Instrument interfacing. (CADA)

6. Robotics.

*7. Molecular modeling:

a. Molecular mechanics using empirical force fields.

b. Quantum mechanics.

B. Classes of Problems.

1. Single molecule calculations:

a. Can give information about the most stable geometry of the molecule; its bond lengths and angles; barriers to internal rotation round single bonds and conformational information; vibrational frequencies; electron distribution; ionization potentials; electron affinities; dipole moments; spin- orbit coupling constants.

b. Such quantities can be computed to an accuracy comparable with experimental accuracy for molecules with up to about twenty atoms.

2. Assemblies of molecules:

a. Large number of molecules in particular solvent molecules; or biomolecules immersed in solution thermodynamic properties can be calculated: free energies, heat capacities, equilibrium constants, binding energies between small molecules and macromolecules.

3. Reactions of molecules:

a. Rate constants for chemical reactions remain perhaps the single biggest unsolved problem.

b. Computer calculations can give an idea of relative rate constants and structures of transition states.

C. Molecular Modeling Overview.

D. Hardware Resources - Requirements increase as follows relative to the number of atoms N in the system:

1. Ab initio QM:

time 5 N⁴

2. Semiempirical QM:

time 5 N² or N³

3. Molecular mechanics:

time 5 N²

E. Method Comparison.

F. Usefulness of Various Methods.

1. Molecular Mechanics Methods:

a. Energy minimization for locating stable conformations.

b. Searching conformation space (conformational analysis).

c. Studying motion of molecules using Molecular Dynamics.

2. Quantum Mechanics Methods:

a. Most stable geometry of small molecules.

b. MO energies and shapes.

c. Heat of formation.

d. Partial atomic charges.

e. Electrostatic potential.

f. Dipole moment.

g. Transition state geometries and energies

h. Bond dissociation energies.

3. Semiempirical Quantum Mechanics: treats up to 100 atoms.

4. Ab initio Quantum Mechanics: treats up to ª 30 atoms.

READING: CS Chem3D Manual - Chapter 9 and Computational Chemistry by Grant and Richards - Chapter 1

II. Quantum Chemistry on Computers.

A. QM Overview.

1. In MM each entity is an atom. In QM each entity is an electron. Thus, QM describes molecules in terms of explicit interactions between electrons and nuclei.

2. QM methods are based on approximate solutions to the Schrodinger equation:

HY = EY

3. Exact solution (i.e. in terms of known analytic functions) is only possible for one-electron atomic species.

e.g. Lowest energy level of Hydrogen-like species:

;

Z = 1 for H

= 2 for He⁺, etc.

= nuclear charge

e.g. 2p atomic orbital:

4. For many-electron atoms we adopt the orbital approximation: each electron is treated separately by placing it into a modified Hydrogen-like atomic orbital. Each e^- interacts with the nucleus and the average field generated by the presence of the other electrons.

e.g. ;

Two-electron wave function = product of two one- electron wave functions

This is what we mean when we write electron configuration of He as:

He = 1s² or 1s^a1s^b, where a, b are two opposite spin directions

or: C = 1s²2s²2p²

Y_He = y₁s(1)y₁s(2) = 1s(1)1s(2)

Z_eff = effective nuclear charge experienced by each electron in the presence of the screening effect of the other e^-.

Table 11.3 Real Hydrogen-like Energy Eigenfunctions.

Here z = Z_eff/n

5. Must enforce antisymmetry condition, since electrons are fermions. Must satisfy the Pauli Principle.

More precisely, Pauli Principle states that Y must change sign if we interchange any pair of electrons:

instead of simply: Y_He = 1s^a(1)1s^b(2)

6. In all future work we will simply write the wave function of a many-electron atom as:

Y = y_1s(1)y_1s(2)y_2s(1)…

keeping in mind that this is not a simple product.

7. Each y atomic Hydrogen-like orbital is a function of r, q, f:

y = (normalization constant) x (exponential function of r) x (spherical harmonics)

8. Important properties of orbitals of a given atom:

- normalization

- orthogonal

If orbital i and j are from different atoms:

- "overlap integral"

9. Molecular Orbitals (MO) make the orbital approximation here also:

Y_molecule = f₁f₂f₃f₄…

and f_i are wave functions of individual electrons in the molecule.

10. Fundamental approximation for molecules is the Born-Oppenheimer approximation:

Due to large mass difference between nuclei and electrons, nuclei are sluggish on the time scale of electron motions. Nuclei are considered motionless relative to electrons.

Therefore, in all computational QM, the wave functions for the electrons (molecular orbitals) are calculated for fixed positions of the nuclei.

11. Each MO is built up approximately by taking a Linear Combination of Atomic Orbitals (LCAO).

i^th MO =

where y_k is an atomic orbital centered on one of the atoms in the molecule.

The sum is typically over all the valence shell atomic orbitals of all the atoms in the molecule.

The larger the set, the more flexible the approximate wave function will be and the more accurate.

12. For example, in a typical quantum calculation on the diatomic molecule HF, the basis set of AO’s used to construct MO’s is:

13. These 5 atomic orbitals are overlapped (i.e. we take linear combinations) in 5 different ways to produce 5 molecular orbitals which will contain the 8 valence electrons of HF. Here is the one of lowest energy:

with E₁ = -49.757 eV

Similar expressions exist for f₂, f₃, f₄ and f₅ with different sets of coefficients and energies.

14. The full printout from a MOPAC semiempirical quantum calculation of HF is included in these notes from which the whole set of coefficients and energies can be extracted. Also given are pictures of these MO wavefunctions in space.

15. The question remains: what procedure is used to determine the values of the weighting coefficients that will produce the most accurate molecular orbitals?

i.e. How are the C_ik determined?

B. Solving the Schrodinger Equation for Molecules.

1. The variation principle forms the basis for making sure we have the "best", most accurate wave function for our molecule.

Approximate MO energy = is always > E_i (true)

So, by varying the coefficients C_ik which make up the MO function f_i until we find the minimum energy possible, we are insured of having the best f_i and energy E_i afforded by the basis set we have chosen. The more robust the basis set, the more accurate the f_i and E_i for the molecule.

2. The actual procedure for the above is described succinctly in "Computational Chemistry" by Grant and Richards, pages 11 - 15. This involves lots of matrix algebra with numbers derived from the following integral expressions:

overlap integrals

Coulomb integrals

3. The S_lk are a snap to solve. The harder problem is the H_lk. Different MO methods vary according to the way they handle this problem.

4. In the Coulomb integrals, H is an effective one- electron Hamiltonian which is complicated by the fact that is accounts for interactions with all the other electrons. H then has within it terms like which gives the electron density of the i^th electron.

Therefore one must "know" the solution f_i before one starts, just to have the effective one-electron Hamiltonians H.

5. This implies an iterative solution called the Hartree- Fock self-consistent field (SCF) iteraction. The Hartree-Fock one electron Hamiltonian is often written as:

Coulomb electronic interactions

exchange terms (due to Pauli)

H^N = attraction to fixed nuclei and simple kinetic energy terms

C. Ab initio Methods.

1. All integrals are actually solved for the basis set employed, then the matrix diagonalization is performed.

The number of integrals needing to be solved increases as N⁴, where N = atoms.

2. Different ab initio methods then differ by the basis set employed.

3. Most intuitively obvious choice for basis set would be Slater-type atomic orbitals. (Hydrogen-like with good z = Z_eff/n values):

y_k= Ce^-zrY_lm(q,f)

Then, e.g. for each Carbon atom we use a 1s, a 2s, and three 2p atomic orbitals.

4. Easier to do integrals if y_k are of a Gaussian form instead of a Slater exponential:

But the true wave function is more nearly like around nuclei, so we might use 3 Gaussian functions combined to approximate the Slater exponential. This is called:

STO - 3G

meaning "Slater-type orbitals - 3 Gaussians" or STO - 4G is even more accurate.

D. Semi-Empirical Methods.

1. Less computationally intensive than ab initio.

2. Many of the integrals are not actually directly solved, but are substituted for by parameters determined empirically, i.e., so as to give good fit to experiments.

3. Particularly the H^N integrals, which treat electron interaction with nuclei, are replaced by H^core, which involves attractions to the inner core rather than a bare nucleus. Then only valence shell electrons are treated explicitly.

4. Some common semi-empirical quantum mechanics methods include:

CNDO complete neglect of differential overlap

INDO intermediate neglect of differential overlap

MINDO

AM1 Austin Model 1

PM3 Parameterized Model revision 3

5. A comparison of these methods is given in the Chem3D manual, p. 125-131.

E. Molecular Geometry Optimization.

1. All QM methods can solve the Schrodinger equation for MO’s and energies (and everything else) for a fixed set of nuclear coordinates input by the user. This is called a single point energy calculation.

2. They are also capable of computing forces on the atoms and moving them in a direction of decreasing total energy until an optimum geometry is reached. This is called geometry optimization, and is done iteratively through many steps.

3. In the HF example, the user inputs starting positions of H and F such that the bond distance was 1.07Å. After 7 cycles of iteration, the energy decreased almost 30 kcal/mol and a final bond distance of 0.827Å was achieved. This is thus the predicted equilibrium bond length of H-F according to the AM1 Hamiltonian.

4. In a geometry optimization of H₂O using the CNDO method, I performed a single-point energy calculation for a series of H-O-H bond angels starting from bent configuration (90°) and going to linear configuration (180°). The following potential energy curve was derived.

Note that CNDO predicts an equilibrium bond angle of 105°, compared with the experimental value of 104°. Further, we see that it costs 40 kcal/mol energy to "straighten out" water to 180°.

F. Molecular Conformation.

1. Most organic and biological molecules have freedom to rotate about bonds and can exist in several conformational states. An example is butane, which rotates around the center C-C bond to produce g₊, trans, and g_- conformations.

2. Conformational analysis is the process of systematically exploring the conformational states of a molecule, at each point doing single-point calculation to generate potential energy surface.

G. Other Quantities Accessible to QM Calculations.

1. Vibrational frequencies may be derived by fitting quadratics to potential energy curve minima. (useful for vibrational spectroscopy)

2. Ionization potentials and electron affinities.

3. Charge distribution and dipole moment.

Mulliken population analysis: e.g. H₂O calculation with AM1 (Dipole moment = 1.66 D)

4. Electrostatic potential maps.

5. Frontier orbitals electron density:

LUMO = lowest unoccupied molecular orbital

HOMO = highest occupied molecular orbital

Used to predict site of reactivity of molecules. Most highly reactive site on a molecule in an electrophilic reaction is the HOMO (the molecule is the e^- donor). For a nucleophilic attack it is the LUMO (the molecule is the e^- acceptor).

READING: CS Chem3D Manual - Chapter 9, p. 125- 131 and Computational Chemistry by Grant and Richards - Chapter 2

III. Molecular Mechanics.

A. The Empirical Energy Function (or Force Field).

1. The fundamental interacting unit is the atom, not individual electrons. (Thousands of atoms can be treated in a calculation.)

2. The potential energy of the collection of atoms can be calculated as a fairly simple function of the atomic coordinates.

3. This function is called the Energy Function, and is derived empirically by giving good fit to experimental spectroscopy data.

4. The Energy Function, sometimes called the Valence Force Field, can be broken down into a sum of important interaction terms. For example:

E_tot = E_b + E_q + E_f + E_nb + E_elec

= bond stretching + bond angle bending + torsion angle rotation + non-bonded interaction + electrostatic interaction

5. Bond stretching energy term:

; the sum is over all covalent bonds.

k_i = Hooke’s law spring constant for bond #i

b_oi = equilibrium bond length for bond #i

b_i = actual current value of bond length i

Therefore: bond stretching is treated as a classical harmonic oscillator term.

6. Bond angle bending term:

; the sum is over all covalent bond angles

k_i = spring constant for angle deformation

q_oi = equilibrium bond angle for angle i

q_i = current value of bond angle i

7. Dihedral (or torsion) angle rotation term:

; the sum is over dihedral angles

V_i = barrier heigth

s = 1 for staggered minima

= -1 for eclipsed minima

n = periodicity (n = 3 for ethane, n = 2 for ethene)

f_i = current value of dihedral angle i

Note: the barriers to rotation in molecules are generally a combination of the dihedral angle term and the 1-4 non-bonded interaction term.

8. Non-bonded interaction terms:

the double sum is over all pairs of atoms separated by more than 2 bonds

Other types of functions are common also.

This is often the time-consuming term in simulating large systems.

9. Electrostatic interaction terms:

q_i = charges

D = Dielectric constant of medium

10. Other terms are added by various simulation packages to improve correlation with experiment:

- Hydrogen bonding

- stretch/bend cross terms

B. Energy Minimization.

1. The primary function of MM calculations is to determine reliable lowest energy geometries for molecules having a large number of atoms. Process is called energy minimization.

2. There are numerous algorithms for taking the starting coordinates and moving the system to a minimum in its energy function:

Steepest Descents, Conjugate Gradient, Newton- Raphson

3. The problem is: one achieves a local minimum in the energy and it may not be the absolute or true minimum energy structure.

This is called the multiple minima problem:

4. One possible solution: put kinetic energy into the model by assigning velocities to the atoms. Allow molecular dynamics to take place (solving Newton’s mechanics) and occasionally stop the trajectory and do an energy minimization. This is called the MD/quenching method.

5. Conformational analysis is readily do-able by MM. This involves styding the range of flexibility available to a molecule along one or several of its coordinates and mapping the potential energy as a function of these coordinates.

For example: Ramachandran plot of peptide dihedral angles

C. Molecular Dynamics.

1. Another primary application of molecular mechanics is molecular dynamics (MD) simulation. Newton’s equations of motion are solved for atomic motion using the force field discussed.

2. Ranges and types of motion can be analyzed.

3. Equilibrium properties of the molecule or system can be time-averaged.

4. MD simulations must have some non-abrupt systematic way of introducing kinetic energy to the atoms. This is the heating/equilibration phase of the simulation, and typically involves "ramping" the kinetic energies up to some value to give a desired temperature.

5. Energy conservation must be maintained also by judiciously rescaling the velocities every so often. This is to offset the cumulative effects of small numerical errors.

6. Small errors also introduce a global translation and rotation of the system which must be subtracted out.