The Schrödinger equation for N electrons and M nuclei of a

molecule

H(r,R) Ψ(r,R,t) = i h 䌍Ψ(r,R,t)/䌍t or

H(r,R) Ψ(r,R) = E Ψ(r,R)

|Ψ(r,R)|2 gives probability density for finding electrons at

r = r1r2 r3 ... rN and nuclei at R1 R2 R3 ...RM .

H contains electronic kinetic energy Te = -h2/2 Σj=1,N me-1 ∇j2

nuclear kinetic energy TM= -h2/2 Σj=1,M mj-1 ∇j2
electron-nuclei Coulomb potentials - Σj=1,MZj Σk=1,N e2/|rk-Rj|
VeM nuclear-nuclear Coulomb repulsions Σj<k=1,M ZjZke2/|Rk-Rj|
and electron-electron Coulomb repulsions Vee= Σj<k=1,Ne2/rj,k

It can contain more terms if, for example, external electric

or magnetic fields are present (e.g., Σk=1,N erkE).
In the Born-Oppenheimer approximation/separation, we
ignore the TM motions of the nuclei (pretend the nuclei are fixed at
specified locations R) and solve
H0 ψΚ(r|R) =EK(R) ψΚ(r|R)
the so-called electronic Schrödinger equation. H0 contains all of
H except TM.
Because H0 is Hermitian, its eigenfunctions form a complete set
of functions of r. So, the full Ψ can be expanded in the ψK :
Ψ(r,R) = ΣK ψK(r,R) χK(R) .
The ψK(r,R) depend on R because H0 does through
- Σj=1,MZj Σk=1,N e2/|rk-Rj|.
This expansion can then be used in
H(r,R) Ψ(r,R) = E Ψ(r,R)
[H0 -h2/2 Σj=1,M mj-1 ∇j2 -E] ΣK ψK(r,R) χK(R) = 0
to produce equations for the χK(R):
0 = [EL(R) -h2/2 Σj=1,M mj-1 ∇j2 -E] χL(R)
+ ΣK< ψL(r,R)| -h2/2 Σj=1,M mj-1 ∇j2 ψK(r,R)> χK(R)
+ ΣK< ψL(r,R)| -h2Σj=1,M mj-1 ∇j ψK(r,R)> ∇j χK(R)
These are the coupled-channel equations.
If we ignore all of the non-adiabatic terms, we obtain a SE
For the vib./rot./trans. Motion
0 = [EL(R) -h2/2 Σj=1,M mj-1 ∇j2 -E] χL(R)
Each electronic state L has its own set of rot./vib. wave functions
and energies
[EL(R) -h2/2 Σj=1,M mj-1 ∇j2 -EL,J,M,ν] χL,J,M, ν(R) = 0
30000

*****
25000 *******
****
* * ****
*
* *
***
*
*****
20000 ***
* *** !!
!! !!
** !!! !!
E in cm-1

* *** !!!
!
** !!!!
* *** !!
! !!!
15000 ** !
!!
! * *** !!
!!!
! * *** !!!
!
* !!
! * **
** !!!!!
* ** !
!
10000 ! * ** !!!!
! * *** !
!!!
* !
! **
! **
! ** *****
***
!!
!
!
!!
!
!! !i
5000 *** !!
!
!
!!
!!
***
*************
!! !!!
!
!!!
!f
!! ! !
!! !!!!
0 !!!!!!! !!
!! !!!

R in Angstroms

The non-adiabatic couplings can induce transitions among these

states (radiationless transitions).
There are major difficulties in solving the electronic SE:

Vee makes the equation not separable- this means ψ is not

rigorously a product of functions of individual electron coordinates.

ψ ≠ φ1 φ2 φ3 φ4 φ5 (e.g., 1sα(1) 1sβ(2) 2sα(3) 2sβ(4) 2p1α(5))

Cusps

The factors (1/rk ∂/∂rk –Ze2/rk) ψ and (1/rk,l ∂/∂rk,l +e2/rk,l) ψ

will blow up unless so-called cusp conditions are obeyed by ψ:

∂/∂rk ψ = Ze2 ψ(as rk→0) and ∂/∂rk,l ψ = - e2 ψ(as rk,l→0).

 Cusp near nucleus Cusp as two electrons
approach

This means when we try to approximately solve the electronic SE,

we should use trial functions that have such cusps. Slater-type
orbitals (exp(-ζr)) have cusps at nuclei, but Gaussians (exp(-αr2))
do not. We rarely use functions with e-e cusps, but we should.
Addressing the non-separability problem:

