
David M. Whisnant and Lisa S. Lever

Computational chemistry is a new field of chemistry that has fully developed in the past ten years or so. Computers have become powerful enough to reliably calculate molecular properties with an accuracy that approaches, or sometimes even exceeds, the quality of experimental results. Computeraided chemistry, if used wisely, can be an important tool in the hands of experimental chemists. These web pages provide the student with an introduction to computational chemistry with information on modern computational methods. You may select a topic from the main menu at right or read consecutively through the pages using the arrow buttons. These reference pages can be printed more concisely using the links below.
 
Computational Results 
Modern computational chemistry puts a broad set of tools at the disposal of chemists. Although they are closely related, we can separate these tools into two classes  molecular modeling and information management.
Molecular modeling involves the development of mathematical models of molecules that can be used to predict and interpret their properties. There are two types of molecular modeling  molecular mechanics and quantum mechanics.
Molecular mechanics is a classical mechanical model that represents a molecule as a group of atoms held together by elastic bonds. In a nutshell, molecular mechanics looks at the bonds as springs which can be stretched, compressed, bent at the bond angles, and twisted in torsional angles. Interactions between nonbonded atoms also are considered. The sum of all these forces is called the force field of the molecule. A molecular mechanics force field is constructed and parameterized by comparison with a number of molecules, for instance a group of alkanes. This force field then can be used for other molecules similar to those for which it was parameterized. To make a molecular mechanics calculation, a force field is chosen and suitable molecular structure values (natural bond lengths, angles, etc.) are set. The structure then is optimized by changing the structure incrementally to minimize the strain energy and spread it over the entire molecule. This minimization is orders of magnitude faster than a quantum mechanical calculation on an equivalent molecule.
Molecular mechanics is a valuable tool for predicting geometries and heats of formation of molecules for which a force field is available. It is a good way to compare different conformations of the same molecule, for instance. Molecular mechanics does have two weaknesses, however. First, force fields are based on the properties of known, similar molecules. If you are interested in the properties of a new type of molecule an appropriate force field probably will not be available for that molecule. Second, because molecular mechanics models look at molecules as sets of springs, they cannot be used to predict electronic properties of molecules, such as dipole moments and spectroscopy. To make predictions about the electronic properties of a molecule, you must use quantum mechanical models.
To make a quantum mechanical model of the electronic
structure of a molecule, we must solve the Schrödinger
equation.
The Hamiltonian operator,
H, depends on the
kinetic and potential energies of the nuclei and electrons in the
atom or molecule. The wavefunction, , will give us information about the probability of
finding the electrons in different places in the molecule. The
energy, E, is related to the energies of individual electrons, which
can be used to help interpret electronic spectroscopy.
Solving the Schrödinger equation is a very difficult problem and
cannot be done without making approximations. Two types of
approximations are the BornOppenheimer approximation and the
independent electron approximation.
In the BornOppenheimer approximation, the positions of
the nuclei are taken to be fixed so that the internuclear distances
are constant. This is a sensible approximation because the massive
nuclei are essentially immobile in comparison with light electrons.
We first choose a geometry (with fixed internuclear distances) for a
molecule and solve the Schrödinger equation for that geometry.
We then change the geometry slightly and solve the equation again.
This continues until we find an optimum geometry with the lowest
energy.
Another approximation (the
independent electron
approximation) commonly made is that the
wavefunction, , can be
written as a product of oneelectron functions. The oneelectron
wavefunctions are called molecular
orbitals 
the molecular equivalent of atomic orbitals. Each
molecular orbital then is expressed as a combination of the atomic
orbitals from the atoms that make up the molecule. For example, the
simplest molecular orbital function for the
H_{2} molecule is
written as
c_{1}1s_{1}
+
c_{2}1s_{2},
where 1s_{i} is a
hydrogen 1s atomic orbital function and
c_{i} is a parameter.
This method is called
LCAOMO theory for
Linear Combination of Atomic Orbitals  Molecular Orbital Theory.
As an example, the two lowest energy molecular orbitals of the
H_{2} molecule are
shown here.
These can be thought of as combinations of 1s
orbitals from the two hydrogen atoms. The molecular orbital on the
left is made by two H atom 1s orbitals combining constructively. This
is a bonding MO
because it helps hold the molecule together. The molecular orbital on
the right is due to the destructive interference of the two 1s
orbitals and is said to be
antibonding.
Although molecular orbitals are written as combinations of atomic
orbitals from the atoms in the molecule,
molecular orbitals are not atomic
orbitals. They are analogous to atomic
orbitals, but instead of being defined for atoms, molecular orbitals
are characteristic of the molecule as a whole.
Because of the large number of particles in a
molecule (benzene, for instance, has 12 nuclei and 78 electrons)
computer programs are used to do the calculations necessary for the
solution of the Schrödinger equation. These calculations involve
an enormous number of difficult integrals for large molecules.
Ab
initio computational methods solve all of
these integrals without approximation. Ab initio methods are the most
reliable for small and mediumsized molecules, but are prohibitively
timeconsuming for large molecules (20 atoms or so for PCs; around
100 atoms if workstations are available ).
For larger molecules,
semiempirical
methods have been developed which ignore or approximate some of the
integrals used in ab initio methods. To compensate for neglecting the
integrals, the semiempirical methods introduce parameters based on
molecular data. Commercial software packages are available for both
ab initio and semiempirical calculations.
Semiempirical molecular orbital theory methods have
been developed which ignore or approximate some of the integrals used
in the solution of the Schrödinger equation. To compensate for
neglecting the integrals, the semiempirical methods introduce
parameters based on molecular data. The most commonly used
semiempirical methods are included in MOPAC (the Molecular Orbital
PACkage) computer programs. The CAChe for Windows software package
includes several computational chemistry tools including a version of
MOPAC. The CAChe MOPAC offers two of the most reliable semiempirical
methods (PM3 and AM1) for predicting heats of formation, ground state
geometries, and ionization potentials. It also includes the ZINDO
program, which does a good job predicting the visibleUV bands for
molecules containing hydrogen and first or second period
elements.
What kind of accuracy do we desire from quantum
mechanical calculations? Ideally a “chemically accurate”
calculation should give the results
below^{2} .
Bond lengths: calculated values within 0. 01  0.02 Å of experiment.
Bond angles: calculated values within 1 2° of experiment.
Semiempirical methods can only be used for elements
for which they have been parameterized, but because these elements
are common ones in organic compounds, semiempirical calculations can
be quite useful. Semiempirical results do not always satisfy the
criteria we set for chemical accuracy. PM3
calculations,^{3} for
instance, generally give bond lengths within ±0.036Å and
bond angles within ±3.9°, not always “accurate”
but still pretty good.
Quantum mechanical calculations also can be used to predict electronic energies and heats of formation. Chemically accurate energies should be within 1 kJ/mol (0.2 kcal/mol), a challenging task.^{4}
The average errors for PM3 calculations of heats of
formation are tabulated
below.^{5}
Type of Compound 

