Vol. 9

COMPUTATIONAL CHEMISTRY

319

COMPUTATIONAL
QUANTUM CHEMISTRY
FOR FREE-RADICAL
POLYMERIZATION

Introduction

Chemistry is traditionally thought of as an experimental science, but recent rapid
and continuing advances in computer power, together with the development of

Encyclopedia of Polymer Science and Technology. Copyright John Wiley & Sons, Inc. All rights reserved.

320

COMPUTATIONAL CHEMISTRY

Vol. 9

efficient algorithms, have made it possible to study the mechanism and kinetics
of chemical reactions via computer. In computational quantum chemistry, one can
calculate from first principles the barriers, enthalpies, and rates of a given chem-
ical reaction, together with the geometries of the reactants, products, and tran-
sition structures. It also provides access to useful related quantities such as the
ionization energies, electron affinities, radical stabilization energies, and singlet–
triplet gaps of the reactants, and the distribution of electrons within the molecule
or transition structure. Quantum chemistry can provide a “window” on the reac-
tion mechanism, and assumes only the nonrelativistic Schr¨odinger equation and
values for the fundamental physical constants.

Quantum chemistry is particularly useful for studying complex processes

such as free-radical polymerization (see R

ADICAL

P

OLYMERIZATION

). In free-radical

polymerization, a variety of competing reactions occur and the observable quanti-
ties that are accessible by experiment (such as the overall reaction rate, the overall
molecular weight distribution of the polymer, and the overall monomer, polymer,
and radical concentrations) are a complicated function of the rates of these individ-
ual steps. In order to infer the rates of individual reactions from such measurable
quantities, one has to assume both a kinetic mechanism and often some additional
empirical parameters. Not surprisingly then, depending upon the assumptions,
enormous discrepancies in the so-called “measured” values can sometimes arise.
Quantum chemistry is able to address this problem by providing direct access
to the rates and thermochemistry of the individual steps in the process, without
recourse to such model-based assumptions.

Of course, quantum chemistry is not without limitations. Since the multi-

electron Schr¨odinger equation has no analytical solution, numerical approxima-
tions must instead be made. In principle, these approximations can be extremely
accurate, but in practice the most accurate methods require inordinate amounts
of computing power. Furthermore, the amount of computer power required scales
exponentially with the size of the system. The challenge for quantum chemists
is thus to design small model reactions that are able to capture the main chemi-
cal features of the polymerization systems. It is also necessary to perform careful
assessment studies, in order to identify suitable procedures that offer a reason-
able compromise between accuracy and computational expense. Nonetheless, with
recent advances in computational power, and the development of improved algo-
rithms, accurate studies using reasonable chemical models of free-radical poly-
merization are now feasible.

Quantum chemistry thus provides an invaluable tool for studying the mech-

anism and kinetics of free-radical polymerization, and should be seen as an impor-
tant complement to experimental procedures. Already quantum chemical studies
have made major contributions to our understanding of free-radical copolymer-
ization kinetics, where they have provided direct evidence for the importance of
penultimate unit effects (1,2). They have also helped in our understanding of
substituent and chain-length effects on the frequency factors of propagation and
transfer reactions (2–5). More recently, quantum chemical calculations have been
used to provide an insight into the kinetics of the reversible addition fragmen-
tation chain transfer (RAFT) polymerization process (6,7). For a more detailed
introduction to quantum chemistry, the interested reader is referred to several
excellent textbooks (8–16).

Vol. 9

COMPUTATIONAL CHEMISTRY

321

Basic Principles of Quantum Chemistry

Ab initio molecular orbital theory is based on the laws of quantum mechanics,
under which the energy (E) and wave function (

) for some arrangement of atoms

can be obtained by solving the Schr¨odinger equation 1 (17).

ˆ

H

 = E

(1)

This is an eigenvalue problem for which multiple solutions or “states” are

possible, each state having its own wave function and associated energy. The low-
est energy solution is known as the “ground state,” while the other higher energy
solutions are referred to as “excited states.” The wave function is an eigenfunction
that depends upon the spatial coordinates of all the particles and also the spin
coordinates. Its physical meaning is best interpreted by noting that its square
modulus is a measure of the electron probability distribution. The term ( ˆ

H) in

equation 1 is called the Hamiltonian operator and corresponds to the total kinetic
( ˆ

T) and potential energy ( ˆ

V) of the system.

ˆ

H

= ˆT + ˆV

(2)

ˆ

T

= −

h

2

8

π

2



i

1

m

i

∂

2

∂x

2
i

+

∂

2

∂y

2
i

+

∂

2

∂z

2
i



(3)

ˆ

V

=



i

<



j



e

i

e

j

r

ij

(4)

In these equations, the kinetic and potential terms are summed over all

particles in the system (nuclei and electrons), and each particle is characterized
by its Cartesian coordinates, its mass (m) and its electric charge (e). In addi-
tion, h is Planck’s constant and r

ij

is the distance separating particles i and j.

It can thus be seen that the only empirical information required to solve the
Schr¨odinger equation are the masses and charges of the nucleus and the electrons,
and the values of some fundamental physical constants. For this reason, quantum-
chemical calculations are referred to as “ab initio” (“without assumptions”)
procedures.

In applying the Schr¨odinger equation to problems of chemical interest, two

major simplifications can be made. Firstly, the Hamiltonian described by equa-
tions 2–4 above is nonrelativistic, and it is valid provided that the velocities of the
particles do not approach the speed of light. This is generally reasonable except
for the inner shell electrons of heavy atoms, and in these cases corrections for
relativistic effects can be made (18). Secondly, the Hamiltonian is further simpli-
fied by neglecting the kinetic energy contribution of the nuclei. This is known as
the Born–Oppenheimer approximation (19), and amounts to assuming that, since
nuclear motion is much slower than electronic motion, the electron distribution
depends only upon the positions of the nuclei and not on their motions. Thus, the
Schr¨odinger equation can be rewritten as an electronic Schr¨odinger equation, as

322

COMPUTATIONAL CHEMISTRY

Vol. 9

follows:

ˆ

H

elec



elec

(r

, R) = E

eff

(R)



elec

(r

, R)

(5)

In this equation,



elec

is the electronic wave function (which depends on

the electronic coordinates r, as well as the nuclear coordinates R) and ˆ

H

elec

is

the electronic Hamiltonian in which the total kinetic energy is summed over all
electrons only. The Born–Oppenheimer approximation is valid provided that the
ratio of electron mass to nuclear mass is sufficiently small.

The objective of computational quantum chemistry is to solve this nonrel-

ativistic, electronic Schr¨odinger equation for E

eff

(R), that is, the energy corre-

sponding to a given arrangement of nuclei. If the Schr¨odinger equation is solved
for all possible arrangements of nuclei in a system, one obtains the potential en-
ergy surface. (In fact, one yields the ground-state potential energy surface and any
number of excited-state surfaces; however, for the remainder of this chapter we
will focus largely on the ground states.) This contains the information required for
the quantitative description of chemical structures and processes. For example,
by identifying the nuclear coordinates corresponding to local minima in the po-
tential energy surface, one can determine the equilibrium geometries of chemical
species. By comparing the energies of alternative local minima, one can determine
the relative energies of alternative conformers and isomers, and thus identify the
preferred (ie, global minimum energy) structure. Alternative local minima on a
potential energy surface may correspond to the reactant(s) and product(s) in a
chemical reaction, and by comparing their total energies (including zero-point vi-
brational energy) one can calculate the reaction enthalpy at 0 K. The transition
structure for a chemical reaction can also be identified from the potential energy
surface as a first-order saddle point, that is, a stationary point in which the energy
is a local maximum in one dimension (corresponding to the reaction coordinate),
and a local minimum in all other dimensions. By comparing the total energies at
the transition structure and the reactants, one can calculate the energy barrier for
the chemical reaction. Having identified the equilibrium geometry of a molecule,
one can calculate the second derivative of the energy at that point in the potential
energy surface, and this yields the vibrational frequencies of the molecule (ie, the
peak positions in its IR and Raman spectra) (see V

IBRATIONAL

S

PECTROSCOPY

). The

vibrational frequencies can be used to calculate the total zero-point vibrational
energy, and also the vibrational contribution to the enthalpy and entropy of a
molecule at any temperature. Finally, the potential energy surface provides the
information required to calculate the rates of chemical reactions.

Solving the Schr¨odinger equation also yields the ground-state wave function

for the given arrangement of nuclei. Since the square modulus of the wave func-
tion is related to the electron probability density, quantum-chemical calculations
thus enable us to determine the distribution of electrons within a molecule, and
hence the charge distribution. In fact, assigning charges to specific atoms within a
molecule raises a number of philosophical problems, as under quantum mechanics
electrons do not “belong” to any specific nucleus but are distributed throughout
the molecule in molecular orbitals. To determine the charge on a specific atom
within a molecule, it is necessary to define a scheme for dividing the molecular
volume into atomic subspaces. A number of such schemes exist, and the main

Vol. 9

COMPUTATIONAL CHEMISTRY

323

ones include the Mulliken population analysis (8), the natural bond orbital analy-
sis (20), and Atoms in Molecules (AIM) theory (21). In very broad terms, the former
two schemes assign electrons to a nucleus if they are located in orbitals centered
on that nucleus, but differ in the manner in which those “orbitals” are defined.
By contrast, AIM theory uses a spatial definition for the atom, and generally pro-
vides the most rigorous and the most physically meaningful charges, though it is
also the most computationally expensive method. A detailed discussion of wave
function analysis methods is provided in Reference (9).

In summary, by solving the Schr¨odinger equation, one has direct access to

electronic-structure information and can calculate from first principles the mech-
anism, kinetics, and thermodynamics of chemical reactions. In principle, such
calculations can be extremely accurate, relying only upon the validity of quan-
tum mechanics, and values for the fundamental physical constants. However, in
practice, the multi-electron Schr¨odinger equation has no analytical solution and
numerical approximations must instead be made. These approximations are a
potential source of error in the calculations, and it is thus important to under-
stand their underlying assumptions. This article is primarily concerned with ab
initio molecular orbital theory, which is one of the principal approaches to solving
the electronic Schr¨odinger equation. In this section, the main principles of ab initio
molecular orbital theory are detailed. Other approaches to solving the Schr¨odinger
equation include density functional theory and semiempirical methods, and these
are also briefly outlined. Finally, the additional theoretical calculations required
in order to use the output of quantum-chemical calculations to obtain the rate
coefficients for chemical reactions are described. The method of molecular me-
chanics (MM)—which is an empirical procedure that is not based on solving the
Schr¨odinger equation—is described elsewhere (22–24).

Ab Initio Molecular Orbital Theory

In ab initio molecular orbital theory, the wave function is approximated using one-
electron functions or “spin orbitals” (

χ). Each spin orbital is a product of a molecu-

lar orbital,

ψ(x, y, z), which depends on the Cartesian coordinates of the electron,

and a spin function,

α or β, which corresponds to the spin angular momentum of

the electron being aligned along the positive or negative z-axes, respectively. The
molecular orbitals (

ψ

i

) are represented mathematically as a linear combination

of a set of N one-electron functions called basis functions (

φ

µ

).

ψ

i

=

N



µ = 1

c

µi

φ

µ

(6)

In this equation, the coefficients c

µi

are called the molecular orbital expan-

sion coefficients, and are optimized during the computational procedure. In chem-
ical terms, one can think of the basis functions as the sets of constituent atomic
orbitals, which mix to form the molecular orbitals of the molecule. In order to
approach the exact solution to the Schr¨odinger equation, an infinite set of ba-
sis functions would be required, as this would introduce sufficient mathematical

324

COMPUTATIONAL CHEMISTRY

Vol. 9

flexibility to allow for a complete description of the molecular orbitals. Of course,
in practice a finite set of basis functions must be chosen, and this introduces a
potential source of error to the calculations.

Having formed a set of N linearly independent molecular orbitals, these

orbitals must then be occupied to form the wave function. In ab initio molecu-
lar orbital theory the wave function is formed as the determinant of a matrix,
which for an n-electron closed-shell system [ie, n is even (If n were odd, and the
system was a doublet species (having one unpaired electron), we could instead
form (n

+ 1)/2 orbitals, one of which would be singly occupied. The treatment

of open-shell systems is discussed in more detail below.)] might be written as
follows.

 = (n!)

− 1/2











ψ

1

(1)

α(1) ψ

1

(1)

β(1) ψ

2

(1)

α(1) ψ

2

(1)

β(1) · · · ψ

n

/2

(1)

α(1) ψ

n

/2

(1)

β(1)

ψ

1

(2)

α(2) ψ

1

(2)

β(2) ψ

2

(2)

α(2) ψ

2

(2)

β(2) · · · ψ

n

/2

(2)

α(2) ψ

n

/2

(2)

β(2)

..

.

..

.

..

.

..

.

..

.

..

.

..

.

ψ

1

(n)

α(n) ψ

1

(n)

β(n) ψ

2

(n)

α(n) ψ

2

(n)

β(n) · · · ψ

n

/2

(n)

α(n) ψ

n

/2

(n)

β(n)











(7)

This is usually abbreviated as follows.

 = |χ

1

(1)

χ

2

(2)

. . .χ

n

(n)

|

(8)

In this determinant, which is known as the Slater determinant (25), the first

row corresponds to all possible assignments of electron 1 to all of the spin orbitals,
the second to all possible assignments of electron 2, and so on. The factor of (n!)

− 1/2

ensures that the total probability of finding an electron anywhere in space is 1.

We have already seen that, from our set of N basis functions (where N should

approach infinity), we can form N linearly independent molecular orbitals. How-
ever, in an n-electron system, we choose only n/2 of these orbitals (if n is even)
to occupy in our Slater determinant. Clearly the most appropriate orbitals to
occupy should be the n/2 lowest energy orbitals, and this is the basis of the
well-known Hartree–Fock theory (26). While Hartree–Fock theory performs ad-
equately in many cases, its use of a single determinant wave function can also
frequently lead to considerable error. The problem is that this approach fails to ac-
count for Coulombic electron correlation. That is, it is assumed that each electron
“sees” the other electrons as an average electric field and thus no instantaneous
electron–electron interactions are included. To approach the exact solution to the
Schr¨odinger equation, the wave function must instead be represented as a linear
combination of (an infinite number of) Slater determinants. In each of the addi-
tional determinants, one or more of the lowest energy orbitals are substituted with
the (previously unoccupied) higher energy orbitals. In chemical terms, one might
think of this as allowing for contributions to the wave function from the various
possible excited configurations. This approach is known as configuration interac-
tion (CI), and when all possible excited configurations are included the method is
known as Full CI.

Under ab initio molecular orbital theory, the exact solution to the Schr¨odinger

equation could be obtained using the Full CI method in conjunction with an infinite
basis set. Since this is impractical, computational methods must place restrictions

Vol. 9

COMPUTATIONAL CHEMISTRY

325

Fig. 1.

Pople diagram (8) illustrating the dependence of the accuracy of a computational

method on the basis set and the treatment of correlation.

on the number of basis functions included in the calculation, and (in general) sim-
plify the method for treating correlation. The combination of method and basis set
chosen for a calculation forms the “level of theory,” and by convention is written as
“method/basis set” (using the standard abbreviations for the particular method
and basis set chosen). By increasing the size of the basis set and/or improving the
method, one can improve the accuracy of the calculation but also its computational
cost (see Fig. 1). Furthermore, the performance of a given level of theory can vary
considerably depending upon the chemistry of the system, and the type of prop-
erties being calculated. The choice of level of theory is therefore very important,
and a qualitative understanding of the approximations made at the various lev-
els of theory is thus helpful in choosing appropriate theoretical procedures. For
the remainder of the present section, the main qualitative features of commonly
used methods and basis sets are described. The accuracy and applicability of these
methods for studying the reactions of relevance to free-radical polymerization is
discussed in a following section.

Basis Sets.

LCAO Scheme.

A basis set is a set of one-electron functions, which are

combined to form the molecular orbitals of the chemical species. This is known
as the Linear Combination of Atomic Orbitals (LCAO) scheme. To approach the
exact solution to the Schr¨odinger equation, an infinite set of basis functions would
be required, as this would introduce sufficient mathematical flexibility to allow
for a complete description of the molecular orbitals. In practical calculations, we
must use a finite number of basis functions, and it is thus important to choose basis
functions that allow for the most likely distribution of electrons within the system.
This is achieved using basis functions that are based on the atomic orbitals of the
constituent atoms of the molecule. For example, if a chemical system contained
an oxygen atom, the chosen basis set would include functions describing each of
the 1s, 2s, and three 2p orbitals of an oxygen atom.

326

COMPUTATIONAL CHEMISTRY

Vol. 9

The basis functions that are typically used in practice are called Gaussian-

type orbitals, and examples of their functional form for s-type and p

y

-type orbitals

are shown below.

g

s

(

α, r) =



2

α

π

3

/4

e

− αr

2

(9)

g

y

(

α, r) =



128

α

5

π

3

1

/4

y e

− αr

2

(10)

The exponent

α determines the size of the orbital, and standard values have

been determined for the different orbitals on the different atoms. Gaussian-type
orbitals such as those shown in equations 9 and 10 are called primitive Gaussians.
Contracted Gaussians are also used, and these consist of linear combinations of
primitive Gaussians:

φ

µ

=



p

d

µp

g

p

(11)

The coefficients in this expression (d

µp

) are fixed for the basis set of a given

atom, and should not be confused with the molecular orbital expansion coefficients
(C

µi

) of equation 6, which are determined for a given system during the ab initio

calculation.

Basis sets thus contain sets of primitive or contracted Gaussians, which

correspond to the atomic orbitals of the constituent atoms. Minimal basis sets
typically contain exactly the number of functions required to accommodate the
electrons in the system while maintaining the overall spherical symmetry. For
example, a minimal basis set would contain a 1s function for each hydrogen atom,
1s, 2s, and three 2p functions for each carbon atom, and so on. However, in cases
where the low lying unoccupied orbitals are frequently involved in bonding (eg,
the 2p orbitals of Li or Be), popular minimal basis sets (such as STO-3G) include
these additional functions as well.