If Vee could be replaced by a one-electron additive potential

VMF = Σj=1,N VMF(rj)

the solutions ψ would be products (actually antisymmetrized

products called Slater determinants) of functions of individual
electron coordinates (spin-orbitals):

ψ = | φ1(r1) φ2 (r2) φ3 (r3) φ4 (r4) φ5 (r5) |

= (N!)-1/2 ΣP=1,N! P φ1(r1) φ2 (r2) φ3 (r3) φ4 (r4) φ5 (r5)

Before considering finding a VMF, let’s examine how important

antisymmetry is by considering two electrons in π and π* orbitals.
Singlet π2 |πα(1) πβ(2)|= 2-1/2 π(1)π(2)[α1β2−β1α2]
Triplet ππ* |πα(1) π∗α(2)|= 2-1/2 [π(1)π∗(2) − π(2)π∗(1)] α1α2
|πβ(1) π∗β(2)|= 2-1/2 [π(1)π∗(2) − π(2)π∗(1)] β1β2
2-1/2 [|πα(1) π∗β(2)|+|πβ(1) π∗α(2)] =
2−1[π(1)π∗(2)α1β2 + π(1)π∗(2)β1α2
−π∗(1)π(2)β1α2 − π∗(1)π(2)α1β2]
Singlet ππ∗ 2-1/2 [|πα(1) π∗β(2)| - |πβ(1) π∗α(2)] =
2−1[π(1)π∗(2)α1β2 − π(1)π∗(2)β1α2
−π∗(1)π(2)β1α2 + π∗(1)π(2)α1β2]
Singlet π∗2 |π∗α(1) π∗β(2)|= 2-1/2 π∗(1)π∗(2)[α1β2−β1α2]
Now think of π = 2-1/2 (L + R) and π∗ = 2-1/2 (L-R)
|πα(1) πβ(2)| = 2-1[|Rα(1) Rβ(2)|+ |Lα(1) Lβ(2)|
+ |Rα(1) Lβ(2)|+ |Lα(1) Rβ(2)|] ionic + diradical
|π∗α(1) π∗β(2)| = 2-1[|Rα(1) Rβ(2)|+ |Lα(1) Lβ(2)|
-|Rα(1) Lβ(2)|- |Lα(1) Rβ(2)|] ionic + diradical

2-1/2 [|πα(1) π∗β(2)| - |πβ(1) π∗α(2)] = 2−3/2[ |Rα(1) Rβ(2)|+

|Lα(1) Rβ(2)|- |Rα(1) Lβ(2)|- |Lα(1) Lβ(2)|]
-2-3/2[|Rβ(1) Rα(2)|+|Lβ(1) Rα(2)|- |Rβ(1) Lα(2)|- |Lβ(1) Lα(2)|]
= 2−1/2[ |Rα(1) Rβ(2)|+|Lα(1) Rβ(2)|] ionic

|πα(1) π∗α(2)| = 2-1[|Lα(1) Rα(2)|- |Rα(1) Lα(2)|]

= |Lα(1) Rα(2)| diradical
 !*2 Ionic

Ionic +
Diradical
Singlet !!*
Triplet !!*
Diradical

!2
0 Twist Angle 90

To adequately describe the π2 bond breaking, we need to mix the

π2 and π∗2 configuration state functions (CSF). This shows how
single configuration functions may not be adequate.
|πα(1) πβ(2)| = 2-1[|Rα(1) Rβ(2)|+ |Lα(1) Lβ(2)|
+ |Rα(1) Lβ(2)|+ |Lα(1) Rβ(2)|] ionic + diradical
|π∗α(1) π∗β(2)| = 2-1[|Rα(1) Rβ(2)|+ |Lα(1) Lβ(2)|
-|Rα(1) Lβ(2)|- |Lα(1) Rβ(2)|] ionic + diradical
So, one must combine 2-1/2{|πα(1) πβ(2)| - |π∗α(1) π∗β(2)|}
to obtain a diradical state and
2-1/2{|πα(1) πβ(2)| + |π∗α(1) π∗β(2)|}
to obtain an ionic state.
Analogous “trouble” occurs whenever one uses a single
determinant HF wave function to describe a bond that one wants to
break:

H2 (σ2) → H(1sA) + H(1sB)

H3C-CH3 (σ2) → H3C• + •CH3

As we will see soon, one can partially solve this problem by using
a so-called unrestricted HF (UHF) wave function, but it has
problems as well.
How does one find a VMF? One way is to postulate that