Type of Compound 

All C, H, N, O 
4.4 
Organic cations 
9.5 
Hydrocarbons 
3.6 
Organics with F, Si, Cl, Br, I 
5.7 
Cyclic Hydrocarbons 
2.4 
Compounds with S 
12.1 
Hydrocarbons, double bonds 
2.8 
Compounds with P 
11.5 
Hydrocarbons, triple bonds 
5.6 
Closed shell anions 
8.8 
Aromatic hydrocarbons 
4.1 
Neutral radicals 
7.4 
Organics with N, O 
5.2 


Semiempirical methods should only
be used to predict the properties of molecules for which reliable
parameters are available. Experience also has shown that
semiempirical methods do not do well for problems involving hydrogen
bonding, transition states, and poorly parameterized
molecules^{6}. If you
are interested in the highest accuracy or in problems such as the
these, ab initio calculations should be used if possible.
When you are comparing semiempirical or ab initio predictions with
experiment you should remember that computational models generally
are of gasphase
molecules. There are ways of computationally
modeling molecules in solution, either at the ab initio or
semiempirical level, but we will not get into this subject
here.
Scaling Vibrational Frequencies
In the last part of the job output from a frequency calculation you will find the predicted vibrational frequencies (cm^{1}) of the normal modes of the molecule. Also supplied are the predicted intensities of the IR and Raman bands corresponding to these normal modes.
1 2 3 B1 B2 A1 Frequencies  1335.5948 1383.4094 1679.4157 4 5 6 A1 A1 B2 Frequencies  2027.8231 3160.8817 3232.9970
Computational results usually have systematic errors. In the case of HartreeFock level calculations, for instance, it is known that calculated frequency values are almost always too high by 10%  12%. To compensate for this systematic error, it is usual to multiply frequencies predicted at the HF/631G(d) level by an empirical factor of 0.893. Similarly, frequencies calculated at the MP2/631G(d) level are scaled by 0.943 ^{8}.
The predicted frequencies after applying the 0.893 scale factor are listed below.
1 2 3 B1 B2 A1 Scaled Frequencies  1193 1235 1450 4 5 6 A1 A1 B2 Scaled Frequencies  1811 2822 2887
Molecular orbitals are not real physical quantities.
Orbitals are a mathematical convenience that help us think about
bonding and reactivity, but they are not physical observables. In
fact, several different sets of molecular orbitals can lead to the
same energy. Nevertheless, they are quite useful. We will use
ethylene as an example to illustrate MO concepts.
We often can classify molecular orbitals as sigma or pi orbitals. A
sigma ( ) orbital has
cylindrical symmetry about the internuclear axis. The hydrogen
orbitals below are sigma orbitals. As we saw earlier a bonding
orbital has a high electron probability density between the nuclei.
An antibonding orbital has a node between adjacent nuclei with lobes
of opposite mathematical characteristics (shown as different
colors).
A pi ( )orbital has
“up and down” properties like an atomic porbital. Bonding
and antibonding orbitals for
ethylene are shown below.
A few acetone molecular orbital surfaces displayed by CAChe following a PM3 calculation help illustrate these concepts.
Acetone molecular orbitals:
The first orbital is a orbital because the lobes are pointing at each other along an internuclear axis and have rotational symmetry about that axis. It is antibonding because lobes of a different color are adjacent to each other. The atomic orbitals that contribute to this molecular orbital are not obvious from the picture.
This is the lowest unoccupied
molecular orbital (LUMO). This is a orbital
because the lobes are perpendicular to an internuclear axis. It is
antibonding because
lobes of a different color are adjacent to each other. This molecular
orbital appears to be predominantly made of a porbital on the
central carbon atom and a porbital on the oxygen atom.
This is the highest
occupied molecular orbital (HOMO). It is a
nonbonding orbital
because the molecular orbital lobes are located on atoms that are not
bonded. A nonbonding orbital does not help hold the molecule
together. Electrons in a nonbonding orbital are a little like lone
pairs of electrons in a Lewis structure. You can see that a porbital
on the oxygen atom is a major contributor to this molecular
orbital.
This is a orbital.
It is bonding because
lobes of the same color are next to each other. As you can see, two
porbitals (one on the carbon atom and one on the oxygen) overlap to
form a molecular orbital that distributes electrons around the
central carbon and the oxygen atom.
This is a
bonding orbital.
It has two portions, each of which appears to be formed by the
overlap of a porbital on a carbon atom and an sorbital on one of
the hydrogen atoms.
You may remember from physics that a distribution of
electric charge creates an electric potential in the surrounding
space. A positive electric potential means that a positive charge
will be repelled in that region of space. A negative electric
potential means that a positive charge will be attracted. A molecule
is a collection of charges that will have an electric potential 
commonly called the “electrostatic potential.” The
electrostatic potential is a physical property of a molecule related
to how a molecule is first “seen” or “felt” by
another approaching
species^{7}. A portion
of a molecule that has a negative electrostatic potential will be
susceptible to electrophilic attack  the more negative the better.
It is not as straightforward to use electrostatic potentials to
predict nucleophilic attack.
At right is an electrostatic potential surface of acetone displayed by CAChe using MM/PM3 geometry with PM3 wavefunction. The surface is color coded according to electrostatic potential (blue is negative and red is positive). What part of the acetone molecule appears to be more susceptible to electrophilic attack?
The electron density surface depicts locations around the molecule where the electron probability density is equal. This gives an idea of the size of the molecule and its susceptibility to electrophilic attack.
Below is an electron density surface of acetone displayed by CAChe using MM/PM3 geometry with PM3 wavefunction. The surface color reflects the magnitude and polarity of the electrostatic potential. Gray, violet and blue colors correspond to a negative electrostatic potential  regions of the molecule susceptible to electrophilic attack.
To make a quantum mechanical model of the electronic
structure of a molecule, we must solve the Schrödinger
equation.
Solving this equation is a very difficult problem and cannot be
done without making approximations. In ab initio methods, no integrals are neglected
in the course of the calculation. One approximation is the BornOppenheimer
approximation in which the positions of the nuclei are fixed so that the
internuclear distances are constant. We first choose a geometry (with fixed
internuclear distances) for a molecule and solve the Schrödinger equation
for that geometry. We then change the geometry slightly and solve the equation
again. This continues until we find an optimum geometry with the lowest energy.
When more than one electron is present, the Schrödinger equation is impossible to solve because of the interelectron terms in the Hamiltonian.
Consider, for instance, the Hamiltonian for the
hydrogen molecule in the BornOppenheimer approximation.
The first two terms are due to the kinetic energy of the electrons.
The last six terms express the potential energy of the system of four
particles.
The potential
energy term due to the repulsion of the electrons makes the
Schrödinger equation impossible to solve.
To produce a solvable Schrödinger equation
we assume that the Hamiltonian is a sum of oneelectron functions,
f_{i}, with an
approximate potential energy that takes the average interaction of
the electrons into account. This leads to a set of oneelectron
equations, called the HartreeFock
equations, where is a oneelectron wavefunction.
The total wavefunction that is a solution to the total
Schrödinger equation, , is
approximated as the product of the solutions to the oneelectron
equations.
This product must be adjusted to satisfy the Pauli
Exclusion principle, but we won't get into that here. If you are
familiar with determinants, it involves writing the wavefunction as a
determinant.
The question remains about the approximate potential energy in the oneelectron functions that take the average interaction of the electrons into account. What is the form of the functions f_{i} in the HartreeFock equations? The most common way of handling this is to define
where
v_{i} is an
average potential
energy due to the interaction of one electron with all the other
electrons and nuclei in the molecule. The average potential depends
on the orbitals, , of the other
electrons, which means we must solve the HartreeFock equations
iteratively.
The iterative solution of the HartreeFock equation is as
follows.
1. Guess reasonable oneelectron orbitals (wavefunctions), , and calculate the average
potential energies,
v_{i}.
2. Using the variation principle, solve the HartreeFock
equations,
to give new oneelectron orbitals, . Use these new orbitals to
calculate new and improved average potential energies,
v_{i}. Because the
solution of the HartreeFock equations depends on the variation
principle, the HartreeFock energy should be higher than the true
energy.
3. Repeat the second step until the oneelectron orbitals and
potential energies don't change (are selfconsistent).
To take the Pauli Principle into account, we must
include electron spin
in our wavefunctions. The orbitals that are calculated by the
HartreeFock method actually are spin
orbitals that are a product of a spatial
wavefunction and a spin function.
In a spin orbital, is the spatial wavefunction describing the probability
of finding the electron in space and or are spin
wavefunctions.
For a closed shell
system, in which all of the electrons are
paired, during the solution of the selfconsistent field equations,
we can restrict the solution so that the spatial wavefunctions for
paired electrons are the same. This is called a
restricted
HartreeFock
(RHF) calculation and
generally is used for molecules in which all
the electrons are paired. When the spin
functions are removed, we are left with a set of spatial orbitals,
each occupied by two electrons.
An example would be the restricted
HartreeFock solution to the Schrödinger equation for the
hydrogen molecule, H_{2}. This would lead
to two spatial orbitals, one occupied by the pair of electrons and
one unoccupied. The orbitals holding electrons are called
occupied orbitals and
the unoccupied orbitals are called virtual
orbitals.
For
open shell systems
that contain unpaired
electrons, the assumption made in the
restricted HartreeFock method obviously won't work. There is more
than one way of handling this type of problem. One way is to
not constrain pairs of
electrons to occupy the same spatial orbital  the
unrestricted
HartreeFock
(UHF) method. In this
method there are two sets of spatial orbitals  those with spin up
() electrons and those with spin
down () electrons.
This leads to two sets of orbitals as pictured at the right and to a
lower energy than if the restricted method were used.
For molecular calculations, the HartreeFock SCF
equations
still cannot be
solved without one further approximation. To solve the equations,
each SCF orbital,, is written as
a linear combination of atomic
orbitals. For instance, for the
H_{2} molecule, the simplest
approximation is to write each spatial SCF orbital as a combination
of 1s atomic orbitals, each centered on one of the protons.
This reduces the problem to solving for the coefficients,
c_{1} and
c_{2}, since the atomic orbitals do not
change.
The set of atomic orbitals that is chosen to represent the SCF
orbitals is called a basis
set. The {1s_{A},
1s_{B}} basis set shown above is a
minimal basis set  the smallest set of orbitals possible that
describe an SCF orbital. Usually, the quality
of a basis set depends on its size. For
instance, a larger basis set, such as
{1s_{A},
1s_{B},
2s_{A},
2s_{B}}would do a better job
approximating the SCF orbital than
{1s_{A},
1s_{B}}.
For manyelectron atoms, we don't know the actual mathematical
functions for the atomic orbitals, so substitutes are used  usually
either Slatertype orbitals (STO) or Gaussiantype orbitals (GTO). We
won't concern ourselves with the exact form of STO and GTO. Suffice
it to say that they are chosen to behave mathematically like the
actual atomic orbitals: stype, ptype,
dtype, and ftype, for instance. A few
commonly used basis sets are listed below. The symbol of the basis
set is given in the left column and the characteristics of the basis
set in the center. At the right is the basis set that would be used
to represent methane. For instance, the STO3G basis set for methane
would be {1s_{H},
1s_{H},
1s_{H},
1s_{H},
1s_{C},
2s_{C},
2p_{xC},
2p_{yC},
2p_{zC}}.
Basis Sets^{8} 
Characteristics 
Basis Set Example (CH_{4})  