The problem with minimal basis sets is that the size and shape of the atomic

orbitals are fixed, and all that can be varied in the quantum chemical calculation
is their contribution to the overall molecular orbitals. However, in reality, the size
and shape of the atomic orbitals (and especially the valence orbitals) can depend
heavily upon the molecular environment. For example, in polar environments, we
might expect a greater degree of asymmetry in the 2p orbitals of oxygen, com-
pared with the same orbitals in nonpolar environments. One way of resolving this
problem would be to create individual basis sets for each atom in each conceiv-
able molecular environment, but this would not only be impractical, it would also
defeat the purpose of performing ab initio calculations. Instead, a practical solu-
tion to this problem is to include additional basis functions of varying sizes and
shapes. By varying their individual contribution to the overall molecular orbitals,
one effectively introduces flexibility to the size and shape of the constituent atomic
orbitals. Some of the main extensions to a minimal basis set are outlined in the
following.

Vol. 9

COMPUTATIONAL CHEMISTRY

327

Fig. 2.

Demonstration of how the mixing of two basis functions of different size is equiv-

alent to introducing a single basis function with variable size. The basis functions are
represented schematically as the boundary surface within which there is a 90% probability
of finding the electron. In the first example, two s-type basis functions are mixed, while in
the second example two p-type functions are mixed.

(1) Double zeta, triple zeta, and quadruple zeta basis sets include additional

sets of each basis function (2, 3, or 4 sets respectively), with each additional
set having a different size. By varying the relative contribution of the al-
ternative sizes to the overall molecular orbital, one effectively introduces a
single function with a variable size (see Fig. 2).

(2) Split-valence basis sets are a simplification to the double, triple, and quadru-

ple zeta basis sets described above. Since the inner shell orbitals are not usu-
ally involved in bonding, and their energies are reasonably independent of
their molecular environment, it is usually only necessary to include the ex-
tra basis functions for the valence orbitals. Basis sets that include different
numbers of basis functions for the inner shell and valence electrons are
known as split-valence basis sets.

(3) Polarization functions are typically included in basis sets to allow for asym-

metry in the electron distribution (see Fig. 3). They also allow for the par-
ticipation of the low lying unfilled atomic orbitals in bonding. Polarization
functions typically consist of basis functions corresponding to the low ly-
ing unfilled atomic orbitals. For example, the basis set for an oxygen atom
might typically include a set of d-type orbitals, while those for hydrogen
might include a set of p-type orbitals. Larger basis sets often include ad-
ditional polarization functions of higher angular momentum. For example,
the 6-311

+ G(3df,2pd) basis set includes a set of f -type functions and three

sets of d-type functions for each first row atom, as well as two sets of p-type
functions and a set of d-type functions for each hydrogen atom.

Fig. 3.

Demonstration of the effect of polarization functions. In the first example an s-type

and a p-type function are mixed, while in the second case a p-type and a d-type function are
mixed. The mixing of higher angular momentum basis functions allows for an asymmetric
distribution of electrons.

328

COMPUTATIONAL CHEMISTRY

Vol. 9

(4) Diffuse functions are basis functions that have very large amplitudes, far

from the atomic nucleus. Sets of diffuse functions are typically added to the
basis set when describing species (such as anions) where the electrons are
not held very tightly to the nucleus.

(5) Effective core potentials (ECPs) are typically used to replace the basis func-

tions of the inner shell electrons for atoms beyond the third row. This re-
duces the computational cost of the calculations, and is possible because the
inner shell electrons on such heavy atoms are relatively unaffected by the
molecular environment. Effective core potentials also include corrections
for relativistic effects, which are significant for the inner shell electrons of
heavy atoms. One of the most commonly used ECPs is called LANL2DZ.

Two commonly used families of extended basis sets are the Pople basis sets

(8) and the Dunning basis sets (27). It is worth making a few brief comments on
their notation.

(1) Pople basis sets have names such as “6-311

+ G(3df,2p).” The “6-311” part

refers to the fact that it is a split valence set with one copy of each ba-
sis function on the inner shell electrons, and three copies on the valence
electrons. The “6” refers to the fact that the basis functions on the inner
shell electrons consist of a contracted Gaussian-type orbital formed from 6
primitive Gaussians. The “311” part refers to the fact that one set of basis
functions on the valence electrons are contracted Gaussians, each formed
from 3 primitive Gaussians, while the other two sets of basis functions are
primitive Gaussians. The “

+” part refers to the fact that one set of s-type

and p-type diffuse functions have been included for each heavy (ie, bigger
than hydrogen) atom. The 6-311

++ G basis set includes an additional set

of s-type diffuse functions on hydrogen atoms as well. The functions in the
bracketed part of the expression (3df,2p) are the polarization functions. The
“3df” part refers to the fact that 3 sets of d-type functions and one set of f -
type functions are included for each heavy atom, and “2p” implies that two
sets of p-type functions are included for each hydrogen atom. Finally, the
notations 6-31G

∗

and 6-31G

∗∗

are also frequently used. The first “

∗

” is an

abbreviation for (d) and indicates that a set of d-type functions are included
for each non-hydrogen atom, while “

∗∗

” stands for (d,p) and indicates that in

addition to the d-type functions for the non-hydrogen atoms, a set of p-type
functions are included for each hydrogen atom.

(2) Dunning basis sets have names such as “cc-pVnZ.” This notation stands for

“correlation-consistent polarized valence n-zeta.” For a double zeta basis
set, n is replaced by a “D,” for a triple zeta basis set, n is replaced by a “T,”
for a quadruple zeta basis set, n is replaced by a “Q,” for a quintuple basis
set we use a “5,” and for a sextuple basis set we use a “6.” When diffuse
functions are included, an “aug” prefix is included in the name, as in “aug-
cc-pVTZ.” The cc-pVTZ basis set generally has a performance similar to
6-311G(2df,p). A special feature of the Dunning basis sets is that they have
been designed so the series DZ, TZ, QZ, 5Z, 6Z

. . . systematically converges

on the infinite basis set limit (27). This feature has been exploited in the

Vol. 9

COMPUTATIONAL CHEMISTRY

329

infinite basis set extrapolation procedure of Martin and Parthiban (28) (see
below).

Plane Waves.

While standard ab initio calculations use the LCAO scheme

almost exclusively, density functional theory (DFT) calculations (discussed later)
sometimes use plane wave basis sets. These are solutions of the Schr¨odinger equa-
tion for a free particle and take the general form,

g

PW

= exp[ikr]

(12)

where the vector 

k is related to the momentum

p of the wave through p =
k. One

of the advantages of plane wave basis sets is that they can handle calculations with
periodic boundary conditions, and they are thus used widely in solid-state physics.
However, a disadvantage is that large basis sets are normally required to achieve
acceptable accuracy and, as a result, plane waves are rarely used in calculations
of molecular systems. For a detailed discussion of plane wave basis sets, and their
application to Car–Parrinello (29) ab initio direct dynamics techniques, the reader
is referred to a review by Bl¨ochl and co-workers (30).

Methods.

Hartree–Fock (HF).

The foundation of ab initio molecular orbital theory

is the Hartree–Fock (HF) method. As we saw above, it is based on a single-
determinant wave function (eq. 7) in which the electrons are assigned to the lowest
energy orbitals. In fact the Slater determinant of equation 7 applies to closed-shell
systems, that is, the n electrons occupy the n/2 orbitals in pairs of opposite spin.
This method is known more specifically as restricted Hartree–Fock (RHF). For
open-shell systems (that is, those having one or more unpaired electrons), two
approaches are possible. In restricted-open-shell Hartree–Fock (ROHF), the de-
terminant is formed from a set of molecular orbitals, which are either doubly or
singly occupied, according to the multiplicity of the species. For example, for a rad-
ical (doublet) species, the determinant would be formed from (n

+ 1)/2 orbitals,

and one of these would be singly occupied. This ensures that there is exactly one
unpaired spin in the system, and the species is indeed a pure doublet. In unre-
stricted Hartree–Fock (UHF), the

α and β spin orbitals are defined and optimized

separately. Thus there would be (n

+ 1)/2 occupied α spin orbitals, and (n −1)/2

occupied

β spin orbitals. The principal differences between the RHF, ROHF, and

UHF theory are illustrated in Figure 4.

The UHF method offers some advantages over the ROHF method. In particu-

lar, the additional freedom in the wave function (with the provision for noninteger
natural orbital occupation numbers) allows it to account in part for nondynamic
electron correlation, and leads to lower energies and a better qualitative descrip-
tion of bond dissociation (9). Through its introduction of spin polarization effects,
the UHF method can also provide a better (though by no means perfect) qualitative
treatment of hyperfine coupling constants (31). However, a disadvantage of UHF
is that the independent optimization of the

α and β spin orbitals can result in

nominally equivalent

α and β orbitals having slightly different eigenvalues. In

other words, a doublet species could have effectively more than one unpaired
electron in the system. This physically unrealistic phenomenon is known as spin

330

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 4.

Electron configuration diagrams highlighting the differences between restricted

Hartree–Fock theory (RHF), restricted open-shell Hartree–Fock theory (ROHF), and un-
restricted Hartree–Fock theory (UHF).

contamination, and can be a particular problem for the transition structures in the
propagation steps of free-radical polymerization. Assessment studies (32,33) for
such reactions have revealed that spin-contaminated UHF wave functions make
poor starting points for M¨oller–Plesset (MP) perturbation theory calculations (see
below). By contrast, MP calculations based on ROHF wave functions (such as
ROMP2) show improved agreement with higher level values. Interestingly, it has
been found that for high level methods such as coupled cluster theory (see below),
the choice of the starting wave function (ie, UHF vs ROHF) makes little difference
to the final calculated energies and geometries for radical reactions (34).

Having constructed a wave function, the remaining unknown parameters in

the Schr¨odinger equation are the molecular orbital expansion coefficients (c

µi

),

as defined in equation 6. Determining the optimum values of these coefficients is
thus the principal task of an ab initio calculation. It is beyond the scope of this
chapter to outline the mathematical equations and numerical algorithms used to
achieve this, but it is worth describing the general approach in qualitative terms.
The basis of this procedure is the variational principle. It merely states that for
any antisymmetric normalized function of the electronic coordinates (eg, a Slater
determinant), the energy of this function is always greater than the expectation
value of the exact wave function for the ground state. In other words, the exact
wave function serves as a lower bound to the energies calculated from our ap-
proximate wave function, and the optimal coefficients c

µi

are merely those that

minimize the energy. This principle leads to the Roothaan–Hall equations (35),
which are solved iteratively until the c

µi

coefficients converge. At convergence

the coefficients are “self-consistent,” and hence the HF theory is also known as
self-consistent field (SCF) theory.

Configuration Interaction (CI).

The description of the wave function using

a single determinant (as in HF theory) fails to take electron correlation into ac-
count. To obtain the exact solution to the Schr¨odinger equation we instead need to
construct the wave function as a linear combination of determinants, with each ad-
ditional determinant corresponding to one of the various possible excited configu-
rations obtained when electrons are promoted to the previously unoccupied higher
energy orbitals (see Fig. 5). The resulting wave function is written as follows.

Vol. 9

COMPUTATIONAL CHEMISTRY

331

Fig. 5.

Electron configuration diagrams showing the configurations corresponding to the

HF wave function, and the various possible excited configurations. For this 4-electron sys-
tem, the combination of single, double, triple, and quadruple excitations constitutes full
CI. However, with an infinite basis set, there would be an infinite number of unoccupied
orbitals (instead of the two shown here) to promote the electrons into, and thus an infinite
number of determinants would be required to obtain the exact solution to the Schr¨odinger
equation.

 = a

0



0

+



s

>0

a

s



s

(13)

The



0

wave function is the HF wave function, while the various



s

determi-

nants correspond to the various excited configurations. The CI method introduces
a further set of unknown parameters into the calculation, the coefficients (a

s

).

These coefficients are optimized as part of the ab initio calculation in order to
minimize the energy, in line with the variational principle. CI methods can be
based on an RHF wave function (RCI), a UHF wave function (UCI), or an ROHF
wave function (URCI).

Full CI is impractical with an infinite basis set (and hence an infinite num-

ber of virtual orbitals), or indeed with a finite basis set and a reasonably small
number of electrons. For example, even for water with the small 6-31G(d) basis
set, the full CI treatment would involve nearly 5

× 10

8

configurations. For this

332

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 6.

Electron configuration diagrams illustrating the lack of size consistency in trun-

cated CI. In the first case, A and B are treated separately by CID, and the treatment thus
considers the double excitations of electrons in each molecule. In the second case, A and
B are calculated as a supermolecule having the A and B fragments at (effectively) infinite
separation. Now the simultaneous excitation of two electrons from each of the A-type and
B-type orbitals constitutes a quadruple excitation, which is not included in the CID method.

reason, methods based on a truncated CI procedure are generally used in prac-
tice. These methods consider a limited number of excited determinants, such as
all possible single excitations (CIS) or all possible single and double excitations
(CISD). Restricting the CI procedure to single, double, and possibly triple excita-
tions is usually a reasonable approximation, since excitations involving one, two,
or three electrons have a considerably higher probability of occurring, and thus
contributing to the wave function, compared with excitations of several electrons
simultaneously.

However, simple truncated CI methods suffer from a lack of size consistency.

That is, the error incurred in calculating molecules A and B separately is differ-
ent from that incurred in calculating a single species, which contains A and B
separated by a large (effectively infinite) distance. This can be seen quite clearly
in the example shown in Figure 6. The lack of size consistency can be a major
problem as it introduces an additional error to calculations of barriers and en-
thalpies in nonunimolecular reactions. This problem is addressed by including
additional terms in the wave function, and the methods based on this approach
include quadratic configuration interaction (QCI) and coupled cluster theory (CC).
These methods are typically applied with single and double excitations (QCISD
or CCSD), and the triple excitations are often included perturbatively, leading to
methods such as QCISD(T) and CCSD(T). When applied with an appropriately
large basis set, these methods usually provide excellent approximations to the
exact solution to the Schr¨odinger equation. However, these methods are still very
computationally expensive.

M ¨oller–Plesset (MP) Perturbation Theory.

By convention, the correlation

energy is simply the difference between the Hartree–Fock energy and the exact
solution to the Schr¨odinger equation. Rather than approximate the exact solu-
tion to the Schr¨odinger equation by attempting to build the exact wave function
through configuration interaction, an alternative (and considerably less expen-
sive approach) is to estimate the correlation energy as a perturbation on the

Vol. 9

COMPUTATIONAL CHEMISTRY

333

Hartree–Fock energy. In other words, the exact wave function and energy are
expanded as a perturbation power series in a perturbation parameter

λ as

follows.

Ψ

λ

= 

(0)

+ λ

(1)

+ λ

2



(2)

+ λ

3



(3)

+ · · ·

(14)

E

λ

= E

(0)

+ λE

(1)

+ λ

2

E

(2)

+ λ

3

E

(3)

+ · · ·

(15)

Expressions relating terms of successively higher orders of perturbation are

obtained by substituting equations 14 and 15 into the Schr¨odinger equation, and
then equating terms on either side of the equation. Having obtained these ex-
pressions, it simply remains to evaluate the first terms in the series, and this is
achieved by taking the



(0)

term as the Hartree–Fock wave function. In practice,

the MP series must be truncated at some finite order. Truncation at the first order
(ie, E

(1)

) corresponds to the Hartree–Fock energy, truncation at the second order

is known as MP2 theory, truncation at the third order as MP3 theory, and so on.
MP methods based on an RHF, UHF, or ROHF wave function are referred to as
RMP, UMP, or ROMP respectively.

When truncated at the second, third, or possibly fourth orders, the MP meth-

ods offer a very cost-effective method for estimating the correlation energy. They
are also size-consistent methods. However, the validity of truncating the series
at some finite order depends on the speed of convergence of the series, and this
will vary considerably depending on how closely the Hartree–Fock energy approx-
imates the exact energy. Indeed in some cases, the MP series can actually diverge,
and the application of MP methods can in such cases increase rather than decrease
the errors in the calculation. As noted above, a relevant example of this problem
occurs in the transition structures for radical addition to alkenes for which UMP2
calculations (based on the spin-contaminated UHF wave function) are frequently
subject to large errors (32,33). Furthermore, when truncated at some finite order,
the MP methods are not variational, and may thus overestimate the correction
to the energy. Hence, although MP procedures frequently provide excellent cost-
effective performance, they must be applied with caution.

Composite Procedures.

The use of CCSD(T) or QCISD(T) methods with a

suitably large basis set generally provides excellent approximations to the exact
solution of the Schr¨odinger equation. However, such methods are computationally
expensive, and in practical calculations smaller basis sets and/or lower cost meth-
ods must be adopted. A major advance in recent years has been the development
of high level composite procedures, which approximate high level calculations
through a series of lower level calculations. Some of the main strategies that are
used are described in the following.

Firstly, it has long been realized that geometry optimizations and frequency

calculations are generally less sensitive to the level of theory than are energy cal-
culations. For example, as will be discussed in a following section, detailed assess-
ment studies (36,37) have shown that even HF/6-31G(d) can provide reasonable
approximations to the considerably more expensive CCSD(T)/6-311

+ G(d,p) level

of theory, for the geometries and frequencies of the species in radical addition to
multiple bonds (such as C C, C C, and C S). By contrast, very high levels of

334

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 7.

Illustration of the relative performance of the high and low levels of theory for

geometry optimizations and energy calculations. The low level of theory shows a very large
error for the absolute energy of structure, a smaller error for the Y X bond dissociation
energy (ie, the well depth), and a very small error for the optimum geometry of the Y X
bond. This reflects the increasing possibility for cancelation of error. In the bond dissociation
energy, errors in the absolute energies of the isolated Y

•

and X

•

species are canceled to some

extent by errors in the Y X energies. In the geometry optimizations, further cancelation
is possible because the position of the minimum energy structure depends on the relative
energies of Y X compounds having very similar Y X bond lengths.