ψ= |φ1 φ2 φ3 ...φN| ,

and write down <ψ H ψ> using the Slater-Condon rules:

<ψ H ψ> = Σk< φk|Te + Ve,n + Vn,n| φk> + 1/2 Σk,l

[< φk(1) φl(2)|e2/r1,2| φk(1) φl(2)> - < φk(1) φl(2)|e2/r1,2| φl(1) φk(2)>]

and observe that Coulomb (J) and exchange (K) interactions among
spin-orbitals arise. If one also minimizes this energy with respect to
the φ’s, one obtains equations h φJ = εJ φJ = [Te + Ve,n + Vn,n] φJ
+ Σk [< φk(1)|e2/r1,2| φk(1)> φJ(2) - < φk(1)|e2/r1,2| φJ(1)> φk(2).

that contain the J and K potentials.

J1,2= ∫ |φ1(r)|2 e2/|r-r’|φ2(r’)|2 dr dr’

K1,2= ∫ φ1(r) φ2(r’) e2/|r-r’|φ2(r) φ1(r’)dr dr’

!2(r')
Overlap region

!1(r)
So, one is motivated to define VMF in terms of the J and K
interactions. This is the Hartree-Fock definition of VMF. It has the
characteristic that < ψ H ψ > = <ψ H0 ψ >, so if H – H0 is viewed
as a perturbation and H0 is defined as

H0 = Te + Ve,n + Vn,n + (J-K)

there is no first-order perturbation correction to the energy. This

choice of H0 forms the basis of Møller-Plesset perturbation theory
(MPn).

It is by making a mean-field model that our (chemists’) concepts

of orbitals and of electronic configurations (e.g., 1s α1s β 2s α 2s β
2p1 α) arise.
Another good thing about HF orbitals is that their energies εK give
approximate ionization potentials and electron affinities
(Koopmans’ theorem). This can be shown by writing down the
energies of two Slater determinants

ψ0 = |φ1 φ2 φ3 ...φN| and ψ- = |φ1 φ2 φ3 ...φN φN+1| ,

using the energy expression Σk< φk|Te + Ve,n + Vn,n| φk> + 1/2 Σk,l

[< φk(1) φl(2)|e2/r1,2| φk(1) φl(2)> - < φk(1) φl(2)|e2/r1,2| φl(1) φk(2)>]

and subtracting the two energy expressions to obtain the energy

difference. You try it as a homework problem and see if you can
show the energy difference is indeed εN+1.
The sum of the orbital energies is not equal to the HF energy:

Ε = Σk< φk|Te + Ve,n + Vn,n| φk> + 1/2 Σk,l

[< φk(1) φl(2)|e2/r1,2| φk(1) φl(2)> - < φk(1) φl(2)|e2/r1,2| φl(1) φk(2)>]

εk = < φk|Te + Ve,n + Vn,n| φk> + Σl

[< φk(1) φl(2)|e2/r1,2| φk(1) φl(2)> - < φk(1) φl(2)|e2/r1,2| φl(1) φk(2)>]

The Brillouin theorem holds:

< |φ1 φ2 φa ...φN| H |φ1 φ2 φm ...φN| > = <φa| Te + Ve,n + Vn,n|φm> + Σl

[< φa(1) φl(2)|e2/r1,2| φm(1) φl(2)> - < φa (1) φl(2)|e2/r1,2| φl(1) φm(2)>]

= <φa|hHF | φm> = 0
The Slater-Condon rules- memorize them
(i) If | > and | ' > are identical, then
<|F+G|>=
Σi < φi | f | φi > +Σi>j [< φiφj | g | φiφj > - < φiφj | g | φjφi >],
where the sums over i and j run over all spin-orbitals in | >;
(ii) If | > and | ' > differ by a single spin-orbital ( φp ≠ φ'p ),
<|F+G|'>=
< φp | f | φ'p > +Σj [< φpφj | g | φ'pφj > - < φpφj | g | φjφ'p >],
where the sum over j runs over all spin-orbitals in | > except φp ;
(iii) If | > and | ' > differ by two ( φp ≠ φ'p and φq ≠ φ'q),
<|F+G|'>=
< φp φq | g | φ'p φ'q > - < φp φq | g | φ'q φ'p >
(note that the F contribution vanishes in this case);
(iv) If | > and | ' > differ by three or more spin orbitals, then
< | F + G | ' > = 0;
(v) For the identity operator I, < | I | ' > = 0 if | > and | ' > differ by one
or more spin-orbitals.
Some single-configuration functions are not single determinants.
There are cases where more than one determinant must be used.
Although the determinant |1sα 1sβ 2sα 2sβ 2pzα 2pyα| is an
acceptable approximation to the carbon 3P state if the 1s and 2s spin-
orbitals are restricted to be equal for α and β spins, the 1S state arising
in this same 1s22s22p2 configuration can not be represented as a single
determinant. The 1S state requires a minimum of the following three-
determinant wave function:

Ψ = 3-1/2 [1sα 1sβ 2sα 2sβ 2pzα 2pzβ|

- 1sα 1sβ 2sα 2sβ 2pxα 2pxβ| - 1sα 1sβ 2sα 2sβ 2pyα 2pyβ| ].

If a state cannot be represented by a single determinant, one should

not use theoretical methods that are predicated on a dominant single
determinant in the expansion of the full wave function.
We have dealt with the non-separability issue, but what about the
cusps? Is doing so necessary? Yes it is!

Example- carbon atom’s total electronic energy is – 1030.080 eV

and J2px,2py = 13 eV, so the J’s (and K’s) are large quantities on a
“chemical” scale of 1 kcal/mol. The Be 1s/1s interaction in the HF
approximation and in reality differ a lot.
So, the electron-electron interactions are large quantities and the
errors made in describing them in terms of the HF mean-field
picture are also large.

Why don’t we use ψ functions that have electron-electron cusps?

Sometimes we do (explicitly correlated wave functions are used in
so-called r-12 methods), but this results in very difficult theories
to implement and very computer-intensive calculations. We’ll
here more later from Martin Head-Gordon about this.
The most common way to improve beyond the HF |φ1 φ2 φ3 ...φN| is
to use trial wave functions of the so-called configuration
interaction (CI) form ψ = ΣL CL1,L2,...LN |φL1 φL2 φL3 ...φLN|. This
makes mathematical sense because the determinants |φL1 φL2 φL3 ...
φLN| form orthonormal complete sets, so ψ can be so expanded.
Physically, what does this mean? Here is a useful identity for two
determinants that one can use to interpret such CI wave functions:

Ψ = C1 | ..φα φβ..| - C2 | ..φ'α φ'β..|

= C1/2 { | ..( φ - xφ')α ( φ + xφ')β..| - | ..( φ - xφ')β ( φ + xφ')α..| }.

with x = (C2/C1)1/2

So a combination of two determinants that differ by doubly

occupied orbital φ being replaced by doubly occupied φ’ is
equivalent to singlet 2-1/2 (αβ - βα) coupled polarized orbital
pairs φ - xφ' and φ + xφ'.
For example π2 → π*2 CI in olefins or 2s2 → 2p2 CI in alkaline
earth atoms produce the following polarized orbital pairs.

!#

! + x!# ! "x!# 2s - a 2p z
left polarized right polarized

2s + a 2p z
2s and 2p z
Placing electrons into different polarized orbital pairs allows
them to avoid one another and thus correlate their motions. This
correlation is how the wave functions attempt to approach the e-e
cusp condition.
Sometimes the CI is essential- for example, to adequately
describe breaking the π bond in the singlet state of an olefin.
However, CI is always important if one wishes to include
electron-electron avoidance that is called dynamical correlation.
In all cases, it is useful to keep in mind the polarized orbital pair
model:

Ψ = C1 | ..φα φβ..| - C2 | ..φ'α φ'β..|

= C1/2 { | ..( φ - xφ')α ( φ + xφ')β..| - | ..( φ - xφ')β ( φ + xφ')α..| }

Let’s get a bit more specific. How does one determine the orbitals φJ
and then how does one determine the CI coefficients CJ?

The orbitals are usually determined by carrying out a HF calculation.

This is not done (except in rare cases) by solving the HF differential
equations on a spatial grid but by expanding the φJ in terms of so-
called atomic orbital (AO) (because they usually are centered on
atoms) basis functions- the LCAO-MO expansion:

φJ = Σµ χµ CJ,µ.

This reduces the HF calculation to a matrix eigenvalue form

Σµ <χν |he| χµ> CJ,µ = εJ Σµ <χν|χµ> CJ,µ

The matrix elements needed to carry out such a calculation are

<χν| he| χµ> = <χν| –h2/2m ∇2 |χµ> + Σa<χν| -Zae2/|ra |χµ>

+ ΣK CK,η CK,γ [<χν(r) χη(r’) |(e2/|r-r’|) | χµ(r) χγ(r’)>

-<χν(r) χη(r’) |(e2/|r-r’|) | χγ(r) χµ(r’)>]and the overlap integrals

<χν|χµ>.

The number of these one- and two electron integrals scales with
the basis set size M as M2 and M4. The computer effort needed to
solve the MxM eigenvalue problem scales as M3. The sum over K
runs over all of the occupied spin-orbitals.
UHF Wavefunctions are not eigenfunctions of S2

<χν| he| χµ> = <χν| –h2/2m ∇2 |χµ> + Σa<χν| -Zae2/|ra |χµ>

+ ΣK CK,η CK,γ [<χν(r) χη(r’) |(e2/|r-r’|) | χµ(r) χγ(r’)>

-<χν(r) χη(r’) |(e2/|r-r’|) | χγ(r) χµ(r’)>].

The matrix elements of the Fock operator are different for an α and
a β spin-orbital because the sum ΣK CK,η CK,γ appearing in these
matrix elements runs over all N of the occupied spin-orbitals. If the
spin-orbital being solved for is of α type, there will be Coulomb
integrals for K = 1sα, 1sβ, 2sα, 2sβ, 2pzα, and 2pyα and exchange
contributions for K = 1sα, 2sα, 2pzα, and 2pyα. On the other hand,
when solving for spin-orbitals of β type, there will be Coulomb
integrals for K = 1sα, 1sβ, 2sα, 2sβ, 2pzα, and 2pyα. but exchange
contributions -only for K =1sβ and 2sβ.
The UHF wave function can be used to describe bond
breaking such as
H2 (σ2) → H(1sA) + H(1sB) and H3C-CH3 (σ2) → H3C• + •CH3
However, the resulting energy curves can have slope jumps.
σ2 σ σ’
Slater-type orbitals (STOs)

χn,l,m (r,θ,φ) = Nn,l,m,ζ Yl,m (θ,φ) rn-1 e-ζr

are characterized by quantum numbers n, l, and m and exponents

(which characterize the radial 'size' ) ζ.

Cartesian Gaussian-type orbitals (GTOs)

χa,b,c (r,θ,φ) = N'a,b,c,α xa yb zc exp(-αr2),

are characterized by quantum numbers a, b, and c, which detail the

angular shape and direction of the orbital, and exponents α which
govern the radial 'size’.

Of course, for both functions, they are also characterized by where

they are located (e.g., a nucleus or bond midpoint).
Slater-type orbitals are similar to Hydrogenic orbitals in the regions
close to the nuclei. Specifically, they have a non-zero slope near the
nucleus on which they are located (i.e., d/dr(exp(-ζr))r=0 = -ζ, so
they can have proper electron-nucleus cusps.
In contrast, GTOs have zero slope near r=0 because
d/dr(exp(-αr2))r=0 = 0.
This characteristic favors STOs over GTOs because we know that
the correct solutions to the Schrödinger equation have such cusps at
each nucleus of a molecule.
However, the multi-center integrals which arise in polyatomic-
molecule calculations cannot efficiently be evaluated when STOs
are employed. In contrast, such integrals can routinely be computed
when GTOs are used. This advantage of GTOs has lead to the
dominance of these functions in molecular quantum chemistry.
To overcome the cusp weakness of GTO functions, it is common to
combine two, three, or more GTOs, with combination coefficients
that are fixed and not treated as LCAO parameters, into new
functions called contracted GTOs or CGTOs. However, it is not
possible to correctly produce a cusp by combining any number of
Gaussian functions because every Gaussian has a zero slope at r = 0
as shown below.
tight Gaussian

orbital with cusp at r = 0

loose Gaussian

medium Gaussian

r
Most AO basis sets contain a mixture of different classes of
functions.
Fundamental core and valence basis functions
Polarization functions
Diffuse functions
Rydberg functions
Minimal basis-the number of CGTOs equals the number of core
and valence atomic orbitals in the atom.
Carbon- one tight s-type CGTO, one looser s-type CGTO and a set
of three looser p-type CGTOs.
Double-zeta (DZ)- twice as many CGTOs as there are core and
valence atomic orbitals.
Carbon- two tight s, two looser s, and two sets of three looser p
CGTOs.
The use of more basis functions is motivated by a desire to provide
additional variational flexibility so the LCAO process can generate
molecular orbitals of variable diffuseness as the local
electronegativity of the atom varies.
Triple-zeta (TZ)- three times as many CGTOs as the number of
core and valence atomic orbitals (extensions to quadruple-zeta and
higher-zeta bases also exist).
Polarization functions- one higher angular momentum than appears
in the atom's valence orbital space.

d-functions for C, N, and O and p-functions for H with exponents (ζ

or α) which cause their radial sizes to be similar to the sizes of the
valence orbitals.
Note- the polarization p orbitals of H are similar in size to the
valence 1s orbital and the polarization d orbitals of C are similar in
size to the 2s and 2p orbitals, not like the valence d orbitals of C.
Polarization functions give angular flexibility to the LCAO process
in forming molecular orbitals between from valence atomic orbitals.
Polarization functions also allow for angular correlations in
describing the correlated motions of electrons.
An example of d polarization functions on C and O

C O C O

Carbon p! and d! orbitals combining to

form a bent ! orbital

C O C O

Oxygen p! and d! orbitals combining to form

a bent ! orbital

C O

! bond formed from C and O bent (polarized) AOs

Valence and polarization functions do not provide enough radial
flexibility to adequately describe very diffuse charge densities.

The diffuse basis functions tabulated on the PNNL web site are
appropriate if the anion under study has its excess electron in a
valence-type orbital (e.g., as in F-, OH-, carboxylates, etc.) but not
for very weakly bound anions (e.g., having EAs of 0.1 eV or less).

For an electron in a Rydberg orbital, in an orbital centered on the

positive site of a zwitterion species, or in a dipole-bound orbital, one
must add to the bases containing valence, polarization, and
conventional diffuse functions yet another set of functions that are
extra diffuse. The exponents of these extra diffuse basis functions
can be obtained by scaling the conventional diffuse functions’
smallest exponent (e.g, by 1/3).
An example of a species needing extra diffuse basis functions-
Arginine anion
aug-cc-pVTZ, cc-pVQZ, pVDZ.
VDZ, VTZ, VQZ or V5Z specifies at what level the valence (V)
AOs are described. Nothing is said about the core orbitals because
each of them is described by a single contracted Gaussian type basis
orbital.
“cc” specifies that the orbital exponents and contraction coefficients
were determined by requiring the atomic energies computed using a
correlated method to agree to within some tolerance with
experimental data. If cc is missing, the AO exponents and
contraction coefficients were determined to make the Hartree-Fock
atomic state energies agree with experiment to some precision.
“p” specifies that polarization basis orbitals have been included in
the basis.
The number and kind of polarization functions differs depending
on what level (i.e., VDZ through V5Z) the valence orbitals are
treated.
For C at the VDZ level, one set of d polarization functions is
added. At the VTZ, two sets of d and one set of f polarization
functions are included. At the VTZ level, three d, two f, and one g
set of polarization functions are present, and at the V5Z, four d,
three f, two g and one h sets of polarization functions are
included. This strategy of building bases has proven especially
useful when carrying out complete-basis extrapolations.
“aug” specifies that (conventional) diffuse basis functions have
been added, but the number and kind depend on how the valence
basis is described. At the pVDZ level, one s, one p, and one d
diffuse function appear; at pVTZ a diffuse f function also is
present; at pVQZ a diffuse g set is also added; and at pV5Z a
diffuse h set is present.
6-31+G** or 3-21G*, 6-311+G*, or 6-31++G
3- or 6- specifies that the core orbitals are described in terms of a
single contracted Gaussian orbital having 3 or 6 terms.
–21 or –31 specifies that there are two valence basis functions of
each type (i.e., the valence basis is of double-zeta quality), one
being a contraction of 2 or 3 Gaussian orbitals and the other (the
more diffuse of the two) being a contraction of a single Gaussian
orbital.
–311 specifies that the valence orbitals are treated at the triple-zeta
level with the tightest contracted function being a combination of 3
Gaussian orbitals and the two looser functions being a single
Gaussian function.
* specifies that polarization functions have been included on the
atoms other than hydrogen; ** specifies that polarization functions
are included on all atoms, including the hydrogen atoms.
+ denotes that a single set of (conventional) diffuse valence basis
AOs have been included; ++ means that two such sets of diffuse
valence basis AOs are present.
Read the supplementary material to see if you understand that a
Carbon Aug-cc-pV5Z basis has 127 contracted basis functions built
from 209 primitive functions.
Keep in mind how things scale with the number of basis functions:
Calculating two-electron integrals<χa(1)χb(2)|1/r1,2 |χc(1)χd(2) > -M4
Solving the HF matrix eigenvalue equations for εk and ϕk- M3
Now that AO bases have been discussed, let’s return to discuss
how one includes electron correlation in a calculation.

There are many ways and each has certain advantages and
disadvantages.
Møller-Plesset perturbation (MPPT)- one uses the single-
configuration (usually single determinant) SCF process to determine
a set of spin-orbitals {φi}. Then, using H0 equal to the sum of the N
electrons’ Fock operators H0 = Σi=1,N F(i), perturbation theory is
used to determine the CI amplitudes for the CSFs. The amplitude for
the reference CSF Φ is taken as unity and the other CSFs'
amplitudes are determined by Rayleigh-Schrödinger perturbation
using H-H0 as the perturbation.

Advantages- Size extensive, no choices of “important” CSFs

needed, decent scaling at low order (M5 for MP2).

Disadvantages- Should not use if more than one determinant is

“important” because it assumes the reference CSF is dominant.
MP2 energy and first-order wave function expressions:

Ψ1 = - Σi<j(occ) Σm<n(virt) [< i,j | e2/r1,2 | m,n > -< i,j | e2/r1,2 | n,m >]

[ εm-εi +εn-εj]-1|Φi,jm,n >

E = <Φ|Η0 + V|Φ + Ψ1> = ESCF - Σi<j(occ) Σm<n(virt)

| < i,j | e2/r1,2 | m,n >

- < i,j | e2/r1,2 | n,m > |2/[ εm-εi +εn -εj ].

Single excitations do not contribute to the first-order wave function

(Brillouin theorem)
Two-electron integral transformation:

< φiφj | e2/r1,2 | φkφl > is what you need

Use φj = Σµ Cj,µ χµ and begin with < χiχj | e2/r1,2 | χkχl >
to form
< χiχj | e2/r1,2 | χkφm> = Σl Cm,l < χiχj | e2/r1,2 | χkχl >. M5
and then
< χiχj | e2/r1,2 | φmφm>, and < χiφj | e2/r1,2 | φkφl >,
and finally < φiφj | e2/r1,2 | φkφl > 4M5 total operation.
Multiconfigurational self-consistent field (MCSCF)- the
expectation value < Ψ | H | Ψ > / < Ψ | Ψ >, with Ψ being a
combination of determinental CSFs, is treated variationally and made
stationary with respect to variations in both the CI and the Cν,i
coefficients giving ΣJ HI,J CJ = E CI and a set of HF-like equations for
the Cν,I.
Advantages- can adequately describe bond cleavage, can give
compact (in CSF-space) description of ψ, can be size extensive if
CSF list is properly chosen, gives upper bound to energy.
Disadvantages- coupled orbital (Ci,µ) and CI optimization is a very
large dimensional optimization with many local minima, so
convergence is often a problem; unless the CSF list is large, not much
dynamical correlation is included.
Configuration interaction (CI)- the LCAO-MO coefficients of
all the spin-orbitals are determined first via a single-
configuration SCF calculation or an MCSCF calculation using a
small number of CSFs. The CI coefficients are subsequently
determined by making stationary the energy expectation value
< Ψ | H | Ψ > / < Ψ | Ψ > which gives ΣJ HI,J CJ = E CI .
Advantages- Energies give upper bounds and are variational (so
lower is better), one can obtain excited states from the CI matrix
eigenvalue problem.
Disadvantages- Must choose “important” CSFs, not size
extensive, scaling grows rapidly as the level of “excitations” in
CSFs increases (M5 for integral transformation; NC2 per
electronic state).
Coupled-Cluster Theory (CC)- one expresses the wave function
as

Ψ = exp(T) Φ,

where Φ is a single CSF (usually a single determinant) used in the

SCF process to generate a set of spin-orbitals. The operator T is
given in terms of operators that generate spin-orbital excitations

T = Σi,m tim m+ i + Σi,j,m,n ti,jm,n m+ n+ j i + ...,

Here m+ i denotes creation of an electron in spin-orbital φm and

removal of an electron from spin-orbital φi to generate a single
excitation. The operation m+ n+ j i represents a double excitation
from φi φj to φm φn.
When including in T only double excitations { m+ n+ j i}, the CC
wave function exp(T) Φ contains contributions from double,
quadruple, sextuple, etc. excited determinants:

exp(T) Φ = {1 + Σm,n,Iij tm,n,i,j m+ n+ j i + 1/2 (Σm,n,Iij tm,n,i,j m+ n+ j i)

( Σm,n,Iij tm,n,i,j m+ n+ j i)

+ 1/6 (Σm,n,Iij tm,n,i,j m+ n+ j i) (Σm,n,Iij tm,n,i,j m+ n+ j i)

(Σm,n,Iij tm,n,i,j m+ n+ j i) + …}Φ.

But note that the amplitudes of the higher excitations are given as
products of amplitudes of lower excitations (unlinked).
To obtain the equations of CC theory, one writes

H exp(T) Φ = Ε exp(T) Φ, then

exp(-T) H exp(T) Φ = Ε Φ, then

uses the Baker-Campbell-Hausdorf expansion

exp(-T) H exp(T) = H -[T,H] + 1/2 [[T,T,H]] - 1/6 [[[T,T,T,T,H]]] +. .

The equations one must solve for the t amplitudes are quartic:

< Φim | H + [H,T] + 1/2 [[H,T],T] + 1/6 [[[H,T],T],T] + 1/24

[[[[H,T],T],T],T] | Φ > = 0;

< Φi,jm,n |H + [H,T] + 1/2 [[H,T],T] + 1/6 [[[H,T],T],T] + 1/24

[[[[H,T],T],T],T] |Φ> =0;

< Φi,j,km,n,p|H + [H,T] + 1/2[[H,T],T] + 1/6 [[[H,T],T],T] + 1/24

[[[[H,T],T],T],T] |Φ> = 0,

The amplitudes of the double excitations that arise in the lowest

approximation are identical to those of MP2
ti,jm,n = - < i,j | e2/r1,2 | m,n >'/ [ εm-εi +εn -εj ].
Summary
1. Basis sets should be used that (i) are flexible in the valence region to allow for
the different radial extents of the neutral and anion’s orbitals, (ii) include
polarization functions to allow for good treatment of electron correlations, and
(iii) include extra diffuse functions if very weak electron binding is
anticipated. For high precision, it is useful to carry out basis set extrapolations
using results calculated with a range of basis sets (e.g., VDZ, VTZ, VQZ).

2. Electron correlation should be included because correlation energies are

significant (e.g., 0.5 eV per electron pair). Correlation allows the electrons to
avoid one another by forming polarized orbital pairs. There are many ways to
handle electron correlation (e.g., CI, MPn, CC, DFT, MCSCF).
3. Single determinant zeroth order wave functions may not be adequate if the spin
and space symmetry adapted wave function requires more than one
determinant. Open-shell singlet wave functions are the most common
examples for which a single determinant can not be employed. In such cases,
methods that assume dominance of a single determinant should be avoided.
4. The computational cost involved in various electronic structure calculations
scales in a highly non-linear fashion with the size of the AO basis, so careful
basis set choices must be made.
Special Tricks for calculating an anion’s energy when it lies above
that of the neutral?

Straightforward variational calculations will “collapse”

To produce a wave function and energy appropriate to
The neutral molecule plus a free electron with low energy.

Collapsed Anion Energies

A+B
E
A + B_

Accurate Anion Energies

RA-B
In the charge-scaling method, one fractionally
increases the nuclear charges on the atoms involved in
the bond, computes the anion-neutral energy
difference as a function of δq, and extrapolates to
δq → 0.

This is essential to do for species such as SO42- or

CO32-, which are not stable as isolated species.

It is also essential when studying σ*-attached states in,

for example, Cl3C-F + e- → Cl3C + F- dissociative
electron attachment or when attaching an electron to a
* orbital of benzene.
Consider calculating the Born-Oppenheimer energies of
various states of O2-. All three lowest states have bond
lengths where the anion is electronically unstable.
In the stabilization method one computes the anion-neutral energy difference
in a series of basis sets whose more diffuse basis functions’ exponents α are
scaled α →η α. Plotting the anion-neutral energy differences vs η produces
a stabilization plot that can be used to determine the metastable state’s energy.
4 Anion State Energy (eV)

E
Valence-localized state's energy
2

1
Resonance State Energy

Continuum states' energies

0
Orbital Scaling Parameter !
These energies grow with η because T scales as η2.
This method requires one to compute the energies of many anion states.
At certain η values, the diffuse basis functions can be combined
to describe the de Broglie λ of the asymptotic ψ and can match
ψ and d ψ /dr throughout.

Valence Radial Basis

Optimal Scaled
Diffuse Radial
Potential Basis

r
Radial Wavefunction
! too large

! too small
The lower-energy curves describe the dominantly-continuum solutions’
variation with η. When one of these solutions gains the proper de Broglie
and can match the energy of the valence-localized state, an avoided
crossing occurs. The energy of this crossing is the resonance energy.
4 Anion State Energy (eV)

E
Valence-localized state's energy
2

1
Resonance State Energy

Continuum states' energies

0
Orbital Scaling Parameter !