STO3G 
A minimal basis set (although not the smallest possible) using three GTOs to approximate each STO. This basis set should only be used for qualitative results on very large systems 
Each H: 1s  
321G 
Inner shell basis functions made of three GTOs. Valence s and porbitals each represented by two basis functions (one made of two GTOs, the other of a single GTO). Use for very large molecules for which 631G is too expensive. 
Each H: 1s, 1s'  
631G(d) 
Inner shell basis functions made of six GTOs. Valence s and p orbitals each represented by two basis functions (one made of three GTOs, the other of a single GTO). Adds six dtype basis functions to nonhydrogen atoms. This is a popular basis set that often is used for medium and large systems. 
Each H: 1s, 2s  
631G(d,p) 
Like 631G(d) except ptype functions also are added for hydrogen atoms. Use when hydrogens are of interest and for final, accurate energy calculations. 
Each H: 1s, 2s,
2p_{x},
2p_{y},
2p_{z} 
Generally, the larger the basis set the more accurate the calculation (within limits) and the more computer time that is required. As an example, consider the calculation of the bond length of HF using different basis sets, as shown below^{8}.
Basis Set 
Bond Length (Å) 
 Error (Å)  

631G(d) 
0.93497 
0.017 
631G(d,p) 
0.92099 
0.003 
631+G(d,p) 
0.94208 
0.025 
631++G(d,p) 
0.92643 
0.009 
6311G(d,p) 
0.91312 
0.004 
6311++G(d,p) 
0.91720 
0.000 
Experimental 
0.917 