theory are required to describe the barriers and enthalpies of these reactions. The
improved performance of low levels of theory in geometry optimizations and fre-
quency calculations can be understood in terms of the increased opportunity for
the cancelation of error, as such quantities depend only upon the relative energies
of very similar structures (see Fig. 7). In contrast, reaction barriers and enthalpies
depend upon the relative energies of the reactants and transition structures or
products, and these can have quite different structures, with different types of
chemical bonds. It is thus possible to optimize the geometry of a compound at a
relatively low level of theory, and then improve the accuracy of its energy using
a single higher level calculation (called a “single point”). Since geometry opti-
mizations and frequency calculations are more computationally intensive than
single-point energy calculations, this approach leads to an enormous saving in
computational cost. By convention, the final composite level of theory is written
as “energy method/energy basis set//geometry method/geometry basis set.”

Secondly, an extension to the above strategy is known as the IRCmax (in-

trinsic reaction coordinate) procedure. It was developed (38) for improving the ge-
ometries of transition structures, though techniques based on the same principle
have also been used to calculate improved imaginary frequencies and tunneling
coefficients (39–41). While low levels of theory are generally suitable for optimiz-
ing the geometries of stable species, the geometries of transition structures are
sometimes subject to greater error at these low levels of theory. To address this
problem, the minimum energy path (MEP) for a reaction is first calculated at a
low level of theory, and then improved via single-point energy calculations at a
higher level of theory. Now, the transition structure is simply the maximum energy
structure along the MEP for the reaction. By identifying the transition structure
from the high level MEP (rather than the original low level MEP), one effectively
optimizes the reaction coordinate at the high level of theory (see Fig. 8).

Vol. 9

COMPUTATIONAL CHEMISTRY

335

Fig. 8.

Illustration of the IRCmax procedure. The minimum energy path (MEP, also

known as the intrinsic reaction coordinate or IRC) is optimized at a low level of theory,
and then improved using high level single-point energy calculations. The improved tran-
sition structure is then identified as the maximum in the high level MEP. This effectively
optimizes the reaction coordinate (often the most sensitive part of the geometry optimiza-
tion) at a high level of theory.

Thirdly, one can improve the single-point energy calculations themselves

using additivity and/or extrapolation procedures. In the former case, the energy
is first calculated with a high level method (such as CCSD(T)) and a small basis
set. The effect of increasing to a large basis set is then evaluated at a lower level
of theory (such as MP2). The resulting basis set correction is then added to the
high level result, thereby approximating the high level method with a large basis
set. The calculation may be summarized as follows.

High Method

/Small Basis Set

+Low Method/Large Basis Set

− Low Method/Small Basis Set

≈High Method/Large Basis Set

(16)

Procedures for extrapolating the energies obtained at a specific level of the-

ory to the corresponding infinite basis set limit have also been devised. The two
main procedures are the extrapolation routine of Martin and Parthiban (18), which
takes advantage of the systematic convergence properties of the Dunning DZ, TZ,
QZ, 5Z,

. . . basis sets, and the procedure of Petersson and co-workers (42), which

is based on the asymptotic convergence of MP2 pair energies. For the mathemat-
ical details of these extrapolation routines, the reader is referred to the origi-
nal references. The Martin extrapolation procedure is easily implemented on a
spreadsheet, while the Petersson extrapolation procedure has been coded into the
GAUSSIAN (43) computational chemistry software package.

336

COMPUTATIONAL CHEMISTRY

Vol. 9

Building on these strategies, several composite procedures for approximating

CCSD(T) or QCISD(T) energies with a large or infinite basis set have been devised.
The main families of procedures in current use are the G3 (44), Wn (28), and CBS
(42) families of methods. These are described in the following.

(1) In the G3 methods, the CCSD(T) or QCISD(T) calculations are performed

with a relatively small basis set, such as 6-31G(d), and these are then cor-
rected to a large triple zeta basis set via additivity corrections, carried out
at the MP2 and/or MP3 or MP4 levels of theory (44). There are many vari-
ants of the G3 methods, depending upon the level of theory prescribed for
the geometry and frequency calculations, the methods used for the basis
set correction, and depending on whether CCSD(T) or QCISD(T) is used
at the high level of theory. Of particular note are the RAD variants (45)
of G3 (such as G3-RAD and G3(MP2)-RAD), which have been designed for
the study of radical reactions. G3 methods include an empirical correction
term, which has been estimated against a large test set of experimental
data, and spin-orbit corrections (for atoms). The G3 methods have been ex-
tensively assessed against test sets of experimental data (including heats
of formation, ionization energies, and electron affinities) and are generally
found to be very accurate, typically showing mean absolute deviations from
experiment of approximately 4 kJ

· mol

− 1

.

(2) In the Wn methods, high level CCSD(T) calculations are extrapolated to

the infinite basis set limit using the extrapolation routine of Martin and
Parthiban (28). Additional corrections are included for scalar relativistic
effects, core-correlation, and spin-orbit coupling in atoms. No additional
empirical corrections are included in this method. The Wn methods are
very high level procedures, and have been demonstrated to display chemical
accuracy. For example, the W1 procedure was found to have a mean absolute
deviation from experiment of only 2.5 kJ

· mol

− 1

for the heats of formation of

55 stable molecules. For the (more expensive) W2 theory, the corresponding
deviation was less than 1 kJ

· mol

− 1

.

(3) In the CBS procedures, the complete basis extrapolation procedure of Pe-

tersson and co-workers is used (42). This calculates the infinite basis set
limit at the MP2 level of theory. This is then corrected to the CCSD(T)
level of theory using additivity procedures, as in the G3 methods. The CBS
procedures also incorporate an empirical correction, and an additional (em-
pirically determined) correction for spin contamination. The accuracy of this
latter term for the transition structures of radical addition reactions has re-
cently been questioned (36,37). Nonetheless, the CBS procedures also show
similar (excellent) performance to the G3 methods, when assessed against
the same experimental data for stable molecules (42).

In summary, using composite procedures, high level calculations can now be

performed at a reasonable computational cost. With continuing rapid increases
in computer power, details on the computational speeds of the various methods
would be rapidly outdated. However, it is worth noting that, at the time of writing,
the most cost-effective G3 procedures can be routinely applied to molecules as big

Vol. 9

COMPUTATIONAL CHEMISTRY

337

as CH

3

SC

•

(CH

2

Ph)SCH

3

, while the state-of-the-art Wn methods are restricted to

smaller molecules, such as CH

3

CH

2

CH(CH

3

)

•

. However, in the near future one

can look forward to applying these methods to yet larger systems. In general, the
composite procedures described above offer “chemical accuracy” (usually defined
as uncertainties of 4–8 kJ

· mol

− 1

), with the best methods offering accuracy in

the kJ range. However, careful assessment studies are nonetheless recommended
when applying methods to new chemical systems. A brief discussion of the per-
formance of computational methods for the reactions of relevance to free-radical
polymerization is provided in a following section.

Multireference Methods.

The post-SCF methods discussed above are all

based on a HF or single configuration starting wave function. At the impracti-
cal limit of performing full CI (or summing all terms in the MP series) with an
infinite basis set, these methods will yield the exact solution to the nonrelativis-
tic Schr¨odinger equation. However, when truncated to finite order, the use of a
single reference wave function can sometimes lead to significant errors. This is
particularly the case in the calculation of diradical species (such as the transition
structures for the termination reactions in free-radical polymerization), excited
states, and unsaturated transition metals. In such situations, the starting wave
function itself should be represented as a linear combination of two or more con-
figurations, as follows.

 =



j

a

j



j

(17)

In this equation, the individual wave functions are formed from the lowest en-

ergy configuration, and various excited configurations of the Slater determinants,
and the a

j

coefficients are optimized variationally. While this method, which is

known as multireference self-consistent field (MCSCF), may seem analogous to
the single-reference CI methods discussed above, there is an important difference
between them. In MCSCF, the molecular orbital coefficients (the c

µi

in eq. 6) are op-

timized for all of the contributing configurations. In contrast, in single-reference
methods, the molecular orbital coefficients are optimized for the Hartree–Fock
wave function, and are then held fixed at their HF values.

The optimization of both the orbital coefficients and the contribution of the

various configurations to the overall wave function can be very computationally
demanding. As a result, MCSCF methods typically only consider a small number
of configurations, and one of the key problems is choosing which configurations
to include. In complete active space self-consistent field (CASSCF), the molecular
orbitals are divided into three groups: the inactive space, the active space, and
the virtual space (see Fig. 9). The wave function is then formed from all possible
configurations that arise from distributing the electrons among the active orbitals
(ie, full CI is performed within the active space). It then remains to decide which
occupied and virtual orbitals should be included in the active space. Where pos-
sible, it is advisable to include all valence orbitals in the active space, together
with an equivalent number of virtual orbitals. However, as with any full CI calcu-
lation, the computational cost rapidly increases with the number of electrons and
orbitals included, and CASSCF calculations are currently limited to active spaces
of approximately 16 electrons in 16 molecular orbitals. Thus, for large chemical

338

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 9.

The partitioning of orbitals between the inactive, active, and virtual spaces in a

CASSCF calculation.

systems, full valence active spaces are not as yet possible, and instead a restricted
number of orbitals must be chosen.

In non-full valence CASSCF, the active space is typically selected on the

basis of “chemical intuition,” and might include orbitals that are directly involved
in the chemical reaction, or are interacting strongly with the reacting orbitals.
For example, in the case of radical–radical reactions, a simplified multireference
approach would be a CAS(2,2) method, in which the active space would consist
of the two singly occupied molecular orbitals. However, such restricted methods
must be used cautiously as they recover correlation in the active space, but not
in the inactive space or between the active and inactive spaces. As a result, such
procedures can sometimes introduce a bias, which, for example, might lead to
an overestimation of the biradical character in systems with nearly degenerate
singlet and triplet states (9).

Multireference methods primarily account for nondynamic electron corre-

lation, which arises from long-range interactions involving nearly degenerate
states. It is still necessary to correct for dynamic correlation, which arises from
short-range electron–electron interactions, and which is primarily addressed in
the single-reference post-SCF methods, such as QCISD or MP2. For example,
in the case of CI-based methods, it would be necessary to consider excitations
from the inactive (as well as active) space orbitals, into all of the virtual orbitals.
Multireference versions of post-SCF methods have been derived, including mul-
tireference CI (eg, MR(SD)CI, which includes all single and double excitations)
and multireference perturbation theory (eg CASPT2, which is a multireference
analogue of MP2). The former of these methods is generally more accurate but

Vol. 9

COMPUTATIONAL CHEMISTRY

339

also more computationally demanding. For more information on multireference
methods, the reader is referred to an excellent review by Schmidt and Gordon
(46).

Semiempirical Methods

Semiempirical methods are often used to study large systems for which ab initio
calculations are not feasible. A number of different procedures are available, with
the main methods being CNDO, INDO, MNDO, MINDO/3, AM1, and PM3. The
latter two procedures are generally the best performing of the current available
methods, and are thus the most popular in current use. Semiempirical methods
are based on ab initio molecular orbital theory, but neglect several of the computa-
tionally intensive integrals that are required in Hartree–Fock theory. Depending
on the procedure, certain interactions between orbitals are either completely ne-
glected or replaced by parameters that are either derived from experimental data
for the isolated atoms or obtained by fitting the calculated properties of molecules
to experimental data. This greatly reduces the computational cost of the calcu-
lations; however, it can also introduce enormous errors. For more detailed infor-
mation on the principles and limitations of semiempirical methods, the reader is
referred to an article by Stewart (47).

In general, the semiempirical methods perform reasonably well, provided

that the species (and properties) being calculated are very similar to those for
which the method was parameterized. However, there are many situations in
which these methods fail dramatically, and hence such methods should be ap-
plied with caution and their accuracy should always be checked against high level
calculations for prototypical reactions. In this context it should be noted that
such testing has already been performed for the case of radical addition to C C
bonds (32). In this work, semiempirical methods were shown to fail dramatically,
and hence (current) semiempirical methods are not generally recommended for
studying the kinetics and thermodynamics of the propagation step in free-radical
polymerization.

Density Functional Theory

Density functional theory (DFT) is a different quantum chemical approach to
obtaining electronic-structure information. The basis of DFT is the Hohenberg–
Kohn theorem (48), which demonstrates the existence of a unique functional for
determining the ground-state energy and electron density exactly. In the ab initio
methods described above, we recall that their objective was to determine the wave
function (an eigenfunction) of a system, which thereby enables the energy (the
eigenvalue) and electron density (the square modulus of the wave function) to
be evaluated. The Hohenberg–Kohn theorem implies that the electronic energy
can be calculated from the electron density and there is thus no need to evaluate
the wave function. This represents an enormous simplification to the calculation
since, in an n-electron system, the wave function is a function of 3n variables,
whereas the electron density is a function of just 3 variables. Unfortunately, the

340

COMPUTATIONAL CHEMISTRY

Vol. 9

Hohenberg–Kohn theorem is merely an existence proof, rather than a constructive
proof, and thus the exact functional for connecting the energy and electron density
is not known. Hence, although in principle DFT can provide the exact solution to
the Schr¨odinger equation, in practice an approximate functional must be used,
and this introduces error to the calculations.

The DFT methods used in practice are based on the equations of Kohn and

Sham (49). They partitioned the total electronic energy into the following terms.

E

= E

T

+ E

V

+ E

J

+ E

X

+ E

C

(18)

In this equation, E

T

is the kinetic energy term (arising from the motions of the

electrons), E

V

is the potential term (arising from the nuclear–electron attraction

and the nuclear–nuclear repulsion), E

J

is the electron–electron repulsion term, E

X

is the exchange term (arising from the antisymmetry of the wave function), and E

C

is the dynamical correlation energy of the individual electrons. The sum of the E

T

,

E

V

, and E

J

terms corresponds to the classical energy of the charge distribution,

while the exchange and correlation terms account for the remaining electronic
energy. The task of DFT methods is thus to provide functionals for the exchange
and correlation terms. As a matter of notation, DFT methods are typically named
as exchange functional–correlation functional, using standard abbreviations for
the various functionals.

Before discussing the functionals themselves, it is worth making a few com-

ments on unrestricted Kohn–Sham theory (50). The effective potential of the Kohn–
Sham equations contains no reference to the spin of the electrons, and the energy is
simply a functional of the total electron density. (It will only become a functional of
the individual spin densities if the potential itself contains spin-dependent parts,
such as it would in the presence of an external magnetic field.). Hence, if the exact
functional were available, there would normally be no need to consider the

α and

β spin densities individually, even for open-shell systems. However, in practice we
must use approximate functionals, and (for open-shell systems) these are gener-
ally more flexible if they explicitly depend on the

α and β spin densities. In an

analogous manner to UHF, unrestricted Kohn–Sham theory allows the

α and β

spin densities to optimize independently, and this allows for a better qualitative
description of bond-breaking processes but leads to physically unrealistic spin
densities and symmetry breaking problems. A more detailed discussion of the ad-
vantages and disadvantages of the unrestricted and spin-restricted theories may
be found in the excellent textbook by Koch and Holthausen (11), while an example
of a practical application of unrestricted Kohn–Sham theory to reactions with bi-
radical transition structures may be found in a paper by Goddard and Orlova (51).
On balance, the unrestricted Kohn–Sham theory is normally preferred for open-
shell systems; however, as always, careful assessment studies are recommended
in order to establish the suitability of any computational method for a specific
chemical problem.

Since the exact functional relating the energy to the electron (or spin) density

is unknown, it is necessary to design approximate functionals, and the accuracy
of a DFT method depends on the suitability of the functionals employed. Many
different functionals for exchange and correlation have been proposed, and it is
beyond the scope of this article to outline their mathematical forms (these may

Vol. 9

COMPUTATIONAL CHEMISTRY

341

be found in textbooks such as Refs. 10 and 11), but it is worth mentioning their
main assumptions. Pure DFT methods may be loosely classified into “local” meth-
ods and “gradient-corrected” methods. The local DFT methods are based on the
local density approximation (LDA, also known as the local spin density approxi-
mation, LSDA), in which it is assumed that the electron density may be treated
as that of a “uniform electron gas.” From this assumption, functionals describing
the exchange (called the “Slater” functional, S) (52) and correlation (Vosko–Wilk–
Nusair, VWN) (53) can be derived, and the resulting method is known as S-VWN.
The treatment of electron density as that of a uniform electron gas is of course an
oversimplification of the real situation, and, while it can often provide reasonable
molecular structures and frequencies, the LDA model fails to provide accurate
predictions of thermochemical properties such as bond energies (for which errors
of over 100 kJ

· mol

− 1

are typical) (54).

Gradient-corrected DFT methods (also sometimes referred to as nonlocal or

semilocal DFT) attempt to deal with the shortcomings of the LDA model through
the generalized gradient approximation (GGA). This corrects the uniform gas
model through the introduction of the gradient of the electron density. In intro-
ducing the gradient, empirical parameters are often incorporated. For example,
the Becke-88 exchange functional (55) was parameterized against the known ex-
change energies of inert gas atoms. This is commonly used in combination with the
“LYP” gradient-corrected correlation functional (56), to give the B-LYP method.
Another example of a gradient-corrected functional is the Perdew-Wang 91 (PW91)
functional, which has both an exchange and a correlation component (57). The
GGA methods show improved performance over the LDA model, especially with
respect to thermochemical properties. In this regard, the errors generally obtained
in standard thermochemical tests of these methods are of the order of 25 kJ

· mol

− 1

(54). However, the GGA methods (as well as the LDA methods) perform poorly for
weakly bound systems (such as those in which Van der Waal’s interactions are
important), and they also perform poorly for reaction barriers (54).

The DFT methods described above are pure DFT methods. Another impor-

tant class of methods is called hybrid DFT. In these methods the exchange func-
tional is replaced by a linear combination of the Hartree–Fock exchange term and
a DFT exchange functional. In addition, the various exchange and correlation
functionals may themselves be constructed as linear combinations of the vari-
ous available methods. For example, the popular hybrid DFT method, B3-LYP, is
defined as follows (58).

E

XC
B3LYP

= E

X
LDA

+ c

0



E

X
HF

− E

X
LDA



+ c

X



E

X
B88

− E

X
LDA



+ E

C
VWN

+ c

C



E

C
LYP

− E

C
VWN



(19)

The coefficients in this expression, c

0

= 0.20, c

X

= 0.72 and c

C

= 0.81, were

obtained by fitting the results of B3-LYP calculations to a test set of experimental
atomization energies, electron affinities, and ionization potentials.

Hybrid DFT methods frequently provide excellent descriptions of the ge-

ometries, frequencies, and even reaction barriers and enthalpies for many chem-
ical systems. However, owing to their empirical parameters, such methods are
increasingly becoming semiempirical in nature and as such can frequently fail
when applied to systems other than those for which they were parameterized. A
good example of this is the hybrid DFT method MPW1K (59). This was fitted to

342

COMPUTATIONAL CHEMISTRY

Vol. 9

a test set of hydrogen abstraction barriers, and performs very well for these re-
actions. However, the same method has recently been shown to have large errors
when applied to the problem of predicting the enthalpies for radical addition to
multiple bonds (36,37). Nonetheless, hybrid DFT methods currently present the
most cost-effective option for studying larger chemical systems but, as always, the
performance of such methods should be carefully assessed for each new chemical
problem.

Calculation of Reaction Rates from Quantum-Chemical Data

Quantum and Classical Reaction Dynamics.

In the quantum-

chemical calculations described above, we solve the electronic Schr¨odinger equa-
tion to determine the energy corresponding to a fixed arrangement of nuclei. If
such calculations are performed for all possible nuclear coordinates in a chemi-
cal system, this yields the potential energy surface. However, as we saw above, in
constructing this potential energy surface we made the Born–Oppenheimer ap-
proximation, and thus ignored the contribution of the motions of the nuclei to the
total kinetic energy. This approximation was appropriate for calculating the elec-
tronic energy at a specific geometry, but is clearly not very useful for studying the
motions of the atoms in chemical reactions. In order to calculate reaction rates,
we must construct a new Hamiltonian in which the kinetic energy of the nuclei is
taken into account. In this Hamiltonian, the potential energy is simply the total
electronic energy, which we obtain from our quantum-chemical calculations. Once
we have formed our new Hamiltonian we can then solve the Schr¨odinger equa-
tion again, this time to follow the motion of the nuclei. This procedure is known
as quantum dynamics, and can in principle yield the exact reaction rates for a
chemical system (within the Born–Oppenheimer approximation).

However, there are several practical limitations to quantum dynamics.

Firstly, we have already seen that, for all but the simplest chemical systems, ob-
taining accurate solutions to the electronic Schr¨odinger equation for a single set
of nuclear coordinates is very computationally intensive. Secondly, to construct a
potential energy surface, these expensive calculations must be repeated for “ev-
ery” possible arrangement of nuclei. Efficient algorithms are available for choos-
ing only those geometries necessary for an adequate description of the chemical
system (60). However, even using these algorithms, large numbers of quantum-
chemical calculations are nonetheless required. For example, approximately 1000
quantum-chemical calculations were required to construct a reliable potential en-
ergy surface for the OH

+ H

2

system (61). Furthermore, the number of data points

required increases substantially with the number of atoms in the system (due to
the increasing dimensionality). Finally, we have the problem of solving the nuclear
Schr¨odinger equation. In practice, this is intractable for all but the simplest sys-
tems, as atoms (being heavier) require even more basis functions than are needed
to solve the electronic Schr¨odinger equation. With current available computing
power, quantum dynamics calculations are thus restricted to very small systems,
such as OH

+ H

2

(61). In this (state-of-the-art) 4-atom calculation, the energies at

approximately 10

7

different points on the potential energy surface were required

Vol. 9

COMPUTATIONAL CHEMISTRY

343

in order to solve the nuclear Schr¨odinger equation, and this requirement scales
exponentially with the number of atoms in the system.

Classical reaction dynamics provides a strategy for calculating the rate co-

efficients of larger chemical systems. Having used quantum-chemical techniques
to calculate the potential energy surface, the motions of the nuclei are studied by
solving the laws of classical or Newtonian dynamics. This is often a reasonable
approximation, since the atoms (being heavier) are considerably less subject to
quantum effects than the electrons. Nonetheless, standard classical reaction dy-
namics calculations are still limited by the need to calculate a full potential energy
surface (including the first and second derivatives at each point). As a result, stan-
dard classical dynamics calculations involving accurate ab initio potential energy
surfaces are also currently restricted to relatively small chemical systems, such
as H

3

C

3

N

3

(62).

An alternative approach to constructing the entire potential energy surface

for a chemical system is provided by direct dynamics. In both standard classical
reaction dynamics and direct dynamics, the basic principle is the same. A starting
arrangement of atoms is adjusted (by a small amount) according to the forces act-
ing on them during a small “step” in time, using the laws of classical mechanics.
This “time step” is then repeated using the force corresponding to the new geom-
etry, and so on. The process is repeated for many thousands of time steps, until a
complete trajectory is mapped out. The process is then repeated for many trajec-
tories until the reaction probability (and other related information) is established
to within an acceptable level of statistical error. As we saw above, in standard
classical reaction dynamics, the force acting on the molecule as a function of ge-
ometry is obtained from the potential energy surface. In direct dynamics, also
known as on-the-fly ab-initio dynamics, this force is calculated (using quantum-
chemical calculations) at each new position (63). The latter approach is simpler,
but less computationally efficient, and is still restricted to relatively small systems
(if accurate levels of theory are used to calculate the forces).

It should be noted that the classical reaction dynamics of much larger sys-

tems can be studied using approximate potential energy surfaces, constructed us-
ing empirical or semiempirical procedures. In particular, the method of molecular
mechanics (MM), which is described elsewhere in this Encyclopedia, is commonly
used to simulate the motion of polymers and proteins. However, the accuracy of
MM simulations are limited by the accuracy of the “force field,” which is the set
of potential functions that are used to govern relative motions of the constituent
atoms. Force fields are typically derived on the basis of empirical and semiempir-
ical information, and are typically only accurate for the type of system for which
they were parameterized. Recently, much effort has been directed at deriving ac-
curate force fields for reacting systems, and prominent examples include ReaxFF
(64) and MMVB (65). However, such force fields are nonetheless approximate,
and only suitable for the types of reactions for which they were designed. Ac-
curate force fields for studying the kinetics of free-radical polymerization do not
currently exist, and instead high level ab initio calculations are necessary in order
to model these reactions accurately.

Transition-State Theory.

To study reactions in larger chemical systems

using accurate ab initio calculations we need a much simpler approach, and this
is provided by transition-state theory (66). In its simplest form, it assumes that, in

344

COMPUTATIONAL CHEMISTRY

Vol. 9

the space represented by the coordinates and momenta of the reacting particles, it
is possible to define a dividing surface such that all reactants crossing this plane
go on to form products, and do not recross the dividing surface. The minimum
energy structure on this dividing plane is referred to as the transition structure
of the reaction. Transition-state theory also assumes there is an internal statisti-
cal equilibrium between the degrees of freedom of each type of system (reactant,
product, or transition structure), and that the transition state is in statistical
equilibrium with the reactants. In addition, it assumes that motion through the
transition state can be treated as a classical translation. From these assumptions,
the following simple equation relating the rate coefficient at a specific tempera-
ture, k(T), to the properties of the reactant(s) and transition state can be derived
(66).

k(T)

= κ(c

◦

)

1

− m

k

B

T

h

Q

‡



reactants

Q

i

e

− E

0

/RT

(20)

In this equation

κ is called the transmission coefficient and is taken to be

equal to unity in simple transition-state theory calculations, but is greater than
unity when tunneling is important (see below), c

◦

is the inverse of the reference

volume assumed in calculating the translational partition function (see below), m
is the molecularity of the reaction (ie, m

= 1 for unimolecular, 2 for bimolecular,

and so on), k

B

is Boltzmann’s constant (1.380658

× 10

− 23

J

· molecule

− 1

· K

− 1

),

h is Planck’s constant (6.6260755

× 10

− 34

J

·s), E

0

(commonly referred to as the

reaction barrier) is the energy difference between the transition structure and
the reactants (in their respective equilibrium geometries), Q

‡

is the molecular

partition function of the transition state, and Q

i

is the molecular partition function

of reactant i.

Transition-state theory thus reduces the problem of calculating the potential

energy surface for “every” possible geometric arrangement of nuclei, to the con-
sideration of a very small number of “special” geometries; namely, the transition
structure and the reactant(s). The transition structure is the minimum energy
structure on the dividing surface between the reactants and products, and must
be located so as to make the “no re-crossing” assumption as valid as possible. In
simple transition-state theory, the transition structure is located as the maximum
energy structure, along the minimum energy path connecting the reactants and
products. This is generally a good approximation for reactions having barriers
that are large compared to k

B

T. However, for reactions with low or zero barriers,

a more accurate approach is required. To this end, in variational transition-state
theory, the transition structure is located as the structure (on the minimum en-
ergy path) that yields the lowest reaction rate. In thermodynamic terms, this may
be thought of as the maximum in the Gibb’s free energy of activation, rather than
the maximum internal energy of activation.

In order to calculate reaction rates via transition-state theory, one needs

to identify the equilibrium geometries of the reactants, and also the transition
structure, and calculate their energies. This information is of course accessible
from quantum-chemical calculations. The molecular partition functions for these

Vol. 9

COMPUTATIONAL CHEMISTRY

345

species are also required. These serve as a bridge between the quantum mechan-
ical states of a system and its thermodynamic properties, and are given by

Q

=



i

g

i

exp



−

ε

i

k

B

T

(21)

The values

ε

i

are the energy levels of a system, each having a number of

degenerate states g

i

, and are obtained by solving the Schr¨odinger equation. In

theory, this equation should be solved for all active modes but in practice the
calculations can be greatly simplified by separating the partition function into
the product of the translational, rotational, vibrational, and electronic terms, as
follows.

Q

= Q

trans

× Q

rot

× Q

vibr

× Q

elec

(22)

This is generally a reasonable assumption, provided that the reaction occurs

on a single electronic surface. Finally, if we assume that reacting species are ideal
gas molecules, analytical expressions for the partition functions are as follows:

Q

trans

= V



2

πMk

B

T

h

2

3

/2

=

RT

P



2

πMk

B

T

h

2

3

/2

(23)

Q

vib

=



i

exp



−

1
2

h

ν

i

kT

×



i

1

1

− exp



−

h

ν

i

kT



(24)

Q

rot, linear

=

1

σ

r



T



r

where



r

=

h

2

8

π

2

Ik

B

(25)

Q

rot, nonlinear

=

π

1

/2

σ

r



T

3

/2

(



r

,x



r

,y



r

,z

)

1

/2

where



r

,i

=

h

2

8

π

2

I

i

k

B

(26)

Q

elec

= ω

0

(27)

In equations 23–27 R is the universal gas constant (8.314 J

· mol

− 1

· K

− 1

);

M is the molecular mass of the species; V is the reference volume, and T and P
the corresponding reference temperature and pressure:

ν

i

are the vibrational fre-

quencies of the molecule; I is the principal moment of inertia of a linear molecule,
while for the nonlinear case I

x

, I

y

, and I

z

are the principal moments of inertia

about axes x, y, and z respectively;

σ

r

is the symmetry number of the molecule

which counts its number of symmetry equivalent forms (67); and

ω

0

is the elec-

tronic spin multiplicity of the molecule (ie,

ω

0

= 1 for singlet species, 2 for doublet

species, etc). The information required to evaluate these partition functions is rou-
tinely accessible from quantum-chemical calculations: the moments of inertia and
symmetry numbers depend on the geometry of the molecule, while the vibrational
frequencies are obtained from the second derivative of the energy with respect to
the geometry.

346

COMPUTATIONAL CHEMISTRY

Vol. 9

A number of additional comments need to be made concerning the use of

equations 23–27. Firstly, in the calculation of the translational partition function
(eq. 23), a reference volume (or equivalently, a temperature and pressure) is as-
sumed. This is needed for the calculation of thermodynamic quantities such as
enthalpy and entropy, but the assumption has no bearing on the calculated rate
coefficient, as the reference volume is removed from equation 20 through the pa-
rameter c

◦

(

= 1/V). Secondly, the vibrational partition function (eq. 24) has been

written as the product of two terms. The first of these corresponds to the zero-point
vibrational energy of the molecule, while the latter corresponds to its additional
vibrational energy at some nonzero temperature T. The zero-point vibrational
energy is often included in the calculated reaction barrier E

0

. When this is the

case, this first term must be removed from equation 24, so as not to count this en-
ergy twice. Thirdly, the external rotational partition function is calculated using
equation 25 if the molecule is linear, and equation 26 if it is not.

It is also worth noting that there is an entirely equivalent thermodynamic

formulation of transition-state theory.

k(T)

= κ

k

B

T

h

(c

◦

)

1

− m

e

S

‡

/R

e

− H

‡

/RT

(28)

A derivation of this expression, which is obtained by noting the relationship

between the thermodynamic properties of a system (eg enthalpy, H, and entropy,
S) and the partition functions, can be found in textbooks on statistical thermo-
dynamics (12–16). The enthalpy of activation (

H

‡

) for this expression can be

written as the sum of the barrier (E

o

), the zero-point vibrational energy (ZPVE),

and the temperature correction (

H

‡

).

H

‡

= E

0

+ ZPVE + H

‡

(29)

The temperature correction (

H) and ZPVE for an individual species can

be calculated from the vibrational frequencies as follows.

ZPVE

= R

1
2



i

h

ν

i

/k

B

(30)

H = R



i

h

ν

i

/k

B

exp(h

ν

i

/k

B

/T) − 1

+

5
2

RT

+

3
2

RT

(31)

In equation 31, the first term is the vibrational contribution to the enthalpy,

the second term is the translational contribution, and the third term is the rota-
tional contribution. The entropy of activation (

S

‡

) is calculated from the vibra-

tional (S

v

), translational (S

t

), rotational (S

r

), and electronic (S

e

) contributions to

the entropies of the individual species, in turn expressed as follows.

S

v

= R



i



h

ν

i

/k

B

T

exp(h

ν

i

/k

B

T)

− 1

− ln(1 − exp(−hν

i

/k

B

T))

(32)

Vol. 9

COMPUTATIONAL CHEMISTRY

347

S

t

= R



ln



2

πMk

B

T

h

2

3

/2

k

B

T

P

+ 1 + 3/2

(33)

S

r, linear

= R



ln



1

σ

r



T



r

+ 1

(34)

S

r, nonlinear

= R



ln



π

1

/2

σ

r



T

3

/2

(



r

,x



r

,y



r

,z

)

1

/2

+ 3/2

(35)

S

e

= R ln(ω

0

)

(36)

The parameters required to evaluate these expressions are the same as those

used in evaluating the partition functions, as described above.

Finally, by evaluating the derivative of (28) with respect to temperature, it

is possible to derive a relationship between the above thermodynamic quantities
and the empirical Arrhenius expression for reaction rate coefficients (15):

k(T)

= Ae

− Ea/RT

(37)

The frequency factor (A) in this expression is related to the entropy of the

system, as follows.

A

= (c

◦

)

1

− m

e

m

k

B

T

h

exp



S

‡

R

(38)

The Arrhenius activation energy is related to the reaction barrier, as follows.

E

a

= E

0

+ ZPVE + H

‡

+ mRT

(39)

From these expressions it can be seen that the so-called temperature-

independent parameters of the Arrhenius expression are in fact functions of tem-
perature, which is why the Arrhenius expression is only valid over relatively
small temperature ranges. It should also be clear that the ZPVE-corrected barrier
(E

0

+ ZPVE), the enthalpy of activation ( H

‡

), and the Arrhenius activation en-

ergy (E

a

) are only equal to each other at 0 K. At nonzero temperatures, these

quantities are nonequivalent and thus should not be used interchangeably.

Extensions to Transition-State Theory.

Many variants of transition-

state theory have been derived, and a comprehensive review of these recent de-
velopments has been provided by Truhlar and co-workers (68). As already noted,
one of the main extensions to transition-state theory is variational transition-state
theory which, in its simplest form, locates the transition structure as that having
the maximum Gibb’s free energy (rather than internal energy). Other variations
of transition-state theory arise through making different assumptions as to the
statistical distribution of the available energy throughout the different molecular
modes, and through deriving expressions for the partition functions for cases other
than ideal gases. In addition, two simple extensions to the transition-state theory

348

COMPUTATIONAL CHEMISTRY

Vol. 9

equations are the inclusion of corrections for quantum-mechanical tunneling, and
the improved treatment of the low frequency torsional modes. Since these are of
importance in treating certain polymerization-related systems, these are briefly
described below.

Tunneling Corrections.

One of the assumptions inherent in simple

transition-state theory is that motion along the reaction coordinate can be con-
sidered as a classical translation. In general, this assumption is reasonably valid
since the reacting species, being atoms or molecules, are relatively large and thus
their wavelengths are relatively small compared to the barrier width. However,
in the case of hydrogen (and to a lesser extent deuterium) transfer reactions, the
molecular mass of the atom (or ion) being transferred is relatively small, and thus
quantum effects can be very important. Corrections for quantum-mechanical tun-
neling are incorporated into the

κ coefficient of equation 20, and are known as

tunneling coefficients.

There is an enormous variety of expressions available for calculating tunnel-