You might notice that although the large basis set, 6311++G(d,p), predicts the correct answer to within 0.001 Å, several others are correct to within 0.01 Å (well within the criteria of chemical accuracy). Although a larger basis set usually gives better results, you often have diminishing returns as you choose larger sets. A point may be reached beyond which the additional computer time is not worth it.
Even with a very large basis set calculation,
HartreeFock results are not exact because they rely on the
independent electron approximation. HartreeFock
SCF Theory is a good baselevel theory that is reasonably good at
computing the structures and vibrational frequencies of stable
molecules and some transition
states^{9}. Electrons
are not independent, though. We say that they are correlated with
each other and that the HartreeFock method neglects
electron correlation.
This means that HartreeFock calculations do not do a good job
modeling the energetics of reactions or bond dissociation. There are
several ways of correcting SCF results to take electron correlation
into account.
One method of taking electron correlation into
account is MøllerPlesset manybody
perturbation theory, which is used after a
RHF or UHF calculation has been made. It is assumed that the
relationship between the exact and HartreeFock Hamiltonians is
expressed by an additional term,
H^{(1)}, so that
H = f_{i} +
H^{(1)}. Calculations
based on this assumption lead to corrections that can improve SCF
results. Various levels of perturbation theory can be applied to the
problem. They are called
MP2,
MP3,
MP4, etc. MP2
calculations are not timeconsuming and usually give quite accurate
geometries and about onehalf of the correlation energy. Because
perturbation theory is not based on the variation principle, the
energy predicted by MP calculations can fall below the actual
energy.
Another important method of correcting for the correlation energy is
configuration interaction
(CI). Conceptually we can think of CI
calculations as using the variation principle to combine various SCF
excited states with the SCF ground state, which lowers its energy. We
won't use CI calculations in our exercises at this level.
When calculating
molecular
orbitals, you should remember that molecular
orbitals are not real physical quantities. Orbitals are a
mathematical convenience that help us think about bonding and
reactivity, but they are not physical observables. In fact, several
different sets of molecular orbitals can lead to the same energy.
Nevertheless, they are quite useful. We will use ethylene as an
example to illustrate MO concepts.
The basis
functions in SCF molecular orbitals are like
atomic orbitals. A RHF/631G(d) calculation on ethylene uses 38 basis
functions (15 for each carbon and 2 for each hydrogen). Since the
molecular orbital wavefunction is expanded in terms of the all the
basis functions,
it might seem that constructing a picture of the
orbital would be difficult. Luckily, most of the coefficients are
zero, so the molecular orbitals are easy to picture. Consider, for
instance, the highest occupied molecular orbital (HOMO) and lowest
unoccupied molecular orbital (LUMO) of ethylene.
HOMO:
The HOMO is a bondingorbital.
LUMO:
The LUMO is an antibonding orbital.
^{1}J. H. Krieger, “Computational Chemistry Impact”, C&E News, 1997, May 12, p 30.; E. K. Wilson, “Computers Customize Combinatorial Libraries”, C&E News, 1998, April 27, p 31.
^{2}J. B. Foresman and Æ. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, Pittsburgh, 199596, p. 118.
^{3}J. Stewart, J. Comp. Chem., 1989, 10, p. 209220.
^{4}P. W. Atkins; R. S. Friedman, Molecular Quantum Mechanics, 3rd Ed., Oxford, New York, 1997, p. 276.
^{5}M. C. Zerner, Rev. Comp. Chem., Vol. 2, VCH, New York, 1991, pp 313365.
^{6}J. B. Foresman and Æ. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, Pittsburgh, 199596, p. 113.
^{7}P. Politzer and J. S. Murray, Rev. Comp. Chem., Vol. 2, VCH, 1991, p. 273.
^{8}J. B. Foresman and Æ. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, Pittsburgh, 199596, p 102.
^{9}J. B. Foresman and Æ. Frisch, Exploring Chemistry with Electronic Structure Methods, Gaussian, Pittsburgh, 1996, p 115.
This page was last updated on March, 2004.