ing coefficients. The most accurate tunneling methods, such as small curvature
tunneling (69), large curvature tunneling (70), and microcanonical optimized mul-
tidimensional tunneling (71), involve solving the multidimensional Schr¨odinger
equation describing motion of the molecules at every position along the reaction
coordinate. To calculate such tunneling coefficients, specialized software (such as
POLYRATE (72)) is used, and additional quantum-chemical data (such as the
geometries, energies, and frequencies along the entire minimum energy path)
are required. As a result, simpler (and hence less accurate) expressions are of-
ten adopted. These are derived by treating motion along the reaction coordinate
as a function of one variable, the intrinsic reaction coordinate, and hence solv-
ing a one-dimensional Schr¨odinger equation. When this is done using the calcu-
lated energies along the reaction path, the procedure is known as zero-curvature
tunneling (73). However, this procedure still entails the numerical solution of the
Schr¨odinger equation, and hence an additional simplification is also often made.
Instead of using the calculated energies along this path, some assumed functional
form for the potential energy is used instead. This is chosen so that the Schr¨odinger
equation has an analytical solution, and thus a closed expression for the tunneling
coefficient can be derived. The derivation of these simple tunneling coefficients is
described by Bell (74), and the main expressions used in practice are as follows.

The simplest tunneling coefficients are based on the assumption that the

change in energy along the minimum energy path can be described by a truncated
parabola. This functional form provides a good description of the energies near
the transition structure (where tunneling is most significant), but a very poor
description elsewhere. The equation for the tunneling coefficient is given as the
following infinite series, which is frequently truncated after the first few terms.

κ =

1
2

u

‡

sin



1
2

u

‡

 − u

‡

y

− u

‡

/2π



y

2

π − u

‡

−

y

2

4

π − u

‡

+

y

3

6

π − u

‡

− · · ·



where

u

‡

=

hv

‡

kT

and

y

= exp



−

2

πV

u

‡

kT

(40)

Vol. 9

COMPUTATIONAL CHEMISTRY

349

Equation 40 is called a Bell tunneling correction (74), and in this expression

V is the reaction barrier and

ν

‡

is the imaginary frequency (as obtained from

the frequency calculation at the transition structure). By taking the first term of
equation 40, expanding it as an infinite series, and then truncating at an early
order, the (even simpler) Wigner (75) tunneling expression is obtained (74).

κ ≈

1
2

u

‡

sin



1
2

u

‡

 ≈ 1 +

u

2

‡

24

+

u

4

‡

5760

+ · · · ≈ 1 +

u

2

‡

24

(41)

A slightly more realistic description of the change in potential energy along

the minimum energy path is provided by the following Eckart function (76):

V(x)

=

Ay

(1

+ y)

2

+

By

(1

+ y)

where

y

= e

x

/

(42)

To ensure that the function passes through the reactants, products and tran-

sition structures, the parameters A and B are defined as the following functions
of the forward (V

f

) and reverse (V

r

) reaction barriers (where the reaction is taken

in the exothermic direction).

A

= (



V

f

+



V

r

)

2

and

B

= V

f

− V

r

(43)

The remaining parameter

is chosen so as to give the most appropriate

fit to the minimum energy path. If this fit is biased toward the points near the
transition structure (where tunneling is most important), it can be calculated
as the following function of the imaginary frequency

ν

‡

(where c is the speed of

light) (39,41):

=

i

2

πcν

‡



1
8

(B

2

− A

2

)

2

A

3

(44)

The

value obtained from this expression is in mass-weighted coordinates,

which enables the reduced mass to be dropped from the standard (76) Eckart
formulae (41), resulting in the following expression for the permeability of the
reaction barrier G(W) as a function of the energy W:

G(W)

= 1 −

cosh(

α − β) + cosh(δ)

cosh(

α + β) + cosh(δ)

where

α =

4

π

2

h

√

2W

, β =

4

π

2

h



2(W

− B), δ =

4

π

2

h



2A

−

h

2

16

π

2

2

(45)

350

COMPUTATIONAL CHEMISTRY

Vol. 9

The Eckart tunneling correction (

κ) is then obtained by numerically inte-

grating G(W) over a Boltzmann distribution of energies, via the formula (74):

κ =

exp(V

F

/k

B

T)

k

B

T



∞

0

G(W) exp(

−W/k

B

T) dW

(46)

Although this expression requires numerical integration, it does not require

sophisticated software and can be easily implemented on a spreadsheet.

Finally, it is worth making a few comments on the use of the tunneling proce-

dures. Firstly in very general terms, tunneling is important for reactions involving
the transfer of a hydrogen or deuterium atom or ion. The importance of tunneling
can also be established through examination of the parameter u

‡

in equation 40, a

value of u

‡

< 1.5 usually indicating negligible tunneling effects (74). Secondly, in

principle, the more accurate multidimensional tunneling coefficient expressions
should always be used. However, in practice, the more convenient one-dimensional
expressions are often adopted. Of these expressions, the Eckart tunneling coef-
ficient is significantly more accurate and should be preferred. For example, the
small curvature tunneling method gives a tunneling coefficient (

κ) of approxi-

mately 10

2

at 300 K, for the hydrogen abstraction reaction between

•

NH

2

and

C

2

H

6

(77). At the same level of theory, the corresponding

κ values for the Wigner,

Bell, and Eckart corrections are approximately 5, 10

3

, and 10

2

respectively, and

hence only the Eckart method yields a tunneling coefficient of the right order of
magnitude for this (typical) reaction. The success of the Eckart tunneling method
has also been noted by Duncan and co-workers (78), who rationalized it in terms
of a favorable cancelation of errors.

Treatment of Low Frequency Torsional Modes.

In the vibrational parti-

tion function (eq. 24), all modes are treated under the harmonic oscillator approx-
imation. That is, it is assumed that the potential field associated with their dis-
tortion from the equilibrium geometry is a parabolic well, as in a vibrating spring
(see Fig. 10a). This is a reasonable assumption for bond stretching motions, but
not for hindered internal rotations (see Fig. 10b). For high frequency modes (

ν >

200–300 cm

− 1

), the contribution of these motions to the overall partition function

is negligible at room temperature (ie Q

vib

,i

≈ 1) and thus the error incurred in

treating these modes as harmonic oscillators is not significant. However, for the
low frequency torsional modes, these errors can be significant and a more rigorous
treatment is often necessary, and this is especially the case for the reactions of
relevance to free-radical polymerization (3–5,79).

The simplest method for treating the hindered internal rotations is to regard

them as one-dimensional rigid rotors. An appropriate rotation angle

θ is identified,

and then the potential V(

θ) is calculated as a function of this angle (eg from 0 to

360

◦

in steps of 10

◦

) via quantum chemistry. In general, it is recommended that

these potentials be obtained as relaxed scans; that is, in calculating the energy
at a specific angle, all geometric parameters other than the rotational angle are
optimized (79). If the rotational potential can be described as a simple cosine
function, the enthalpy and entropy associated with the mode can be obtained
directly from the tables of Pitzer and co-workers (80). In order to use these tables,

Vol. 9

COMPUTATIONAL CHEMISTRY

351

Fig. 10.

Typical potentials associated with (a) a harmonic oscillator and (b) a hindered

internal rotation.

one calculates two dimensionless quantities:

x

=

V

k

B

T

and

y

=

σ

int

h



8

π

3

I

m

k

B

T

(47)

In these equations, V is the barrier to rotation,

σ

int

is the symmetry number

associated with the rotation (which counts the number of equivalent minima in
the potential energy curve), and I

m

is the reduced moment of inertia associated

with the rotation. This latter parameter is given by the following formula:

I

m

= A

m



1

−



i

= x,y,z

A

m

λ

2
mi



I

i

(48)

In this equation, A

m

is the moment of inertia of the torsional coordinate

itself, I

i

is the principal moment of inertia of the whole molecule about axis i, and

352

COMPUTATIONAL CHEMISTRY

Vol. 9

λ

mi

is the direction cosine between the axis of the top and the principal axis of

the whole molecule. More information on the calculation of reduced moments of
inertia can be found in Reference (81). When the rotational potential cannot be
fitted with a simple cosine function (as in Fig. 10b), the partition function (and
hence the enthalpy and entropy) is obtained instead by (numerically) solving the
one-dimensional Schr¨odinger equation.

−

h

2

8

πI

m

∂

2



∂θ

2

+ V(θ) = ε

(49)

This yields the energy levels of the system (

ε

i

), which are then used to eval-

uate the partition function via the following equation.

Q

int rot

=

1

σ

int



i

exp



−

ε

i

k

B

T

(50)

Having obtained the partition function (or equivalently, the enthalpy and

entropy) associated with a low frequency torsional mode, this is used in place of
the corresponding harmonic oscillator contribution for that mode.

The above treatment of hindered rotors assumes that a given mode can be

approximated as a one-dimensional rigid rotor, and studies for small systems
have shown that this is generally a reasonable assumption in those cases (82).
However, for larger molecules, the various motions become increasingly coupled,
and a (considerably more complex) multidimensional treatment may be needed in
those cases. When coupling is significant, the use of a one-dimensional hindered
rotor model may actually introduce more error than the (fully decoupled) harmonic
oscillator treatment. Hence, in these cases, the one-dimensional hindered rotor
model should be used cautiously.

Software

There are a large number of software packages available for performing compu-
tational chemistry calculations. Some of the programs available include ACES II
(83), GAMESS (84), GAUSSIAN (43), MOLPRO (85), and QCHEM (86). Other
programs, such as POLYRATE, (72), have been designed to use the output of
quantum-chemistry programs to calculate reaction rates and tunneling coeffi-
cients. Whereas computational chemistry software has traditionally been oper-
ated on large supercomputers, versions for desktop personal computers are in-
creasingly becoming available. In addition, programs for visualizing the output
of computational chemistry calculations such as Spartan (87), Molden (88), CS
Chem3D (89), Molecule (90), Jmol (91), Gauss View (43), and MacMolPlt (92) are
also available. These programs allow one to visualize the geometry and electronic
structure of the resulting molecule, and animate its vibrational frequencies. Many
of these programs also have built-in computation engines. Computational chem-
istry is thus increasingly becoming accessible to the nonspecialist user, which
brings with it its own problems (see also M

OLECULAR

M

ODELING

).

Vol. 9

COMPUTATIONAL CHEMISTRY

353

Accuracy and Applicability of Theoretical Procedures

By solving the Schr¨odinger equation exactly, quantum chemistry can, in princi-
ple, provide accurate electronic-structure data. However, in practice, approximate
numerical methods must be adopted, and these can introduce error to the calcula-
tions. As we have already seen, an enormous number of approximate methods are
available, and these range from the accurate but computationally expensive to the
cheap but potentially nasty. Furthermore, the performance of a particular method
varies considerably depending upon the chemical system and the property being
calculated. It is therefore very important that computational chemistry studies
are accompanied by rigorous assessments of theoretical procedures. In such “cali-
bration” studies, prototypical systems are calculated at a range of levels of theory,
and the results are compared both internally (against the highest level proce-
dures) and externally (against reliable experimental data) in order to identify
those methods which offer the best compromise between accuracy and expense.

In the present section, the main conclusions from recent assessment stud-

ies for the reactions of importance to free-radical polymerization are outlined.
In presenting such studies, it must be acknowledged that, with continuing rapid
increases in computer power, some of these results will soon be outdated. In par-
ticular, as computer power increases, the need to rely upon lower levels of theory
for large polymer-related systems will diminish. Instead, the higher levels of the-
ory outlined below will be able to be used routinely. Nonetheless, with increasing
computer power, the temptation to apply existing levels of theory to yet larger
systems will no doubt ensure that the main conclusions of these studies retain
some relevance into the near future.

Radical Addition to C C Bonds.

Radical addition to C C bonds are of

importance for free-radical polymerization as this reaction forms the propagation
step, and thus influences the reaction rate and molecular weight distribution in
both conventional and controlled free-radical polymerization, and the copolymer
composition and sequence distribution in free-radical copolymerization. Numer-
ous studies have examined the applicability of high level theoretical methods for
studying radical addition to C C bonds in small radical systems (32,33,37,93,94).
The most recent study (37) included W1 barriers and enthalpies, and geometries
and frequencies at the CCSD(T)/6-311G(d,p) level of theory, and is the highest level
study to date. The main conclusions from this study, and (where still relevant) the
previous lower level studies, are outlined below.

Geometry optimizations are generally not very sensitive to the level of the-

ory, with even the low cost HF/6-31G(d) and B3-LYP/6-31G(d) methods providing
reasonable approximations to the higher level calculations (37). In the latter case,
there is a small error arising from the tendency of B3-LYP to overestimate the
forming bond length in the transition structures, and this can be reduced us-
ing an IRCmax technique (94). Alternatively, the error in the B3-LYP transition
structures is also reduced when the larger 6-311

+ G(3df,2p) basis set is used

(37). In addition, the UMP2 method should generally be avoided for these reac-
tions, as they are subject to spin-contamination problems (32,33,37). Frequency
calculations are also relatively insensitive to the level of theory, especially when
the frequencies are scaled by their appropriate scale factors. (Scale-factors for the
most commonly used levels of theory may be found in Reference (95). In particular,

354

COMPUTATIONAL CHEMISTRY

Vol. 9

the B3-LYP/6-31G(d) level of theory provides excellent performance for frequency
factors, temperature corrections, and zero-point vibrational energy calculations,
and would be a suitable low cost method for studying larger systems (37).

Barriers and enthalpies are very sensitive to the level of theory. Where possi-

ble, high level composite procedures should be used for the prediction of absolute
reaction barriers and enthalpies, and of these methods the “RAD” variants of G3
provide the best approximations to the higher level Wn methods (when the lat-
ter cannot be afforded) (37). It should also be noted that the (empirically based)
spin-correction term in the CBS-type methods appears to be introducing a consid-
erable error to the predicted reaction barriers for these reactions and, until this
is revised, these methods should perhaps be avoided for these reactions (37).

When composite methods cannot be afforded, the use of RMP2/6-

311

+ G(3df,2p) single points provides reasonable absolute values and excellent

relative values for the barriers and enthalpies of these reactions (37). In contrast,
the hybrid DFT methods such as B3-LYP and MPW1K show considerable error in
the reaction enthalpies, even when applied with large basis sets. However, they do
provide reasonable addition barriers, owing to the cancelation of errors in the early
transition structures for these reactions (37). Interestingly, for the closely related
radical addition to C C bonds, the situation is reversed and the B3-LYP methods
perform well and the RMP2 methods perform poorly (37), and this highlights the
importance of performing assessment studies before tackling new chemical prob-
lems. Finally, it should be stressed that semiempirical methods do not provide an
adequate description of the barriers and enthalpies in these reactions (32).

For rate coefficients, the importance of treating the low frequency torsional

modes in radical addition reactions as hindered internal rotations has been in-
vestigated in a number of assessment studies (37,79,93). For small systems such
as methyl addition to ethylene and propylene, the errors are relatively minor
(less than a factor of 2) (37). However, for reactions of substituted radicals (such
as n-alkyl radicals (93) and the ethyl benzyl radical (79)), the errors are some-
what larger (as much as a factor of 6), as there are additional low frequency
torsional modes to consider. Nonetheless, the errors are still relatively small, and
the harmonic oscillator approximation might be expected to provide reasonable
“order-of-magnitude” estimates of rate coefficients.

Radical Addition to C S Bonds.

Radical addition to C S bonds, and

the reverse

β-scission reaction, forms the key addition and fragmentation steps of

the RAFT polymerization process (96). Ab initio calculations have a role to play in
elucidating the effects of substituents on this process, and in providing an under-
standing of the causes of rate retardation (6,7). A detailed assessment of theoreti-
cal procedures has been recently carried out for this class of reactions (36), and the
main conclusions are similar to those for addition to C C bonds (37), as outlined
above. In general, low levels of theory, such as B3-LYP/6-31G(d), are suitable for ge-
ometry optimizations and frequency calculations, provided an IRCmax procedure
is used to correct the transition structures. However, high level composite methods
are required to obtain reliable absolute barriers and enthalpies, though reason-
able relative quantities can be obtained at the RMP2/6-311

+ G(3df,2p) level. As

in the case of addition to C C bonds, the spin correction term in the CBS-type
methods appears to require adjustment, and the RAD variants of G3 should be
preferred when the higher level Wn calculations are impractical.

Vol. 9

COMPUTATIONAL CHEMISTRY

355

Hydrogen Abstraction.

Hydrogen abstraction reactions are important

chain-transfer processes in free-radical polymerization. In particular, hydrogen
abstraction by the propagating polymer radical from transfer agents (such as
thiols), monomer, dead polymer, or itself (ie intramolecular abstraction) can affect
the molecular weight distribution, the chemical structure of the chain ends, and
the degree of branching in the polymer. The accuracy of computational procedures
for studying hydrogen abstraction reactions has received considerable attention,
and the results of some of the most recent and extensive studies (59,94,97–99) are
summarized below.

Geometry optimizations are relatively insensitive to the level of theory; how-

ever, there are some important exceptions. In particular, the HF and MP2 meth-
ods should be avoided for spin-contaminated systems (99). Moreover, the B3-LYP
method does not describe the transition structures very well for a number of hy-
drogen abstraction reactions (59,97). However, improved performance is obtained
using newer hybrid DFT methods such as MPW1K (59) and KMLYP (97), and
these methods are suitable as low cost methods, when high level procedures can-
not be afforded.

Barriers and enthalpies are more sensitive to the level of the theory, and,

where possible, high level composite procedures should be used. In particular, the
“RAD” variants of G3 provide an excellent approximation to the higher level Wn
methods, and would provide an excellent benchmark level of theory when the lat-
ter could not be afforded (99). As in the case of the addition reactions, the spin
contamination correction term in the CBS-type methods appears to be introduc-
ing a systematic error to the predicted reaction barriers and enthalpies and, un-
til this is revised, this method should perhaps be avoided for spin-contaminated
reactions (99). When composite methods cannot be afforded, methods such as
RMP2, MPW1K, or KMLYP have been shown to provide good agreement with the
high level values (59,97,99), with a procedure such as MPW1K/6-311

+ G(3df,2p)

providing the best overall performance. By contrast, the popular B3-LYP method
performs particularly poorly for reaction barriers and enthalpies (59,94,97–99),
and should thus be generally avoided for abstraction reactions. Interestingly, it
has been noted that the errors in B3-LYP increase with the increasing polarity of
the reactants (98), which suggests that assessment studies based entirely on rela-
tively nonpolar reactions (such as CH

3

•

+ CH

4

) may lead to the wrong conclusions.

As noted in the previous section, tunneling is significant in hydrogen abstrac-

tion reactions, and hence accurate quantum-chemical studies of these systems
require the calculation of tunneling coefficients. The accuracy of tunneling coeffi-
cients is profoundly affected by both the tunneling method and the level of theory
at which it is applied. A systematic comparison of the various tunneling methods
for the hydrogen abstraction reactions of relevance to free-radical polymerization
does not appear to have been performed. However, in the example provided in the
previous section, it was seen that the Eckart method was capable of providing
the tunneling coefficients of the right order of magnitude (when compared with
the more accurate multidimensional methods), while the Wigner and Bell meth-
ods respectively underestimated and overestimated the tunneling coefficients by
an order of magnitude. Hence, when multidimensional tunneling methods are
not convenient, the Eckart tunneling method should be preferred as the best
low cost method. An assessment of the effects of level of theory on the tunneling

356

COMPUTATIONAL CHEMISTRY

Vol. 9

coefficients, as calculated using this Eckart method, has recently been published
(41). It was found that errors in the imaginary frequency at the HF level (with
a range of basis sets) leads to errors in the calculated tunneling coefficients of
several orders of magnitude (compared to high level CCSD(T)/6-311G(d,p) calcu-
lations). The B3-LYP and MP2 methods performed significantly better, showing
errors of a factor of 2–3. However, even better performance could be obtained
by correcting the HF values to the CCSD(T)/6-311G(d,p) level via an IRCmax
procedure.

Applicability of Chemical Models

Assuming high levels of theory are used, quantum-chemical calculations might be
expected to yield very accurate values of the rate coefficients for the specific chem-
ical system being studied. With current available computing power, this would in
all likelihood consist of a small model reaction in the gas phase. If, for exam-
ple, this information is then to be used to deduce something about solution-phase
polymerization kinetics, the effects of the solvent and (in most cases) the effects of
chain length need to be considered. Unfortunately, the treatment of these effects
at a high level of theory is not generally feasible with current available computing
power, and hence the neglect of these effects (or their treatment at a crude level
of theory) remains a potential source of error in quantum-chemical calculations.
In this section, these additional sources of error are briefly discussed.

Solvent Effects.

The presence of solvent molecules may affect the poly-

merization process in a variety of ways (100). For example, if polar interactions
are significant in the transition structure of the reaction, the presence of a high
dielectric constant solvent may stabilize the transition structure and lower the re-
action barrier. Solvents may also affect the reactivity of the reacting radicals, and
even the mechanism of the addition or transfer reaction, through some specific
interaction such as hydrogen bonding or complex formation. In addition, the pref-
erential sorption of the monomer or solvent around the reacting polymer radical
may lead to the effective concentrations available for reaction being different to
those in the bulk solution, resulting in a difference between the observed and pre-
dicted rate coefficients. Solvent effects such as these result in the experimentally
measured rate coefficients for a free-radical polymerization varying according to
the solvent type.

Over and above these system-specific solvent effects, there is a more general

“entropically based” difference between the rate coefficients for gas-phase and
solution-phase systems. Whereas in the gas-phase an isolated molecule might be
expected to have translational, rotational, and vibrational degrees of freedom, in
the solvent phase the translational and rotational degrees of freedom are effec-
tively “lost” in collisions with the solvent molecules. In their place, it is necessary
to consider additional vibrational degrees of freedom involving a solute-solvent
“supermolecule” (101). Since the vibrational, translational, and rotational modes
make different quantitative contributions to the enthalpy and entropy of acti-
vation, significant differences might be expected between the gas and solution
phases. For bimolecular reactions this effect can be considerable, because the main
contribution to the entropy of activation is the six rotational and translational

Vol. 9

COMPUTATIONAL CHEMISTRY

357

degrees of freedom in the reactants, which are converted to internal vibrations
in the transition structure. In contrast, for unimolecular reactions, the entropi-
cally based gas-phase/solution-phase difference is generally much smaller, as the
rotational and translational modes are similar for the transition structure and re-
actant molecule, and thus their contribution largely cancels from the reaction rate.

The treatment of solvent effects varies according to their origin. The influ-

ence of the dielectric constant on polar reactions can be dealt with routinely using
various continuum models (102), implemented using standard computational soft-
ware, such as GAUSSIAN (43). However, it should be stressed that these models
do not account for the entropically based gas-phase/solvent-phase difference, nor
do they deal with direct solvent interactions in the transition structure. When
direct interactions involving the solvent are important, it is necessary to include
solvent molecules in the quantum-chemical calculation. In theory, one should in-
clude many hundreds of solvent molecules but in practice one includes a small
number of molecules, and combines this with a continuum model (102). However,
even with this simplification, the additional solvent molecules increase the com-
putational cost of the calculations, and it is not currently feasible to apply these
methods (at any reasonable level of theory) for polymer-related systems. Even
when additional solvent molecules are included in the ab initio calculations, vari-
ous extensions to transition-state theory are required in order to model the rates of
solution-phase reactions (68,101). Unfortunately, the existing models are compli-
cated to use and require additional parameters which are not readily accessible
for polymerization-related systems. The development of simplified yet accurate
models for dealing with solvent effects is an on-going field of research.

While strategies for calculating solution-phase rate coefficients exist, with

the current available computing power these methods are not generally feasible
for polymer-related systems. Instead, the following practical guidelines for dealing
with solvent effects are suggested. Firstly, when the solvent participates directly
in the reaction, the inclusion of the interacting solvent molecule in the gas-phase
calculation is essential for gaining a mechanistic understanding of the reaction.
Secondly, when polar interactions are expected to be important, the use of a con-
tinuum model is recommended, especially if the results are to be used to interpret
the polymerization process in polar solvents. Thirdly, for bimolecular reactions,
if the a priori prediction of absolute rate coefficients is required, a consideration
of the entropically based gas-phase/solution-phase difference is necessary. This
“entropic” solvent effect might be estimated by comparing corresponding exper-
imental solution- and gas-phase rate coefficients for that class of reaction. For
example, solution-phase experimental values for radical addition reactions gen-
erally exceed the corresponding gas-phase values by approximately one order of
magnitude (103). One might also benchmark the gas-phase calculations by cal-
culating the rate coefficient for a similar reaction, and comparing the calculated
result with reliable experimental data. Fourthly, provided that specific interac-
tions are not important, one might expect that solvent effects should largely cancel
from relative rate coefficients, and hence the gas-phase values should generally
be suitable for studying substituent effects and solving mechanistic problems. Fi-
nally, when specific interactions are important, simple gas-phase calculations are
still useful, as they can provide complementary information about the underlying
influences on the mechanism in the absence of the solvent.

358

COMPUTATIONAL CHEMISTRY

Vol. 9

Chain-Length Effects.

The other simplification that is necessary in order

to use high level ab initio calculations on polymer-related systems is to approxi-
mate the propagating polymer radical (which may be hundreds or thousands of
units long) as a short-chain alkyl radical. Provided that the reaction is chemi-
cally controlled, this is not an unreasonable assumption. In chemical terms, the
effects of substituents decrease dramatically as they are located at positions that
are increasingly remote from the reaction center. For example, the terminal and
penultimate units of a propagating polymer radical are known to affect its re-
activity and selectivity in the propagation reaction (104); however, substituent
effects beyond the penultimate position are rarely invoked in copolymerization
models. In order to include the most important substituent effects, it is generally
recommended that propagating radicals be represented as

γ -substituted propyl

radicals (as a minimum chain length). For some systems this is not currently
possible without resorting to a low (and thus inaccurate) level of theory. In those
cases, the possible influence of penultimate unit effects must be taken into account
when interpreting the results of the calculations.

The entropic influence of the chain length on the reaction rates extends

slightly beyond the penultimate unit. For example, Deady and co-workers (105)
showed experimentally that there was a chain-length effect on the propagation
rate coefficient of styrene, which converged at the tetramer stage (ie an octyl
radical). Heuts and co-workers (3) have explored this chain-length dependence
theoretically, and suggested that it arises predominantly in the translational and
rotational partition functions. More specifically, they suggest that there is a small
effect of mass that can be modeled by including an unrealistically heavy isotope
of hydrogen at the remote chain end. For example, in the propagation of ethylene,
a model such as X (CH

2

)

n

CH

2

•

could be used, and in this model X is set as a hy-

drogen atom that happens to have a molecular mass of 9999 amu! They also noted
that there is an effect of chain length on the rotational entropy (and especially the
hindered internal rotations), which required the more subtle modeling strategy
of using slightly longer alkyl chains (ie n

> 1). Nonetheless, using their “heavy

hydrogen” approach, their calculated frequency factors converged to within a fac-
tor of 2 of the long chain limit at even the propyl radical stage (ie n

= 2). More

recently it was shown that the consideration of the propagating radical as a sub-
stituted hexyl radical (without a heavy hydrogen at the remote chain end) was
also sufficient to reproduce experimental values for the frequency factors of prop-
agation reactions (106). For short-chain branching reactions, it has been shown
that inclusion of just one methyl group beyond the reaction center is sufficient
for modeling the long-chain reactions, provided that the additional methyl group
is substituted with a heavy (ie 9999 amu) hydrogen atom (5). Thus, in general,
it appears that small model alkyl radicals are capable of providing a reasonable
description of polymeric radicals in chemically controlled reactions.

Finally, it is worth noting that the chain-length effects on the propagation

steps amount to approximately an order of magnitude difference between the first
propagation step and the long chain limit, with the small radical additions having
faster rate coefficients. This chain-length error is of the same magnitude and acts
in the opposite direction to the gas-phase/solvent-phase difference in bimolecu-
lar reactions, and hence substantial cancelation of error might be expected in
these cases. Indeed an (unpublished) high level G3(MP2)-RAD calculation of the

Vol. 9

COMPUTATIONAL CHEMISTRY

359

propagation rate coefficient for methyl acrylate at 298 K produced the value of 2.0
× 10

4

L

· mol

− 1

· s

− 1

, which is in remarkable agreement with the corresponding ex-

perimental value (also 2.0

× 10

4

L

· mol

− 1

· s

− 1

at ambient pressure) (107), despite

the fact that both the medium and chain-length effects were ignored in the calcula-
tion. Hence, careful efforts to correct for chain-length effects but not solvent effects
(or vice versa) may actually introduce greater errors to calculated rate coefficients.

Applications of Quantum Chemistry in Free-Radical Polymerization

Quantum chemistry provides a powerful tool for studying kinetic and mechanistic
problems in free-radical polymerization. Provided a high level of theory is used,
ab initio calculations can provide direct access to accurate values of the barriers,
enthalpies, and rates of the individual reactions in the process, and also provide
useful related information (such as transition structures and radical stabilization
energies) to help in understanding the reaction mechanism. In the following, some
of the applications of quantum chemistry are outlined. This is not intended to
be a review of the main contributions to this field, nor is it intended to provide
a theoretical account of reactivity in free-radical polymerization (108). Instead,
some of the types of problems that quantum chemistry can tackle are described,
with a view to highlighting the potential of quantum-chemical calculations as a
tool for studying free-radical polymerization (see R

ADICAL

P

OLYMERIZATION

).

A Priori Prediction of Absolute Rate Coefficients.

The a priori pre-

diction of accurate absolute rate coefficients is perhaps the most demanding task
in computational chemistry. For example, high level ab initio calculations (at the
G3(MP2)-RAD level as a minimum) are required for the calculation of accurate
reaction barriers (and enthalpies) in radical addition reactions and, with current
available computing power, these can be performed routinely on systems of up
to 12–14 non-hydrogen atoms. This allows for the most common polymerization
substituents to be included, and often allows for penultimate unit effects to be
taken into account. However, it does not allow for the accurate treatment of sol-
vent or chain-length effects. Of course, with continuing rapid increases in com-
puter power, these problems should soon be overcome. In addition, as was noted
above, for a bimolecular reaction, the chain-length effects and solvent effects act
in opposite directions and may in fact largely cancel, leading to surprisingly accu-
rate results. Certainly, with current available computing power, the prediction of
rate coefficients to within 1–2 orders of magnitude (or better) is possible for most
chemically controlled polymerization-related reactions, and further increases in
accuracy should be attainable in the near future.

The prediction of absolute rate coefficients via quantum chemistry is partic-

ularly useful when direct experimental measurements are not possible. A good
example of this is in the RAFT polymerization process (96), in which the experi-
mentally observable quantities (such as the overall polymerization rate, the over-
all molecular weight distribution, and the concentrations of the various species)
are a complicated function of the rates of the various individual reactions. In order
to measure the rates of these individual steps, it is necessary to relate their rate co-
efficients to the observable quantities via some kinetic-model-based assumptions.
Depending upon the assumptions, enormous discrepancies arise in the estimated

360

COMPUTATIONAL CHEMISTRY

Vol. 9

rate coefficients. For example, alternative experimental measurements for the
fragmentation rate in cumyl dithiobenzoate mediated polymerization of styrene
at 60

◦

C differ by six orders of magnitude (109,110), and this difference arises, at

least in part, in the model-based assumptions of the alternative experimental tech-
niques (For a discussion of problems measuring fragmentation rate coefficients,
the reader is referred to Reference 111). Recent ab initio calculations of the frag-
mentation rate coefficient for a model of this system were able to provide direct
evidence in support of one measurement (and hence one set of assumptions) over
the other, and in doing so provide an insight into the causes of rate retardation in
the RAFT process (6).

As accurate calculations become routine, there will be many other applica-

tions for the a priori prediction of absolute rate coefficients. The calculation of
propagation rate coefficients (which can be measured experimentally via pulsed
laser polymerization or PLP (112)) is currently used to benchmark theoretical pro-
cedures, but accurate calculations may also be helpful as a stand-alone technique
for toxic and hazardous monomers, and in other cases where PLP is difficult (due
to problems such as the monomer absorbing at the wavelength of the laser or reac-
tions such as chain transfer broadening the molecular weight distribution). It has
also been noted that ab initio calculations may provide the best means of study-
ing long- and short-chain branching in free-radical polymerization (5). Quantum
chemistry will also be helpful in extracting the rate coefficients of the individual
reactions in other complicated reaction schemes, such as free-radical copolymer-
ization, and the various types of controlled radical polymerization systems.

Prediction of Relative Rate Coefficients: Discriminating Models and

Mechanisms.

High level ab initio calculations can already predict relative rate

coefficients with remarkable accuracy. This is because accurate relative values
of quantities such as barriers and enthalpies can generally be obtained at lower
levels of theory than corresponding absolute values, because of substantial cancel-
lation of error. In addition, one might generally expect that chain-length effects
and solvent effects should be reasonably consistent for a series of similar reac-
tions, and thus cancel from the comparative values. The prediction of relative
rate coefficients is important for modeling and hence optimizing various aspects
of the polymerization process, and some of these applications are outlined below.
The prediction of relative rate coefficients is also important for understanding the
effects of substituents on the various individual reactions, and these applications
are discussed in the following section.

Copolymerization.

In free-radical copolymerization (qv), the relative rates

of addition of the various types of radical to the alternative monomers are called
reactivity ratios, and are the key parameters governing the composition and se-
quence distribution of the resulting copolymer. Experimentally, these parameters
are “measured” for a given system by treating them as adjustable parameters in
a least-squares fit of some assumed kinetic model to experimental values of the
composition, sequence distribution, and/or propagation rate coefficients. However,
there are two important problems with this approach (104). Firstly, it is not always
(if ever) clear which copolymerization model (eg the terminal model, the implicit
or explicit penultimate model, etc) is appropriate for a given system. If a phys-
ically incorrect model is chosen, then the estimated parameters will lack their
assumed physical meaning, and will thus be unsuitable for mechanistic studies,

Vol. 9

COMPUTATIONAL CHEMISTRY

361

or for predicting other copolymerization properties. Secondly, even if the correct
model is chosen, there are often more adjustable parameters than are numerically
needed to fit the data, and this can result in large and highly correlated uncertain-
ties in the estimated parameters. The problems in the experimental estimation of
reactivity ratios are discussed in more detail in the article on copolymerization
(qv).

Quantum chemistry is able to address these problems, as it allows the re-

activity ratios to be calculated directly from the rate coefficients for the vari-
ous types of reaction, without making kinetic-model-based assumptions. Further-
more, quantum chemistry can assist in determining which kinetic model is most
appropriate for a given system. For example, by measuring the effects of sub-
stituents at the penultimate position of the propagating radical, quantum chem-
istry can be used to determine whether a terminal or penultimate model is ap-
propriate for a given system. In this respect, quantum chemistry has already
made an important contribution to this field by providing evidence for penulti-
mate unit effects in the barriers (1,106) and frequency (2,106) factors of a variety
of different systems. Such calculations can thus simplify the model discrimina-
tion process by ruling out either the terminal- or penultimate-based models for a
given system, and by providing estimates for the main model parameters (ie the
reactivity ratios). In addition, quantum chemistry can provide a simple and cost-
effective method for screening the reactivity ratios for a wide variety of monomer
pairs (including hazardous and yet-to-be synthesized monomers). This may al-
low the identification of monomer pairs with suitable reactivity ratios for a given
application.

Chain Transfer.

The relative rate of abstraction to propagation for a spe-

cific propagating radical is known as the chain-transfer constant, and is a key
parameter governing the molecular weight of the resulting polymer. High level
ab initio calculations allow chain-transfer constants to be predicted for a given
propagating radical and transfer agent, and thus provide an effective means of
screening large numbers of transfer agents. This in turn may assist in the selec-
tion of suitable transfer agents for a given polymerization system, and can also
aid in our understanding of reactivity in the chain-transfer processes. Already
ab initio calculations have been used to study intramolecular chain-transfer pro-
cesses in free-radical polymerization (5,106), in order to provide estimates of the
short-chain branching ratios. Another study investigated the kinetics and ther-
modynamics of the hydrogen abstraction by and from the monomer in ethylene
polymerization (113), and demonstrated that abstraction from the monomer was
the kinetically (but not thermodynamically) preferred process. There is also an
enormous body of literature concerning ab initio calculations of hydrogen abstrac-
tion reactions in general, and in applications of relevance to other fields such as
biochemistry (114,119).

Controlled Radical Polymerization.

In recent years, the field of free-radical

polymerization has been revolutionalized by the development of techniques for
controlling the molecular weight and architecture of the resulting polymer, in-
cluding nitroxide-mediated polymerization (NMP) (120), atom-transfer polymer-
ization (ATRP) (121), and reversible addition fragmentation chain transfer (RAFT)
polymerization (96) (see L

IVING

R

ADICAL

P

OLYMERIZATION

). In order to control

polymerization, these processes aim to minimize the influence of bimolecular

362

COMPUTATIONAL CHEMISTRY

Vol. 9

termination processes via a delicate balance of the relative rates of two or more
competing reactions. Quantum chemistry can assist in optimizing these processes
by providing reliable values for the rate coefficients of the competing reactions,
across a wide range of systems. More generally, quantum chemistry can help to
provide a detailed understanding of the kinetics and mechanisms of the individual
reactions, which can in turn allow for the design of improved control agents.

Quantum chemistry has already been used to study the RAFT process (6,7).

As noted above, high level ab initio molecular orbital calculations have been used
to obtain direct measurements of the rates of addition and fragmentation in model
dithiobenozate-mediated systems, which in turn provided evidence that slow frag-
mentation is responsible for rate retardation in these systems (6). More recently,
ab initio calculations have revealed a new and unexpected side reaction (

β-scission

of the alkoxy group) in certain xanthate-mediated systems, which may provide an
explanation for the experimentally observed inhibition in those systems (7). The
ab initio calculations also indicated that fragmentation was considerably faster
in the xanthate-mediated systems (compared with other dithioester systems), be-
cause of the stabilizing influence of the alkoxy group on the thiocarbonyl product
(7). In addition to these ab initio calculations, semiempirical methods have been
used to survey the transfer enthalpies for a series of RAFT agents in styrene,
methyl methacrylate, and butyl acrylate polymerization (122). Despite the low
level of theory used, and the limited accuracy of the quantitative predictions at
this level, the calculations were nonetheless shown to have some qualitative value
in determining which transfer agents would show the best molecular weight con-
trol in RAFT polymerization.

Understanding Reactivity via the Curve-Crossing Model.

We have

already seen that quantum chemistry can be used to calculate the absolute
(and hence relative) rate coefficients for the individual reactions in free-radical
polymerization, and is well suited to compiling systematic surveys of the rate
coefficients for homologous series of reactions. When such calculations are
combined with qualitative theoretical models, quantum chemistry can help to
provide an understanding of trends in reactivity. It is beyond the scope of this
chapter to provide an account of the main qualitative models and concepts
in general theoretical chemistry. However, it is worth introducing one such
model—the curve-crossing model (123)—as this model has been used extensively
to provide a qualitative rationalization of the effects of substituents in some of the
key radical reactions that occur in free-radical polymerization, including radical
addition to alkenes (103) (ie the propagation step), and hydrogen abstraction
(115,116,124) (ie, chain-transfer processes). The main qualitative features of the
model and its application to free-radical polymerization are described below; for
a more detailed description, the reader is referred to the original references 123,
and also the excellent book by Pross (125).

The curve-crossing model (also referred to as the valence-bond state corre-

lation model, the configuration mixing model, or the state correlation diagram)
was developed by Pross and Shaik (123) as a unifying theoretical framework for
explaining barrier formation in chemical reactions. It is largely based on valence
bond (VB) theory (126,127), but also incorporates insights from qualitative molec-
ular orbital theory (128). To understand the curve-crossing model, it is helpful to
think of a chemical reaction as being composed of a rearrangement of electrons,

Vol. 9

COMPUTATIONAL CHEMISTRY

363

accompanied by a rearrangement of nuclei (ie a geometric rearrangement). We can
then imagine holding the arrangement of electrons constant in its initial configu-
ration (which we call the reactant VB configuration), and examining how the en-
ergy changes as a function of the geometry. Likewise, we could hold the electronic
configuration constant in its final form (the product VB configuration), and again
examine the variation in energy as a function of the geometry. If these two curves
(energy vs geometry) are plotted, we form a “state correlation diagram.” The over-
all energy profile for the reaction, which is also plotted, is formed by the resonance
interaction between the reactant and product configurations (and any other low
lying configurations). State correlation diagrams allow for a qualitative explana-
tion for how the overall energy profile of the reaction arises, and can then be used
to provide a graphical illustration of how variations in the relative energies of the
alternative VB configurations affect the barrier height. This in turn allows us to
rationalize the effects of substituents on reaction barriers, and to predict when
simple qualitative rules (such as the Evans–Polanyi rule (129)) should break down.

This procedure is best illustrated by way of an example, such as the case of

radical addition to alkenes. For this type of reaction, the principal VB configura-
tions that may contribute to the ground-state wave function are the four lowest
doublet configurations of the three-electron–three-center system formed by the
initially unpaired electron at the radical carbon (R) and the electron pair of the
attacked

π bond in the alkene (A) (103).

The first configuration (RA) corresponds to the arrangement of electrons in

the reactants, the second (RA

3

) to that in the products, and the others (R

+

A

−

and

R

−

A

+

) to possible charge-transfer configurations. The state correlation diagram

showing (qualitatively) how the energies of these configurations vary as a function
of the reaction coordinate is provided in Figure 11 (103). Let us now examine how
this plot could be made for a specific system.

The anchor points on these diagrams are generally accessible from quantum-

chemical calculations. For example, the energy difference between the RA config-
uration at the reactant geometry, and the RA

3

configuration at the product geom-

etry, is simply the energy change of the reaction. The energy difference between
the RA and RA

3

configurations at the reactant geometry is simply the vertical

singlet-triplet gap of the alkene. At the product geometry, the RA

− RA

3

energy

difference is also an excitation energy, but this time between the ground state of
the doublet product and an excited doublet state (In fact, studies of radical addi-
tion to alkenes have ignored the influence of this excitation gap, without detriment
to the predictive value of the results (103). This is probably due to the fact that
the transition structure is very early in these reactions. However, in studies of
hydrogen abstraction reactions, the RA and RA

3

gaps at both the reactant and

product geometries can be important. These are calculated as the singlet-triplet
gap of the closed-shell substrate in each case.) The charge-transfer configurations
can be anchored at the reactant geometry, where they are given as the energy for

364

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 11.

State correlation diagram for radical addition to alkenes showing the variation

in energy of the reactant (RA), the product (RA

3

), and the charge-transfer configurations

(R

+

A

−

and R

−

A

+

) as a function of the reaction coordinate. The dashed line represents the

overall energy profile of the reaction.

complete charge transfer between the isolated reactants. For example, the energy
difference between the R

+

A

−

and the RA configuration at the reactant geometry

would be given as the energy change of the reaction R

+ A → R

+

+ A

−

. It can be

seen that the energy change of this reaction is simply the difference between the
ionization energy (R

→ R

+

+ e−) of the donor species and electron affinity (A

−

→

A

+ e−) of the acceptor.

Although the anchor points in the diagram are obtained quantitatively, we

generally interpolate the intervening points on the VB configuration curves qual-
itatively, on the basis of spin pairing schemes and VB arguments (123). At this
point it should be stressed that the overall energy profile for the reaction is of
course quantitatively accessible from our quantum-chemical calculations. The ob-
jective of the curve-crossing model analysis is not to generate the overall reaction
profile but to understand how it arises—and a qualitative approach to generating
the VB configuration curves is generally adequate for this purpose. If we consider
first the product configuration, its energy is lowered during the course of the re-
action because of bond formation between the radical and attacked carbon. At the
same time, the relative energy of the reactant configuration increases because the
π bond on the attacked alkene is stretched, and this is not compensated for by
bond formation with the attacking radical. The energies of the charge-transfer

Vol. 9

COMPUTATIONAL CHEMISTRY

365

configurations are initially very high in energy, but are stabilized by Coloumb
attraction as the reactants approach one another.

The overall energy profile for the reaction can be formed from the resonance

interaction of these contributing configurations. In the early stages of the reac-
tion, the reactant configuration is significantly lower in energy than the others and
dominates the ground-state wave function. However, in the vicinity of the tran-
sition structure, the reactant and product configurations have similar energies,
and thus significant mixing is possible. This stabilizes the wave function, with the
strength of the stabilizing interaction increasing with the decreasing energy dif-
ference between the alternate configurations. It is this mixing of the reactant and
product configurations which leads to the avoided crossing, and accounts for bar-
rier formation. Beyond the transition structure, the product configuration is lower
in energy and dominates the wave function. The charge-transfer configurations
generally lie significantly above the ground-state wave function for most of the
reaction. However, in the vicinity of the transition structure, they can sometimes
be sufficiently low in energy to interact. In those cases, the transition structure is
further stabilized, and (if one of the charge transfer configurations is lower than
the other) the mixing is reflected in a degree of partial charge transfer between
the reactants. Since the charge distribution within the transition structure is ac-
cessible from quantum-chemical calculations, this provides a testable prediction
for the model.

Using this state-correlation diagram, in conjunction with simple VB argu-

ments, the curve-crossing model can be used to predict the influence of various
energy parameters on the reaction barrier. For radical addition to alkenes (103),
the barrier depends mainly on the reaction exothermicity (which measures the
energy difference between the reactant and product configurations at their opti-
mal geometries), the singlet-triplet gap in the alkene (which measures the energy
difference between the reactant and product configurations at the reactant ge-
ometry), and the relative energies of the possible charge-transfer configurations.
The effects of individual variations in these quantities are illustrated graphically
in Figure 12. It can be seen that the barrier height is lowered by an increase in
the reaction exothermicity, a decrease in the singlet-triplet gap, or a decrease in
the relative energy of one or both of the charge-transfer configurations (provided
that these are sufficiently low in energy to contribute to the ground-state wave
function).

A strategy for understanding the effects substituents in the barriers of rad-

ical reactions, such as addition, is to calculate these key quantities [ie, the re-
action exothermicity, the singlet-triplet excitation gap of the closed-shell sub-
strate(s), and the energy for charge transfer between the reactants], and look
for relationships between these quantities and the barrier heights. In this way,
one could establish, for example, the extent of polar interactions in a particular
class of reactions. As noted at the beginning of this section, a number of such
studies have already been performed for the key reactions in free-radical poly-
merization. The curve-crossing analysis of radical addition reactions, which is
reviewed in detail elsewhere (94,99), indicate that, in the absence of polar in-
teractions, the barrier height depends on the reaction exothermicity, in accor-
dance with the Evans–Polanyi rule (129). However, for combinations of electron-
withdrawing and electron-donating reactants, polar interactions are significant,

366

COMPUTATIONAL CHEMISTRY

Vol. 9

Fig. 12.

State correlation diagrams showing separately the qualitative effects of (a) in-

creasing the reaction exothermicity, (b) decreasing the singlet–triplet gap, and (c) decreas-
ing the energy of the charge-transfer configuration. For the sake of clarity the adiabatic
minimum energy path showing the avoided crossing, as in Figure 9, is omitted from (a)
and (b).

and cause substantial deviation from Evans–Polanyi behavior. More recently, the
curve-crossing model has been used to explain the relative reactivity of the C C,
C O, and C S bonds (which is of relevance to RAFT polymerization) (130), and
to examine why alkynes are less reactive to addition than alkenes (131). In these
cases the differing singlet-triplet gaps of the alternative substrates are also im-
portant in governing their relative reactivities. Curve-crossing studies have also
been applied to various types of hydrogen abstraction reactions (115,116,124),
and, depending upon the substituents, the singlet-triplet gaps (in this case of both
the reactant and product substrates), exothermicities, and polar interactions have
all been found to be important in governing reactivity in these reactions.

The Future

There are many more potential applications for quantum chemistry in free-radical
polymerization. As computer power increases, one will be able to calculate rate
coefficients for yet larger systems, and at much greater accuracy. In this way,
quantum chemistry has a role to play in modeling the kinetics and mechanism of
polymerization processes, and predicting the properties of the resulting polymer
(such as the molecular weight distribution, the copolymer composition, and the
degree of branching). Quantum chemistry may also help in designing improved
agents or processes for controlling polymerization, and in identifying their mech-
anisms. Of course, quantum chemistry may also be used effectively to study other
types of polymerization processes, such as ring-opening processes. Provided care

Vol. 9

COMPUTATIONAL CHEMISTRY

367

is taken to ensure the accuracy of the methods used, quantum-chemical methods
provide a powerful tool for studying free-radical polymerization, and should be
seen as a valuable complement to experimental approaches.

BIBLIOGRAPHY

1. M. L. Coote, T. P. Davis, and L. Radom, Theochem 461/462, 91–96 (1999);

Macro-

molecules 32, 2935–2940 (1999); Macromolecules 32, 5270–5276 (1999).

2. J. P. A. Heuts, R. G. Gilbert, and I. A. Maxwell, Macromolecules 30, 726 (1997).
3. J. P. A. Heuts, R. G. Gilbert, and L. Radom, Macromolecules 28, 8771–8781

(1995).

4. D. M. Huang, M. J. Monteiro, and R. G. Gilbert, Macromolecules 31, 5175 (1998).
5. J. S.-S. Toh, D. M. Huang, P. A. Lovell, and R. G. Gilbert, Polymer 42, 1915–1920

(2001).

6. M. L. Coote and L. Radom, J. Am. Chem. Soc. 125, 1490 (2003).
7. M. L. Coote and L. Radom, Macromolecules 37, 590 (2004).
8. W. J. Hehre, L. Radom, P. v. R. Schleyer, and J. A. Pople, Ab Initio Molecular Orbital

Theory, John Wiley & Sons, Inc., New York, 1986.

9. F. Jensen, Introduction to Computational Chemistry, John Wiley & Sons, Inc., New

York, 1999.

10. R. G. Parr and W. Yang, Density Functional Theory of Atoms and Molecules, Oxford

University Press, New York, 1989.

11. W. Koch and M. C. Holthausen, A Chemist’s Guide to Density Functional Theory,

Wiley-VCH, Weinheim, 2000.

12. S. W. Benson, Thermochemical Kinetics, John Wiley & Sons, Inc., New York, 1976.
13. D. A. McQuarrie, Statistical Mechanics, Harper & Row, New York, 1976.
14. R. G. Gilbert and S. C Smith, Theory of Unimolecular and Recombination Reactions,

Blackwell Scientific, Oxford, 1990.

15. J. I. Steinfeld, J. S. Francisco, and W. L. Hase, Chemical Kinetics and Dynamics, 2nd

ed., Prentice-Hall, Englewood Cliffs, N.J., 1999.

16. P. W. Atkins, Physical Chemistry, 6th ed., W. H. Freeman and Co., San Francisco,

2000.

17. E. Schr¨odinger, Ann. Physik 79, 361 (1926).
18. P. Pyykk¨o, Chem. Rev. 88, 563–594 (1988).
19. M. Born and J. R. Oppenheimer, Ann. Physik 84, 457 (1927).
20. A. E. Reed, L. A. Curtiss, and F. Weinhold, Chem Rev. 88, 899 (1988).
21. R. F. W. Bader, Atoms in Molecules: A Quantum Theory, Clarendon Press, Oxford,

1990.

22. C. J. Cramer, Essentials of Computational Chemistry: Theories and Models, John

Wiley and Sons, Inc., New York, 2002.

23. E. G. Lewars, Computational Chemistry: Introduction to the Theory and Applications

of Molecular and Quantum Mechanics, Kluwer Academic Publishers, Dordrecht, the
Netherlands, 2003.

24. A. K. Rappe and C. J. Casewit, Molecular Mechanics Across Chemistry, University

Science Books, Herndon, Va., 1997.

25. J. C. Slater, Phys. Rev. 34, 1293 (1929); J. C. Slater, Phys. Rev. 35, 509 (1930).
26. D. R. Hartree, Proc. Cambridge Philos. Soc. 26, 376 (1928); V. Fock, Z. Phys. 61, 126

(1930).

27. A. K. Wilson, T. van Mourik, and T. H. Dunning Jr., J. Mol. Struct. (Theochem) 388,

339–349 (1996), and references cited therein.

368

COMPUTATIONAL CHEMISTRY

Vol. 9

28. J. M. L. Martin and S. Parthiban, in J. Cioslowski, ed., Quantum Mechanical Predic-

tion of Thermochemical Data, Kluwer Academic Publishers, Dordrecht, the Nether-
lands, 2001, pp. 31–65, and references cited therein.

29. R. Car and M. Parrinello, Phys. Rev. Lett. 55, 2471 (1985).
30. P. E. Bl¨ochl, P. Margl, and K. Schwarz, in B. B. Laird, R. B. Ross, and T. Ziegler,

eds., Chemical Applications of Density Functional Theory, American Chemical Society,
Washington, D.C., 1996, pp. 54–69.

31. D. M. Chipman, Theor. Chim. Acta 82, 93–115 (1992).
32. M. W. Wong and L. Radom, J. Phys. Chem. 99, 8582–8588 (1995).
33. M. W. Wong and L. Radom, J. Phys. Chem. 102, 2237–2245 (1998).
34. Y. Y. Chuang, E. L. Coitino, and D. G. Truhlar, J. Phys. Chem. A 104, 446 (2000).
35. C. C. J. Roothaan, Rev. Mod. Phys. 23, 69 (1951); G. G. Hall, Proc. R. Soc. (London)

A 205, 541 (1951).

36. M. L. Coote, G. P. F. Wood, and L. Radom, J. Phys. Chem. A 106, 12124–12138 (2002).
37. R. G´omez-Balderas, M. L. Coote, D. J. Henry, and L. Radom, J. Phys. Chem. A 108,

2874–2883 (2004).

38. D. K. Malick, G. A. Petersson, and J. A. Montgomery, J. Chem. Phys. 108, 5704–5713

(1998).

39. V. D. Knyazev and I. R. Slagle, J. Phys. Chem. 100, 16899–16911 (1996);

V. D.

Knyazev, A. Bencsura, S. I Stoliarov, and I. R. Slagle, J. Phys. Chem. 100, 11346–
11354 (1996).

40. M. Schwartz, P. Marshall, R. J. Berry, C. J. Ehlers, and G. A. Petersson, J. Phys. Chem.

A 102, 10074–10081 (1998).

41. M. L. Coote, M. A. Collins, and L. Radom, Mol. Phys. 101, 1329–1338 (2003).
42. J. A. Montgomery Jr., M. J. Frisch, J. W. Ochterski, and G. A. Petersson, J. Chem.

Phys. 112, 6532–6542 (2000), and references cited therein.

43. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman,

J. A. Montgomery Jr., T. Vreven, K. N. Kudin, J. C. Burant, J. M. Millam, S. S. Iyengar,
J. Tomasi, V. Barone, B. Mennucci, M. Cossi, G. Scalmani, N. Rega, G. A. Petersson,
H. Nakatsuji, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T.
Nakajima, Y. Honda, O. Kitao, H. Nakai, M. Klene, X. Li, J. E. Knox, H. P. Hratchian,
J. B. Cross, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O. Yazyev, A. J.
Austin, R. Cammi, C. Pomelli, J. W. Ochterski, P. Y. Ayala, K. Morokuma, G. A. Voth,
P. Salvador, J. J. Dannenberg, V. G. Zakrzewski, S. Dapprich, A. D. Daniels, M. C.
Strain, O. Farkas, D. K. Malick, A. D. Rabuck, K. Raghavachari, J. B. Foresman, J.
V. Ortiz, Q. Cui, A. G. Baboul, S. Clifford, J. Cioslowski, B. B. Stefanov, G. Liu, A.
Liashenko, P. Piskorz, I. Komaromi, R. L. Martin, D. J. Fox, T. Keith, M. A. Al-Laham,
C. Y. Peng, A. Nanayakkara, M. Challacombe, P. M. W. Gill, B. Johnson, W. Chen, M.
W. Wong, C. Gonzalez, and J. A. Pople, Gaussian 03, Revision B.03, Gaussian, Inc.,
Pittsburgh Pa., 2003. Available at http://www.gaussian.com/

44. L. A. Curtiss and K. Raghavachari, Theor. Chem. Acc. 108, 61–70 (2002), and refer-

ences cited therein.

45. D. J. Henry, M. B. Sullivan, and L. Radom, J. Chem. Phys. 118, 4849–4860 (2003).
46. M. W. Schmidt and M. S. Gordon, Annu. Rev. Phys. Chem. 49, 233–266 (1998).
47. J. J. P Stewart, J. Compu. Aided Mol. Des. 4, 1–105 (1990).
48. P. Hohenberg and W. Kohn, Phys. Rev. B 136, 864 (1964).
49. W. Kohn and L. J. Sham, Phys. Rev. A 140, 1133 (1965).
50. J. P. Perdew, A. Savin, and K. Burke, Phys. Rev. A 51, 4531–4540 (1995).
51. J. D. Goddard and G. Orlova, J. Chem. Phys. 111, 7705–7712 (1999).
52. J. C. Slater, Quantum Theory of Molecular Solids, Vol. 4: The Self-Consistent Field for

Molecular Solids, McGraw-Hill, Inc., New York, 1974.

53. S. H. Vosko, L. Wilk, and M. Nusair, Can. J. Phys. 58, 1200 (1980).

Vol. 9

COMPUTATIONAL CHEMISTRY

369

54. W. Kohn, A. D. Becke, and R. G. Parr, J. Phys. Chem. 100, 12974–12980 (1996).
55. A. D. Becke, Phys. Rev. A 38, 3098 (1988).
56. C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988); B. Miehlich, A. Savin,

H. Stoll, and H. Preuss, Chem. Phys. Lett. 157, 200 (1989).

57. J. P. Perdew, K. Burke, and Y. Wang, Phys. Rev. B 54, 16533 (1996), and references

cited therein.

58. A. D. Becke, J. Chem. Phys. 98, 5648 (1993).
59. B. J. Lynch, P. L. Fast, M. Harris, and D. G. Truhlar, J. Phys. Chem. A 104, 4811–4815

(2000);

B. J. Lynch and D. G. Truhlar, J. Phys. Chem. A 105, 2936–2941 (2001);

B. J. Lynch and D. G. Truhlar, J. Phys. Chem. A 106, 842–846 (2002).

60. M. A. Collins, Theor. Chem. Acc. 108, 313–324 (2002).
61. D. H. Zhang, M. A. Collins, and S.-Y. Lee, Science 290, 961–963 (2000).
62. K. Song and M. A. Collins, Chem. Phys. Lett. 335, 481 (2001).
63. A. Sanchez-Galvez, P. Hunt, M. A. Robb, M. Olivucci, T. Vreven, and H. B. Schlegel,

J. Am. Chem. Soc. 122, 2911–2924 (2000); M. Sharma, Y. Wu, and R. Car, J. Chem.
Phys. 95, 821–829 (2003).

64. A. C. T. van Duin, S. Dasgupta, F. Lorant, and W. A. Goddard III, J. Phys. Chem. A

105, 9396–9409 (2001).

65. M. Garavelli, F. Ruggeri, F. Ogliaro, M. J. Bearpark, F. Bernardi, M. Olivucci, and

M. A. Robb, J. Comput. Chem. 24, 1357–1363 (2003).

66. H. Eyring, J. Chem. Phys. 3, 107 (1935).
67. For a more detailed definition of this term, see for example: A. J. Karas, R. G. Gilbert,

and M. A. Collins, Chem. Phys. Lett. 193, 181–184 (1992).

68. D. G. Truhlar, B. C. Garrett, and S. J. Klippenstein, J. Phys. Chem. 100, 12771 (1996).
69. R. T. Skodje, D. G. Truhlar, and B. C. Garrett, J. Phys. Chem. 85, 3019 (1981).
70. B. C. Garrett, D. G. Truhlar, A. F. Wagner, and T. H. Dunning Jr., J. Chem. Phys. 78,

4400 (1983).

71. Y. P. Liu, D. H. Lu, A. Gonzalez-Lafont, D. G. Truhlar, and B. C. Garrett, J. Am. Chem.

Soc. 115, 7806 (1993).

72. J. C. Corchado, Y.-Y. Chuang, P. L. Fast, J. Vill `a, W.-P. Hu, Y.-P. Liu, G. C. Lynch,

K. A. Nguyen, C. F. Jackels, V. S. Melissas, B. J. Lynch, I. Rossi, E. L. Coiti ˜

no,

A. Fernandez-Ramos, J. Pu, T. V. Albu, R. Steckler, B. C Garrett, A. D. Isaacson, and
D. G. Truhlar, POLYRATE 9.1, University of Minnesota, Minneapolis, Minn., 2002.
Available at http://comp.chem.umn.edu/polyrate/

73. A. Kuppermann and D. G. Truhlar, J. Am. Chem. Soc. 93, 1840 (1971); B. C. Garrett,

D. G. Truhlar, R. S. Grev, and A. W. Magnuson, J. Phys. Chem. 84, 1730 (1980).

74. R. P. Bell, The Tunnel Effect in Chemistry, Chapman and Hall, London, 1980.
75. E. Wigner, J. Chem. Phys. 5, 720 (1937).
76. C. Eckart, Phys. Rev. 35, 1303 (1930).
77. Y.-X. Yu, S.-M. Li, Z.-F. Xu, Z.-S. Li, and C.-C. Sun, Chem. Phys. Lett. 302, 281 (1999).
78. W. T. Duncan, R. L. Bell, and T. N. Truong, J. Comput. Chem. 19, 1039 (1998).
79. V. Van Speybroeck, D. Van Neck, M. Waroquier, S. Wauters, M. Saeys, and G. B.

Martin, J. Phys. Chem. A 104 10939–10950 (2000).

80. K. S. Pitzer and W. D. Gwinn, J. Chem. Phys. 10, 428–440 (1942);

J. C. M. Li and

K. S. Pitzer, J. Phys. Chem. 60, 466–474 (1956).

81. A. L. L. East and L. Radom, J. Chem. Phys. 106, 6655 (1997).
82. V. Van Speybroeck, D. Van Neck, and M. Waroquier, J. Phys. Chem. A 106, 8945–8950

(2002).

83. J. F. Stanton, J. Gauss, J. D. Watts, M. Nooijen, N. Oliphant, S. A. Perera, P. G. Szalay,

W. J. Lauderdale, S. A. Kucharski, S. R. Gwaltney, S. Beck, A. Balkov `a, D. E. Bern-
holdt, K. K. Baeck, P. Rozyczko, H. Sekino, C. Hober, and R. J. Bartlett. Integral pack-
ages included are VMOL (J. Alml¨of and P. R. Taylor), VPROPS (P. Taylor), ABACUS

370

COMPUTATIONAL CHEMISTRY

Vol. 9

(T. Helgaker, H. J. Aa. Jensen, P. Jørgensen, J. Olsen, and P. R. Taylor); Quantum
Theory Project, University of Florida, Available at http://www.qtp.ufl.edu/Aces2/

84. M. W. Schmidt, K. K. Baldridge, J. A. Boatz, S. T. Elbert, M. S. Gordon, J. H.

Jensen, S. Koseki, N. Matsunaga, K. A. Nguyen, S. J. Su, T. L. Windus, M. Dupuis,
and J. A. Montgomery, J. Comput. Chem. 14, 1347–1363 (1993).

Available at

http://www.msg.ameslab.gov/GAMESS/GAMESS.html

85. R. D. Amos, A. Bernhardsson, A. Berning, P. Celani, D. L. Cooper, M. J. O. Deegan,

A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, P. J. Knowles, T. Korona, R. Lindh, A. W.
Lloyd, S. J. McNicholas, F. R. Manby, W. Meyer, M. E. Mura, A. Nickla

β, P. Palmieri,

R. Pitzer, G. Rauhut, M. Sch ¨

utz, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni,

T. Thorsteinsson, and H.-J. Werner, Available at http://www.molpro.net/

86. J. Kong, C. A. White, A. I. Krylov, C. D. Sherrill, R. D. Adamson, T. R. Furlani,

M. S. Lee, A. M. Lee, S. R. Gwaltney, T. R. Adams, C. Ochsenfeld, A. T. B. Gilbert, G. S.
Kedziora, V. A. Rassolov, D. R. Maurice, N. Nair, Y. Shao, N. A. Besley, P. E. Maslen,
J. P. Dombroski, H. Daschel, W. Zhang, P. P. Korambath, J. Baker, E. F. C. Byrd,
T. Van Voorhis, M. Oumi, S. Hirata, C.-P. Hsu, N. Ishikawa, J. Florian, A. Warshel,
B. G. Johnson, P. M. W. Gill, M. Head-Gordon, and J. A. Pople, J. Comput. Chem. 21,
1532–1548 (2000). Available at http://www.q-chem.com/

87. Wavefunction, Inc., Available at http://www.wavefun.com/
88. G. Schaftenaar, Available at http://www.cmbi.kun.nl/

∼schaft/molden/molden.html

89. Cambridgesoft, Available at http://www.cambridgesoft.com/
90. N. J. R. van Eikema Hommes, Available at http://www.ccc.uni-erlangen.de/molecule/
91. Free Open Source Software, Available at http://jmol.sourceforge.net/
92. B. M. Bode and M. S. Gordon, J. Mol. Graphics Mod. 16, 133–138 (1998). Available

at http://www.msg.ameslab.gov/GAMESS/GAMESS.html

93. J. P. A. Heuts, R. G. Gilbert, and L. Radom, J. Phys. Chem. 100, 18997–19006 (1996).
94. M. Saeys, M.-F. Reyniers, G. B. Marin, V. van Speybroeck, and M. Waroquier, J. Phys.

Chem. A 107, 9147–9159 (2003).

95. A. P. Scott and L. Radom, J. Phys. Chem. 100, 16502–16513 (1996).
96. J. Chiefari, Y. K. B. Chong, F. Ercole, J. Krstina, J. Jeffery, T. P. T. Le, R. T. A.

Mayadunne, G. F. Meijs, C. L. Moad, G. Moad, E. Rizzardo, and S. H. Thang, Macro-
molecules 31, 5559 (1998).

97. J. K. Kang and C. B. Musgrave, J. Chem. Phys. 115, 11040 (2001).
98. H. Basch and S. Hoz, J. Phys. Chem. A 101, 4416 (1997).
99. M. L. Coote, J. Phys. Chem. A 108, 3865–3872 (2004).

100. M. L. Coote, T. P. Davis, B. Klumperman, and M. J. Monteiro, J. Macromol. Sci., Rev.

Macromol. Chem. Phys. C 38, 567 (1998), and references cited therein.

101. M. Ben-Nun and R. D. Levine, Int. Rev. Phys. Chem. 14, 215–270 (1995); G. A. Voth

and R. M. Hochstrasser, J. Phys. Chem. 100, 13034–13049 (1996).

102. J. R. Pliego Jr. and J. M. Riveros, J. Phys. Chem. A 105, 7241–7247 (2001).
103. H. Fischer and L. Radom, Angew. Chem., Int. Ed. Engl. 40, 1340–1371 (2001).
104. M. L. Coote and T. P. Davis, Prog. Polym. Sci. 24, 1217 (1999), and references cited

therein.

105. M. Deady, A. W. H. Mau, G. Moad, and T. H. Spurling, Makromol. Chem. 194, 1691

(1993).

106. J. Filley, T. McKinnon, D. T. Wu, and G. H. Ko, Macromolecules 35, 3731–3738 (2002).
107. M. Buback, C. H. Kurz, and C. Schmaltz, Macromol. Chem. Phys. 199, 1721–1727

(1998).

108. J. P. A. Heuts, in K. Matyjaszewski and T. P. Davis, eds., Handbook of Radical Poly-

merization, John Wiley & Sons, Inc., New York, 2002, pp. 1–76.

109. C. Barner-Kowollik, J. F. Quinn, D. R. Morsley, and T. P. Davis, J. Polym. Sci., Part A:

Polym. Chem. 39, 1353 (2001).

Vol. 9

CONTROLLED RELEASE FORMULATION, AGRICULTURAL

371

110. Y. Kwak, A. Goto, Y. Tsuji, Y. Murata, K. Komatsu, and T. Fukuda, Macromolecules

35, 3026 (2002).

111. C. Barner-Kowollik, M. L. Coote, T. P. Davis, L. Radom, and P. Vana, J. Polym. Sci.,

Part A: Polym. Chem. 41, 2828–2832 (2003), and references cited therein.

112. M. L. Coote, M. D. Zammit, and T. P. Davis, Trends Polym. Sci. (Cambridge, U.K.) 4,

189 (1996).

113. J. P. A. Heuts, A. Pross, and L. Radom, J. Phys. Chem. 100, 17087–17089 (1996).
114. D. L. Reid, D. A. Armstrong, A. Rauk, and C. von Sonntag, Phys. Chem. Chem. Phys.

5, 3994–3999 (2003).

115. L. Song, W. Wu, K. Dong, P. C. Hibberty, and S. Shaik, J. Phys. Chem. A 106, 11361–

11370 (2002).

116. S. Shaik, W. Wu, K. Dong, L. Song, and P. C. Hiberty, J. Phys. Chem. A 105, 8226

(2001).

117. S. L. Boyd and R. J. Boyd, J. Phys. Chem. A 105, 7096 (2001).
118. R. Sumathi, H. H. Carstensen, and W. H. Green Jr., J. Phys. Chem. A 105, 8969–8984

(2001).

119. R. Sumathi, H. H. Carstensen, and W. H. Green Jr., J. Phys. Chem. A 105, 6910–6925

(2001).

120. C. J. Hawker, A. W. Bosman, and E. Harth, Chem. Rev. 101, 3661–88 (2001).
121. J. S. Wang and K. Matyjaszewski, J. Am. Chem. Soc. 117, 5614–5615 (1995).
122. S. C. Farmer and T. E. Patten, J. Polym. Sci., Part A: Polym. Chem. 40, 555 (2002).
123. A. Pross and S. S. Shaik, Acc. Chem. Res. 16, 363 (1983);

A. Pross, Adv. Phys. Org.

Chem. 21, 99–196 (1985); S. S. Shaik, Prog. Phys. Org. Chem. 15, 197–337 (1985).

124. A. Pross, H. Yamataka, and S. Nagase, J. Phys. Org. Chem. 4, 135 (1991).
125. A. Pross, Theoretical and Physical Principles of Organic Reactivity, John Wiley &

Sons, Inc., New York, 1995.

126. W. Heitler and F. London, Z. Phys. 44, 455 (1927).
127. L. Pauling, The Nature of the Chemical Bond, Cornell University Press, Ithaca, N.Y.,

1939.

128. R. B. Woodward and R. Hoffman, Angew Chem. 81, 797 (1969); R. B. Woodward and

R. Hoffman, Angew Chem. Int. Ed. Engl. 8, 781 (1969).

129. M. G. Evans, Disc. Faraday Soc. 2, 271–279 (1947); M. G. Evans, J. Gergely, and E. C.

Seaman,
J. Polym. Sci. 3, 866–879 (1948).

130. D. J. Henry, M. L. Coote, R. G´omez-Balderas, and L. Radom, J. Am. Chem. Soc. 126,

1732–1740 (2004).

131. R. G´omez-Balderas, M. L. Coote, D. J. Henry, H. Fischer, and L. Radom, J. Phys.

Chem. A 107, 6082–6090 (2003).

M

ICHELLE

L. C

OOTE

Australian National University

CONDUCTIVE POLYMER COMPOSITES.

See Volume 5.

CONFORMATION AND CONFIGURATION.

See Volume 2.

372

CONTROLLED RELEASE FORMULATION, AGRICULTURAL

Vol. 9

CONTACT LENSES.

See H

YDROGELS

.

CONTROLLED MICROSTRUCTURE.

See S

INGLE

S

ITE

C

ATALYSTS

.

CONTROLLED RADICAL POLYMERIZATION.

See L

IVING

R

ADICAL

P

OLYMERIZATION

.