Martin Head-Gordon

Chemical Accuracy from Electronic

Structure Calculations 1
Modern Electronic Structure Methods:
2. Wavefunction-based theories.
MARTIN HEAD-GORDON,
Department of Chemistry, University of California, and
Chemical Sciences Division,
Lawrence Berkeley National Laboratory
Berkeley CA 94720, USA
This presentation will probably
involve audience
discussion, which will create
action items. Use
PowerPoint to keep track of
these action items during
your presentation
In Slide Show, click on the
right mouse button
Select Meeting Minder
Select the Action Items tab
Type in action items as they
come up
Click OK to dismiss this box
This will automatically create
an Action Item slide at the
end of your presentation
with your points entered.
Outline
1. Overview of standard wavefunction-based methods
Mean field, perturbation theory, coupled cluster theory
2. Fast methods for dynamic correlations
3. Fast methods for static correlations (with applications)
Complexity
Many-body Schrdinger equation is a partial differential
equation in 3n unknowns the positions of the electrons.
Exact (brute force) solution will scale approximately
exponentially with the number of electrons.
When done in a given basis of 1-electron expansion
functions, this is full configuration interaction (FCI).
Largest FCIs involve many-body expansions containing
billions of terms, for molecules with 2 or 3 of atoms (!).
Approximations are imperative. Accuracy vs feasibility.
Fundamental approximations
(1) The one-body problem (atomic orbital basis)
One-electron functions are expanded in a finite basis.
We use atom-centered functions with some maximum
angular momentum (L) on atoms of a given period.
Standardized basis sets of increasing L are available.
(2) The n-body problem (electron correlation method)
More about this in a minute mean field is the simplest
possible treatment.
Well-defined electronic structure models are completely
specified by these two approximations.
Branches of the family tree
Wavefunction-based electronic structure theory:
Minimize the energy by varying the wavefunction
Tremendously complicated unknown function:
Modeling the wavefunction yields model chemistries
Density functional theory
The unknown is very simple:
Hohenberg-Kohn theorem guarantees that:
True functional is unknown and probably unknowable
Modeling the functional gives DFT model chemistries.
Wave function approaches to electronic structure
Hartree-Fock (MO) theory (mean field): A
3
-A
4
cost.
HF: 99%! Correlation energy is roughly 1 eV per pair.
Perturbative treatment of the electron correlations: A
5
MP2: 80% of the last 1%. Pair correlations (doubles).
Treat single & double substitutions self-consistently: A
6
CCSD: 95% of the last 1%.
Correct for the triple substitutions perturbatively: A
7
CCSD(T): > 99% of the last 1%. Chemical accuracy.
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 2
State of the art thermochemistry
Up to 4 first row atoms: 0.3 kcal/mol
Martin, Taylor, Dunning, Helgaker, Klopper
include relativistic effects, core correlation, anharmonicity
CCSD(T) with extrapolation to the complete basis set limit
limited by (?): extrapolation, sometimes CCSD(T) itself.
Up to 10 or 15 first row atoms: 1 kcal/mol
Curtiss, Raghavachari, Pople: G2/G3 methods
CBS methods of Petersson et al
based on methods like CCSD(T), with either additivity
corrections for basis set, or extrapolation.
Mean field (Hartree-Fock) theory
Each electron moves in the average field of all others
Wavefunction is an antisymmetrized product of orbitals
The energy is minimized with respect to variations of the
orbitals, giving the Hartree-Fock equations:
F is the Fock operator, which is a one-particle Hamiltonian
with electron-repulsions averaged:
Chemistry with the Hartree-Fock model
Size-consistent; exact for isolated electrons
Absolute energies: about 99% of exact
The remaining energy is the correlation energy
Atomization energies: a disaster
The correlation energy of all separated electrons is the error
Thus roughly 1eV error (ie. 100 kJ mol) per separated pair
Optimized molecular geometries: are reasonable
Bondlengths are systematically slightly too short
Vibrational frequencies are systematically about 10% high
Dodging the correlation bullet
Look at energy changes where correlation errors cancel.
Energy change for isogyric reactions: a bit better
Isogyric reactions conserve the total number of electron pairs
Example: H
2
+ Cl !" HCl + H
Energy change for isodesmic reactions: much better!
Isodesmic reactions conserve numbers of each type of electron
pair (different bond pairs, lone pairs, etc)
Example: H
2
C=CH#CH
3
+ CH
4
!" H
2
C=CH
2
+ C
2
H
6
Schematic view of the Hamiltonian matrix
Group determinants by substitution level from HF
0, s(ingles), d(oubles), t(riples), etc
Look at the coupling of HF (0) to higher substitutions
HF orbitals are optimized to zero coupling to singles
HF wavefunction cannot couple to triples (or higher)
Double substitutions are the leading correction to HF
Double substitutions describe pair correlations
Correlated fluctuations of 2 electrons
In general, they are quartic in number.
Empty levels
Occupied levels
Correlations decay rapidly with separation.
Exponentially between functions of the same electron
Algebraically (R
-3
) between the 2 electrons ("dispersion)
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 3
Truncated configuration interaction
The exact (full CI) wavefunction is a linear superposition
As discussed before, this is exponentially expensive!
Truncated configuration for electron correlation?
For instance include all single + double substitutions (CISD)
Not useful because it fails to be size-consistent
CISD is exact for 1 H
2
molecule (only 2 electrons)
But is not exact for 2 H
2
s (that would require CISDTQ)
( ) 1 2 3

exact n HF
C C C C ! = + + + + "
Mller-Plesset perturbation theory
Apply perturbation theory to electronic structure, with:
The zero order Hamiltonian is the Fock operator:
Since the canonical MOs are eigenvectors of f:
Then the zero order problem is:
1
st
and 2
nd
order expansions
Expansion to first order gives the Hartree-Fock energy
Leading estimate of correlation energy is from 2
nd
order:
Recall from the Hamiltonian that the coupling is just doubles
MP2 correlation energy expression
The 2
nd
order expression in terms of doubles (d) is:
Or in terms of the occupied (i,j) and empty (a,b) orbitals:
MP2 correlation energy is negative definite
Computational complexity is 5
th
order in molecule size.
MP2 as a theoretical model chemistry
MP2 is size-consistent. Proper treatment of dispersion.
For closed shell molecules, accuracy is much better than HF
Bond lengths are (properly) longer, frequencies lower.
Directly reflects role of antibonding orbitals
For open shell molecules, accuracy is erratic.
MP2 is sensitive to spin-contamination in UHF.
Efficient computational algorithms are available:
Cubic memory (OVN), (larger) cubic disk space
Cost is asymptotically 5
th
order, but often lower.
Higher order corrections: MP3 and MP4
MP3 energy introduces coupling between doubles:
MP3 is A
6
in computational complexity
MP4 includes singles, doubles, triples & quadruples:
The MP4 triples are A
7
to compute
The triples energy contribution is chemically significant
MP3 and MP4 energies are erratic in accuracy
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 4
Coupled cluster theory
The exact wavefunction may be written as the
exponential of the sum of substitution operators
Truncation of T defines limited coupled cluster models.
Simplest models are CCD (doubles) and CCSD
Truncated CC models are clearly size-consistent: e.g.

!
exact
= exp

T
1
+

T
2
+

T
3
++

T
n
( )
"

!
CCSD
= exp

T
1
+

T
2
( )
"
Finding the coupled cluster amplitudes
Variational CC energy is intractable (can you see why?)
Instead we define CC equations by subspace projection:
Project with reference (0), and all substitutions in T, to give
enough equations to determine the energy, and T amplitudes.
This gives tractable equations, but the resulting energy
does not obey a variational bound.
The quality of the results will be best for cases when the single
reference is the dominant part of the wavefunction
( ) ( ) ( ) ( )

exp exp exp exp
VCC
E T H T T T = ! ! ! !

!
p

H " E
( )
exp

T
( )
! = 0
Energy & amplitude equations for CCSD
Project with the reference (0) to determine the energy
Project with singles (s) & doubles (d) for amplitudes:
Non-linear equations with A
6
computational complexity.
( ) ( )

exp 0 H E T ! " ! =

!E = "

H exp

T
( )
" = "

H 1+

T
1
+

T
2
+
1
2

T
1
2
( )
"

!
s

H 1+

T
1
+

T
2
+
1
2

T
1
2
+

T
1

T
2
+
1
3!

T
1
3
( )
!
C
= 0

!
d

H 1+

T
1
+

T
2
+
1
2

T
1
2
+

T
1

T
2
+
1
3!

T
1
3
+
1
2

T
2
2
+
1
2

T
1
2

T
2
+
1
4!

T
1
4
( )
!
C
= 0
CCSD as a theoretical model chemistry
Exact for a pair of electrons
Size-consistent
Analysis in terms of the Moller-Plesset series:
Complete in second order
Complete in third order
Missing the triples in fourth order
Contains some fifth and higher order terms
A significant improvement over MP2 for problem cases (radicals!)
Not much improvement for closed shell molecules
Occasional non-variational pathologies (E.g. stretched bonds)
Not chemically accurate for relative energies
Reflects the fact that the triples are chemically significant
Accounting for the triples in coupled cluster theory
CCSDT is the next complete truncation
Solving for the triples requires A
6
storage (vs A
4
for CCSD)
And A
8
computation (vs A
6
for CCSD).
Hence it is tempting to develop a perturbative approach
MP4 triples require A
4
storage and A
7
computation
The best standard method for this is CCSD(T)
First solve iteratively for CCSD (A
6
computation, A
4
storage)
Then non-iteratively account for triples (one A
7
step)
What is CCSD(T)?
Begin with the MP4 triples expression:
Simplest correction to CCSD is to just add it on
Better is to use converged doubles:
Best is to use converged S+D:
The energy expression is chosen to get the effect of singles on
triples correct in 5
th
order (cf doubles on triples in 4
th
order)
This defines the CCSD(T) method
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 5
CCSD(T) as a model chemistry
The gold standard of modern quantum chemistry for
closed shell molecules!
Structures accurate to about 0.002
Vibrational frequencies accurate to about 2%
Chemical accuracy (< 2 kcal/mol error) in bond energies!!)
Reducing the need for large basis sets
G2 and G3 methods use energy additivity approximations
Dunning cc-pVnZ basis sets can be extrapolated towards the
basis set limit using 2 point fits. Useful for energies.
Limitations of CCSD(T)
Can be erratic for radicals (reflects problems of HF)
Reduce by using optimized orbitals (or even KS orbitals)
Exhibits pathologies in bond-breaking
Do not use restricted orbitals when they are unstable!
Succeeds partly by cancellation of errors
Slightly overestimates the triples, which compensates for
neglect of connected quadruples.
Computational cost is high.
Wave function approaches to electronic structure
Hartree-Fock (MO) theory (mean field): A
3
-A
4
cost.
HF: 99%! Correlation energy is roughly 1 eV per pair.
Perturbative treatment of the electron correlations: A
5
MP2: 80% of the last 1%. Pair correlations (doubles).
Treat single & double substitutions self-consistently: A
6
CCSD: 95% of the last 1%.
Correct for the triple substitutions perturbatively: A
7
CCSD(T): > 99% of the last 1%. Chemical accuracy.
Outline
1. Overview of standard wavefunction-based methods
2. Fast methods for dynamic correlations
Auxiliary basis expansions
Local correlation methods
Opposite-spin second order correlation
Towards linear scaling
3. Fast methods for static correlations (with applications)
Auxiliary basis expansions
Early contributions by Whitten, Dunlap, Baerends, Almlof, others.
Auxiliary basis or resolution of the identity (RI) or density fitting (DF)
Popularized by Ahlrichs and co-workers for DFT, MP2.
Demonstrated efficiency
Developed standardized auxiliary basis sets: 3-4 times the AO basis size.
Replace 4-center integrals by (inexact) expansions:
Coefficients C minimize the Coulomb deviation of the fit:
ia jb ! C
ia
K
K L C
jb
L
KL
"
= B
ia
K
B
jb
K
K
"

C
ia
K
= ia L L K
!1
L
"

B
ia
K
= ia L L K
!1/ 2
L
"
Benefits of RI-MP2
The prefactor of the MP2 calculation is greatly reduced.
Far fewer integrals & better scaling with angular momentum.
The asymptotic rate-determining step is actually more expensive!
Timings for alanine tetrapeptides on 2 GHz G5. Tight cutoffs (12/9).
cc-pVDZ cc-pVQZ
------------------------------------------------------------------------
SCF (min) 55 714
------------------------------------------------------------------------
MP2 (min) 63 902
------------------------------------------------------------------------
RI-MP2 (min) 4 30
------------------------------------------------------------------------
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 6
Accuracy of the auxiliary basis expansions
(1) G2 dataset
Atomization energies of 148 neutral molecules: MP2 vs RI-MP2
RMS deviation = 0.10 kcal/mol
Max deviation = 0.27 kcal/mol
(2) Relative conformational energies of the tetrapeptide
27 conformations. Compare MP2 and RI-MP2
RMS deviation = 0.003 kcal/mol
Max deviation = 0.007 kcal/mol
Errors are chemically insignificant at the MP2 level.
Large molecule regime for RI-MP2
When SCF cost is smaller than RI-MP2.
M
2
(or less) cost for SCF vs M
5
for RI-MP2
DZP basis sets: 30-50 non-H atoms & up
depends on globularity, AO basis, cutoffs, etc.
cc-pVDZ alanine polypeptides;
2GHz G5; tight cutoffs (12/9).
tetrapeptide octapeptide hexadecapeptide
--------------------------------------------------------------------------------------------------
SCF (min) 55 390 1779
--------------------------------------------------------------------------------------------------
RI-MP2 (min) 4 120 3935
--------------------------------------------------------------------------------------------------
What is a local correlation model?
A local correlation model approximates a standard
correlation method (well look at MP2), computing only
some of the amplitudes, based on spatial truncation.
Convention wavefunction-based models. One can
(1) Make models for the correlation energy, stopping at some
order in excitations (CCSD) or perturbation theory (MP2).
(2) Or, could select retained excitations based on some cutoff
criterion such methods are today not widely used.
A similar choice will confront us in local correlation
Local correlation methods as model chemistries
Local correlation models involve strong approximations.
Often recover 98 or 99% of the correlation energy
Local correlation models thus define modified
theoretical model chemistries that must be tested.
They should have the features of a good theoretical
model chemistry (not all do, however!).
Efficient (a good tradeoff between accuracy and cost)
Size-consistent
No molecule-dependent adjustable parameters
Continuous potential energy surfaces, etc.
Pulay-Szabo Local Electron Correlation Method
First and most widely used approach to reduce the
unphysically high cost of wavefunction-based methods:
[Link], Chem. Phys. Lett. 100, 151 (1983).
[Link] & [Link], Ann. Rev. Phys. Chem. 44, 213 (1993)
Successfully adopted and extended by other groups:
Murphy, Friesner et al, JCP 103, 1481 (1995).
Reynolds, Martinez, Carter, JCP 105, 6455 (1996).
Schutz, Hetzer, Werner, JCP 111, 5691 (1999).
Pulay-Szabo local model for double substitutions
(1) Localize the occupied orbitals
symmetry-equivalent atoms may not be equivalent: benzene
(2) Project atomic orbitals into the virtual space
(3) Allow only substitutions from a pair of occupied
orbitals into all pairs of virtual orbitals that strongly
overlap one or the other of the occupied orbitals.
potential energy surfaces will not be strictly continuous.
Computational complexity is greatly reduced.
98% (or more) of correlation energy is recovered.
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 7
Status of the Pulay-Saebo model
Linear scaling achieved with the work of Schutz/Werner
Long-range dispersion truncated
Extension to CCSD and triples completed.
CCSD is used for strong pairs, MP2 for weak pairs
Very weak pairs are further approximated or truncated.
Some drawbacks for chemical applications due to:
Discontinuous potential energy surfaces
Use of adjustable numerical parameters
Failure for delocalized electronic structure
Atomic local correlation models
We re-examined the 2 steps to make local correlation models:
(1) Choose suitable local 1 particle representations
Make them atom centered!
A minimal atom-centered basis for occupied space
Projected atomic orbitals for the virtual space
(2) Truncate the n-particle spaces in a physical way.
Satisfy model chemistry criteria with atomic truncations
Parameter free models, continuous potential surfaces
TRIM local correlation model
Exact (tetra-atomics in molecules)
Tri-atomics in molecules (TRIM)
Diatomics in molecules (DIM)
Atoms in molecules (AIM)
Fractional correlation recovery in linear alkanes
6-31G* basis.
TRIM yields
99.7% of
correlation energy
DIM yields 95%
of correlation
energy
TRIM
DIM
Torsional barriers via local MP2 methods
12 conjugated ethylene and benzene derivatives
R-OH, R-CHO, R-C
2
H
3
, R-CFO, R-NH
2
, RNO
2
Delocalization and barriers are tests for local correlation
6-311G** basis, frozen core.
Deviations from untruncated MP2 (in kcal/mol)
(magnitude of barriers ranges from 1 to 8 kcal/mol)
TRIM-MP2 energy
The working expression for the TRIM energy is:
where the DIM energy itself is:
Indices i,a represent canonical occupieds and virtuals
Indices P,Q are recanonicalized atomic single substitutions, with
corresponding excitation energies E
P
TRIM reduces the number of doubles to cubic.
Rate-determining steps are 4th order, not 5th order
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 8
Auxiliary basis TRIM-MP2 (Rob Distasio Jr)
RI-TRIM greatly reduces the TRIM prefactor.
Design goals: quadratic memory, cubic disk; cubic I/O
Sketch of algorithm: N
ASIS
~N
aux
~3n
AO
cost
(1) Prepare the EPAOs, recanonicalization transform: N
3
(2) Prepare B coeffs in MO and AO reps; save: N
3
+ ovN
2
+ ovnN
(3) Make TRIM (& DIM) K integrals; save: 2ovN
2
(4) Make TRIM (& DIM) J integrals; add to K: 2ovN
2
Speedup vs RI-MP2: about o/30
Auxiliary basis TRIM-MP2 (Rob Distasio Jr)
Quadratic memory; cubic disk; cubic I/O
Implemented by Rob Distasio Jr.
cc-pVDZ alanine polypeptides; 2GHz G5; tight cutoffs (12/9).
tetrapeptide octapeptide hexadecapeptide
-------------------------------------------------------------------------------------------
SCF (min) 55 390 1779
-------------------------------------------------------------------------------------------
RI-MP2 (min) 4 120 3935
-------------------------------------------------------------------------------------------
RI-LMP2 (min) 3 53 850
-------------------------------------------------------------------------------------------
TRIM accuracy (Rob Distasio Jr)
Compare against MP2 for 27 tetrapeptide conformations
RMS (kcal/mol) Max (kcal/mol)
------------------------------------------------------------------------
cc-pVDZ 0.19 0.43
------------------------------------------------------------------------
cc-pVTZ 0.16 0.36
------------------------------------------------------------------------
cc-pVQZ 0.09 0.20
------------------------------------------------------------------------
(TQ)extrap 0.04 0.08
------------------------------------------------------------------------
Atomic truncations for triples (Paul Maslen)
Exact (6 atoms coupled together)
A
7
cost scaling
Singly ionic (4 atom clusters)
A
5
cost scaling
Covalent (3 atom clusters)
A
4
cost scaling
Fractional triples energy recovery (Paul Maslen)
Poly-ynes, 6-31G* basis, frozen core, MP4 triples.
50
60
70
80
90
100
2 4 6 8
poly-yne chain length
%

t
r
i
p
l
e
s

r
e
c
o
v
e
r
y
Singly ionic
Covalent
Scattering picture of triple substitutions
Why is covalent triples truncation much worse (< 70%) than it
was for doubles (>90%)?
No 3-body terms in H, so triples result from 2 binary collisions
Second electron is scattered twice
(before antisymmetrizing)
Hence less local than doubles
Local triples models should take this into account
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 9
Singly antisymmetrized triples energy I
Antisymmetrize once or twice?
To preserve the scattering picture of triples once is plenty!
Non-locality of the second pair (bj) is now preserved in the bra
E
(T )
= !
1
3
T2" J { }
aibjck

A
3
T2" J { }
aibjck
D
aibjck aibjck
#
E =

A!

H

A! = !

H

A

A! = !

H

A!
T2! J { }
aibjck
= T
aidj
J
dbck
" T
aibl
J
ljck
l
#
d
#
Singly antisymmetrized triples energy II
To change basis, use Almlof/Haser Laplace transform approach:
Do some further rearrangements to get:
The bra can be kept non-local in the doubly scattered variables (bj)
The ket carries the rest of the antisymmetrizer.
We can express everything in the projected minimal and projected
AO representations and make the local models
1
x
= exp !"x ( )
0
#
$
dx % w
L
exp !"
L
x ( )
L
N
L
&
E
(T )
= !
1
3
T2
L
" J
L
{ }
aibjck
L
N
L
#
$
%
&
'
(
)

A
3
T2" J { }
aibjck
aibjck
#
Spin components of the MP2 energy
Opposite spin (OS) correlation is alpha-beta
No OS correlation at all in mean field
Largest component of short-range correlation
Basis set convergence is slow (L
-4
)
Same spin (SS) correlation is alpha-alpha and beta-beta
Correlation from Fermi statistics is included in mean-field
Smaller in magnitude
Basis set convergence is slightly less slow
E
MP2
= E
SS
+ E
OS
Scaling the MP2 spin components
a = 1/3 (small part gets smaller)
b = 6/5 (dominant part gets larger)
Results are significantly improved over MP2!
RMS error (vs QCISD(T)) for 41 mostly isogyric reaction energies
MP2: 4.4 kcal/mol
SCS-MP2: 2.2 kcal/mol
See Stefan Grimme, J. Chem. Phys. 2003, 118, 9095-9102.
E
SCS!MP2
= aE
SS
+ bE
OS
Scaled opposite spin (SOS) MP2 correlation?
2 reasons to examine this simpler variant:
(1) To see how essential it is to employ 2 parameters
(2) If not essential, no SS correlation makes fast methods simpler!
(The SS term contains the indirect exchange contributions)
What should the scaling factor be?
Lower bound ~ 1 (simply neglect SS)
Best guess ~ 6/5 + (1/3)(1/3) ~ 1.3 (mimics SCS at short-range)
Upper bound = 2 (correct at long-range)
E
SOS!MP2
= cE
OS
Tests of the scaling factor (Rohini Lochan)
For the same test set of reaction energies
Conclusions:
(1) Even no scaling is better than MP2 (remarkable!)
(2) Scale factor of 1.2 or 1.3 looks best
(3) Quality of results is close to SCS-MP2
We adopt the value of 1.3
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 10
Fast evaluation of the opposite spin MP2 energy
Using an auxiliary basis, the SOS-MP2 energy can be evaluated
with no 5th order steps.
This is done using the Laplace formulation in an auxiliary basis.
No other approximation (e.g. local model) is required.
The working equations are: (with rate-determining step ovN
2
)

E
MP2
OS
= !
1
2
w
z
IA| JB ( )
2
AB
vir
"
IJ
occ
"
z
"
= !
1
2
w
z
Y
KM
K L Y
LN
M N
MN
aux
"
KL
aux
"
z
"

Y
KL
= C
ia
K
C
ia
L
ia
!
Timing comparison (Yousung Jung)
Timings in minutes (375 MHz IBM Power3)
Basis is 6-31G*; all electrons correlated.
4th order algorithm wins cleanly for bigger systems.
Auxiliary basis (N
5
) Laplace (N
4
) Error (a.u.)
C20H42 15 19 0.000000
C30H62 89 83 0.000002
C40H82 364 243 0.000003
C50H102 977 569 0.000004

ArNe dimer (CP-cc-pVQZ) (Rohini Lochan) Modified opposite spin scaling
Issue: At long range, same-spin and opposite spin contributions to
the MP2 energy are equal.
But SOS scales by only 1.3, not 2.
Solution: Define modified integrals for opposite spin correlation:
Choose K to get the long-range right.
The free parameter ! can then be adjusted.

ia jb ( )
!
= dr d " r #
i
r ( ) $ $
#
a
r ( )
1+ Kerf ! r % " r ( )
r % " r
&
'
(
(
)
*
+
+
#
j
" r ( )#
b
" r ( ).

K = 2 !1
77 atomizations as a function of ! (Rohini Lochan)

30 barrier heights as a function of ! (Rohini Lochan)

Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 11
ArNe dimer again (CP-cc-pVQZ) (Rohini Lochan) Numerical data for alkane chains (Yousung Jung)
g r
1
, r
2
( ) = ! r
1
" r
2
( )
C50
1.0E-13
1.0E-11
1.0E-09
1.0E-07
1.0E-05
1.0E-03
1.0E-01
1.0E+01
0 10 20 30 40 50
Position of carbons
F
i
t

c
o
e
f
f
i
c
i
e
n
t
s
Coulomb
Overlap
Very slow long-range
decay in the Coulomb fit!
Initial decay is rapid for
Coulomb & overlap fits

C
ia
K
= ia L
g
L K
g
!1
L
"
g r
1
, r
2
( ) = r
1
! r
2
!1
Ahlrichs VDZ(d) basis; fairly tight cutoffs (10/8)
Plot max fit cofficient for each C for ss distribution on C
1
Long-range Coulomb fit coeffs decay algebraically
This is clear when plotted directly rather than logarithmically
Is this unphysical? Should non-overlapping functions play a role?
Yet chemical results with the Coulomb metric are better.
Long range decay of Coulomb fit (C50)
y = 0.0001x
-1.2429
R
2
= 0.9646
1.0E-06
1.0E-05
0 10 20 30 40 50
Position of carbons
A 2-site model (Alex Sodt)
Fit an s-type Gaussian with q=1 using 2 auxiliary basis functions.
Both have q=1 also, and are s functions.
1st is on the same site but has wrong exponent; 2nd is on a remote site.
The Coulomb fitting equations become:
with solutions:
In the long-range limit (r large):
On-site Gaussian does not contain correct charge
Off-site Gaussian partly compensates for this
A non-local correction to partly fix a local deficiency
i r
!1
r
!1
i
"
#
$
%
&
'
c
1
c
2
"
#
$
%
&
'
=
j
r
!1
"
#
$
%
&
'
c
1
= ijr
2
!1 ( ) / i
2
r
2
!1 ( )
c
2
= r(i ! j) / i
2
r
2
!1 ( )
c
2
! (i " j) / i
2
r ( ) c
1
!i / j
Long-range decay is affected by spatial packing
Analytical theory: Peter Gill (to be published)
Numerical experiments: Alex Sodt
2-d graphitic sheets 3-d diamond chunks
Attenuated Coulomb operator for fitting
Want good chemistry of Coulomb fit and fast decay of overlap fit.
Fit operator is Coulombic at short range and zero at long range?
The attenuated Coulomb operator:
First look at decay as a function of attenuation (!)
g r
1
, r
2
( ) = erfc ! r
1
" r
2 ( ) / r
1
" r
2
w = 0.01
1.00E-13
1.00E-11
1.00E-09
1.00E-07
1.00E-05
1.00E-03
1.00E-01
0 10 20 30 40 50
Position of carbons
F
i
t

c
o
e
f
f
i
c
i
e
n
t
s
Attenuated Coulomb
Overlap
w = 0.04
1.00E-13
1.00E-11
1.00E-09
1.00E-07
1.00E-05
1.00E-03
1.00E-01
0 10 20 30 40 50
Position of carbons
F
i
t

c
o
e
f
f
i
c
i
e
n
t
s
Attenuated Coulomb
Overlap
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 12
Effect of attenuation on quality of results
Atomization energies of 148 molecules from the G2 database.
Deviations from conventional MP2 calculations; VDZ(d) basis
Small ! values give the Coulomb metric chemistry (best)
Large ! values give the overlap metric chemistry (worse), but locality.
The largest ! value that does not degrade the chemistry is 0.2 or so.
0.01 0.1 0.2 0.3 0.4 1 100
a
MAE 0.08 0.08 0.08 0.11 0.14 0.27 0.52
b
MAXE 0.27 0.28 0.39 0.53 0.62 1.01 2.24
c
RMS 0.10 0.09 0.11 0.16 0.19 0.35 0.68

Overlap-like Coulomb-like
Sparsity in 3-center integrals
Systems are simple linear alkanes, with the VDZ(d) basis.
Compare 3 metrics (!=0.1 for attenuated), with 2 tolerances.
Quadratic growth is evident with the Coulomb metric
(mn|K): Cutoff = 1.0E-6
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
(mn|K): Cutoff = 1.0E-9
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
10
#9
10
#6
Sparsity in the fit coefficients, C.
The fit coefficients are dimensionless:
Linear growth is evident at low precision for all metrics.
At high precision, Coulomb metric is not linear.
C(mn,K): Cutoff = 1.0E-6
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
C(mn,K): Cutoff = 1.0E-9
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
10
#6
10
#9

C
ia
K
= ia L
g
L K
g
!1
L
"

ia jb ! C
ia
K
K L C
jb
L
KL
"
Sparsity in the 3-center B factors.
B factors are defined so their product yields integrals:
Thus the B factors contain some long-range terms.
This accounts for quadratic scaling at high precision. There is still useful
some useful sparsity in the B factors at low precision.
B(mn,K): Cutoff = 1.0E-9
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
B(mn,K): Cutoff = 1.0E-6
0.0E+00
1.0E+08
2.0E+08
3.0E+08
4.0E+08
5.0E+08
0 10 20 30 40 50 60 70
Chain length
N
u
m
b
e
r

o
f

s
i
g
n
i
f
i
c
a
n
t

e
l
e
m
e
n
t
s
Coulomb
Attenuated Coulomb
Overlap
10
#6
10
#9

B
ia
K
= C
ia
L
L K
1/ 2
L
!

ia jb ! B
ia
K
B
jb
K
K
"
Exploiting sparsity in SOS-MP2 (Yousung Jung)
Forming the Y matrix (at each quadrature point) is the bottleneck:
Employ attenuated Coulomb metric (with !=0.1) for sparse C.
Apply sparsity to evaluate Y. Exact code: no local model.
Also employ sparsity to perform transformations for C

Y
KL
= C
ia
K
C
ia
L
ia
!
Timings for a full code (Yousung Jung)
Timings (minutes) on a 2 GHz Apple xserve (1 cpu).
VDZ(d) basis.
SCF threshold of 10
-8
. Linear scaling algorithms enabled.
Drop tolerance of 10
-6
for evaluation of sparse C, Y,
C30H62 C60H122 C90H182
t (SCF) 22 74 130
t (RI-MP2) 20 463 3052
t (SOS-MP2) 20 235 1033
t (local SOS) 13 60 160

Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 13
Outline
1. Overview of standard wavefunction-based methods
2. Fast methods for dynamic correlations
3. Fast methods for static correlations
Lessons from variational coupled cluster calculations
Perfect and imperfect pairing
The worlds first indefinitely stable singlet diradical?
The phenalenyl dimer
Variational coupled cluster doubles
Troy Van Voorhis
Failures of restricted CCD and CCSD are well-known
for bond-breaking problems.
Is the problem mainly neglect of higher excitations?
Is the problem nonvariational solution for the amplitudes
Or, perhaps more likely, are both these issues fatal?
Performing variational coupled cluster doubles (VCCD)
calculations will answer these questions.
VCCD energy is an expansion up to n-fold excitations.
Hence it is a restricted form of full CI (FCI) wavefunction.
Perform benchmark calculations using a special FCI code.
Water double dissociation (Troy Van Voorhis)
CCD
FCI
VCCD
Are
doubles
enough?
STO-3G
Water double dissociation (Troy Van Voorhis)
CCD
FCI
VCCD
VDZ basis
The larger VCCD
error in the VDZ basis shows
the value of using a valence
active space...
STO-3G N
2
dissociation (Troy Van Voorhis)
CCD
FCI
VCCD
How
important are
quadruples and
hextuples?
VDZ basis N
2
dissociation (Troy Van Voorhis)
CCD
FCI
VCCD
VCCD/VDZ
shows much larger errors
than VCCD/STO-3G
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 14
Main benchmark conclusions (Troy Van Voorhis)
(1) VCCD yields dramatically improved results relative to
CCD for 4 and 6 electron bond-breaking problems.
The main problem with restricted CCD, CCSD in bond-
breaking is nonvariational failure.
Searching for novel doubles-based coupled cluster
methods for bond-breaking is potentially fruitful.
(2) VCCD performs noticeably better in the minimal basis
than in the VDZ basis, relative to FCI.
Traditional ideas about partitioning the correlation energy
into two parts are going to be important.
Quadratic coupled cluster method
Troy Van Voorhis
It is clearly desirable to improve CCD in a way which
puts in more physics.
Look for an intermediate step towards VCCD which
does not greatly increase the computational cost.
Previous efforts in this direction include:
extended coupled cluster (ECC) method of Arponen
UCC and XCC methods of Bartlett and co-workers
suggestions of Kutzelnigg
Usual coupled cluster energy can be written as a
variational principle:
where we have defined:
We shall investigate a quadratic generalization:
which will be minimized with respect to T, $ amplitudes
Quadratic coupled cluster energy functional
Troy Van Voorhis
Water double dissociation (Troy Van Voorhis)
QCCD and VCCD
curves coincide
STO-3G
CCD fails
STO-3G N
2
dissociation (Troy Van Voorhis)
CCD
VCCD
QCCD
QCCD is far
better than CCD, but
there is noticeable
error relative to
VCCD
Valence space electron correlation
Divide the correlation problem into two parts:
Low energy, long-wavelength static correlations
Associated with near-degeneracies between orbitals (bond-
breaking).
Important for qualitative description of molecular
chemistry and potential surfaces
High energy, short wavelength dynamic correlations:
Associated with atomic-like correlations
Dynamic correlation is important for quantitative
prediction of reaction energies. Also dispersion forces.
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 15
Defining the valence correlation
(1) full valence active space
N(active) = M (valence minimal basis)
Roos, Ruedenberg
(2) perfect pairing active space
N(active) = n (number of valence electrons)
One correlating orbital, j*, for each occupied orbital, j
Bonds and lone pairs are correlated.
Electron pair becomes the basic quantity.
2 Valence LCC models for pair correlations
(1) Simplest model is perfect pairing
Only intrapair correlations, only a linear number of variables
Hurley, Pople, Lennard-Jones model (1950s)
Popularized by Goddard as GVB-PP
Coupled cluster version first explored by Cullen (1996)
Exact for the case of 1 electron pair!
1. Intrapair correlations
perfect pairing
2 Valence LCC models for pair correlations
(2) Next model must be imperfect pairing!
Includes interpair correlations that dont transfer electrons.
A CC analog of GVB-RCI (but size-consistent!)
T. Van Voorhis and MHG, Chem. Phys. Lett. 317, 575 (2000).
T. Van Voorhis and MHG, J. Chem. Phys. 115, 7814 (2001).
2. +Covalent interpair correlations
imperfect pairing
Efficient algorithms for valence LCC models
Linear or quadratic number of variables.
M
3
computation. Great improvement over M
6
!!
Outer loop is over orbital iterations (until converged).
Make 6%O Coulomb and exchange matrices (O=# of valence pairs) this is the
rate-determining step cubic computation.
From these matrices build required 2-electron integrals
Solve the amplitude equations (essentially free at present).
Form the orbital gradient and update orbitals.
T. Van Voorhis, MHG, J. Chem. Phys. 117, 9190 (2002).
Alex Sodt, Greg Beran, MHG (in preparation).
The auxiliary basis implementation!!
What is a singlet diradical?
Like aromaticity, it is not defined by a single quantity
Experimental perspective:
Low singlet-triplet gap
High reactivity
Low-lying electronic excitations (e.g. infra-red)
Wavefunction perspective:
Highest occupied MO (HOMO) and lowest unoccupied MO (LUMO) are
almost equal in energy
Instead of the HOMO having 2 electrons and the LUMO 0, we approach
having 1 electron in each antiferromagnetically coupled.
This means there are strong correlations in the wavefunction
Occupation number scale for diradical character
Measure diradical character by 100%n(LUMO)%
n(LUMO) is the occupation number of the (nominal) LUMO
Example of H2 dissociation:
At equilibrium, n(LUMO) & 0.05 (i.e. 5% diradical)
Wavefunction is well-described by 1 determinant
Towards dissociation, n(LUMO) " 1 (i.e. 100% diradical)
Wavefunction consists of 2 determinants
Characterize different diradicals on this common scale.
What wavefunctions can we use to get occupation numbers?
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 16
A new stable singlet diradical
Scheschkewitz et al Science, 2002, 295, 1880-1881
Planar BPBP ring
Stable at room temp.
B-B distance of 2.6
No ESR signal.
Compare against Si(100) surface
Use Si
9
H
12
cluster to model a dimer on Si(100).
48 valence electrons (24 pairs).
We optimized the structure to test whether the dimer buckles (it does not).
HOMO: 1.68 electrons LUMO: 0.32 electrons
How diradicaloid is this molecule? (Yousung Jung)
Valence active space coupled cluster calculations
94 active electrons (47 pairs), 6-31G* basis
HOMO: 1.83 electrons LUMO: 0.17 electrons
Stability comes from reduced (17%) diradical character
Y. Jung and MHG, ChemPhysChem 4, 522 (2003)
Origin of the stability (Yousung Jung)
3 step orbital interaction view:
1. Atomic pz orbitals on the two B atoms interact to make:
SYM (bonding) orbital.
ASYM (antibonding) orbital.
2. Coupling between SYM and lower P-H '-bonding orbitals
This will destabilize SYM and decrease the gap.
3. Coupling between SYM and higher P-H '* anti-bonding levels
This will stabilize SYM and increase the gap.
Which interaction is strongest and controls the chemistry?
1. Effect of the initial B-B splitting (Yousung Jung)
Compare (BPBP) ring with (BNBN) analog.
For simplicity replace methyl groups by hydrogens
Shorter B-B distance in (BNBN)
2.04 vs 2.60
Should imply a larger gap and less diradical character
Observe the opposite!
44% diradical character vs 22%.
The more diradicaloid compound has the shorter bondlength!
Conclude that initial splitting is substantially perturbed by
interactions with the neighboring groups.
2. Coupling to " bonding orbitals (Yousung Jung)
Investigate the (BSBS) compound: no valence antibonding
orbitals of appropriate symmetry.
Coupling will be to filled lone pairs on the S atoms only.
Ordering reverses! Bonding (SYM) level becomes the LUMO!
HOMO (1.68 electrons) LUMO (0.32 electrons)
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 17
3. All interactions (revisited) (Yousung Jung)
Coupling to both bonding and antibonding levels is
evident graphically in the HOMO see the lack of
amplitude at the phosphorus atoms.
(along HBBH axis) (above HBBH axis)
3. Coupling to bonding & antibonding levels
BPBP ring BNBN ring
Phenalenyl dimer: a strong # stacking complex
C
13
H
9

: a very stable radical in

solution and in the solid state.
C
13
H
9

dimerizes to form a
stable ( stacking complex
A crystal structure of the
dimer of the tri-t-Bu
derivative has been obtained
experimentally.
Phenalenyl dimer: the 12 center 2-electron bond?
HOMO LUMO
Phenalenyl dimer: computational approaches
Previous theory: Y. Takano, H. Isobe, T. Kubo, Y. Morita, K. Yamamoto,
[Link], T. Takui and K. Yamaguchi, J. Am. Chem. Soc. 124, 11122-
11130 (2002).
Restricted HF, DFT methods inapplicable
Unrestricted DFT unbound
Unrestricted HF enormous spin contamination
So: we decided to use perfect pairing (with just 1 pair),
and then correct via 2
nd
order perturbation theory.
DFT potential surface (B3LYP/6-31G*)
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
2 4 6 8 10 12 14 16 18 20
Separation ()
R
e
l
a
t
i
v
e

E
n
e
r
g
y

(
k
c
a
l
/
m
o
l
)

Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 18
Phenalenyl dimer potential surface (Yousung Jung)
-30.0
-20.0
-10.0
0.0
10.0
20.0
30.0
40.0
50.0
1.0 3.0 5.0 7.0 9.0
R ()
R
e
l
a
t
i
v
e

e
n
e
r
g
y

(
k
c
a
l
/
m
o
l
)PP(1)
CP-MRMP2
MRMP2
MRMP2 results (Yousung Jung)
PP(1) calculations show no minimum
Covalency by itself is not responsible for net binding
PP(1) + 2
nd
order perturbation theory finally, a bond!!
Equilibrium distance (calc) 3.1 Angstroms
(expt) 3.1 Angstroms
Binding energy (theory using counterpoise correction)
(calc) -10.8 kcal/mol
(expt) -8.8 kcal/mol
-9.5 kcal/mol
Character of the chemical bond (Yousung Jung)
Diradicaloid character:
LUMO occupation number (ON) at R=3.1A: 0.25 electrons
i.e. 25% diradicaloid character.
Dispersion-assisted 12-center 2-electron bond:
Covalency alone cannot overcome repulsions to give net binding.
This weak and diradicaloid chemical bond is assisted by dispersion
This interplay between weak covalency and dispersion makes the bond
in the phenalenyl dimer distinctive relative to weaker ( complexes
0.95 0.32 0.25 0.15 ON(LUMO)
5.0 3.3 3.1 2.6 R ()
Phenalenyl dimer cation
Compare the 12-center 2-electron bond (neutral) against the 12-
center 1-electron bond (cation)
Experimental dimer binding energies:
J. Kochi, unpublished (private communication).
Neutral: -8.8 kcal/mol (ES) 9.5 kcal/mol (ESR)
Cation: -5.3 kcal/mol (ES) -4.1 kcal/mol (ESR)
2-electron bond roughly twice as strong as the 1-electron bond
Theoretical methods:
Use ROHF based MP2 (RMP2) with counterpoise, for cation.
This is approximately comparable to PP(1) + MP2 for neutral.
DFT surface for cation (B3LYP/6-31G*)
Similar to
H
2
+
Self-
interaction
problem
-16.5
-15.5
-14.5
-13.5
-12.5
-11.5
-10.5
-9.5
-8.5 2 4 6 8 10 12 14 16 18 20
Separation ()
R
e
l
a
t
i
v
e

E
n
e
r
g
y

(
k
c
a
l
/
m
o
l
)

Phenalenyl dimer potential surface (David Small)
Martin Head-Gordon
Chemical Accuracy from Electronic
Structure Calculations 19
Features of the 12-center 1-electron bond
The cation dimer is bound at the ROHF level
Recall the neutral was unbound at the corresponding PP level!
We have lost roughly half of the covalent interaction
But gained electrostatics: the effect of charge-quadrupole interactions and induced
moments (cf. Mg-benzene, etc)
Conclusion: electrostatics in the cation outweighs covalency!
With dispersion, the cation dimer is more strongly bound than the neutral by
roughly a factor of 2.
Consistent with the effect of dispersion being similar in both systems.
Inconsistent with experiment: cation dimer is less strongly bound by roughly a
factor of 2.
Reconciling theory and experiment
Experiments are in solution, theory is gas phase.
Solvent is dichloromethane ($
r
! 9 ).
Simplest possible solvation treatment is the Born model
A cation dimer radius of 5 gives
The monomer is smaller and thus more strongly solvated.
If we assume monomer solvation to be roughly half that of the
dimer, we can reconcile theory with experiment.
Need experimental gas phase studies of the dimer cation!
Sigma dimerization of phenalenyl (David Small)
A short summary of another fascinating channel:
Normal sigma C-C bond strength ~ 84 kcal/mol
Calculated sigma dimerization energy ~ 16 kcal/mol
Why so weak? Bond formation vs loss of resonance.
Implies sp
2
" sp
3
conversion costs ~ 34 kcal/mol
Summary: dispersion assisted diradicaloid bonds
Stacking interactions are a fascinating interplay between:
Filled orbital Coulomb and exchange repulsions
Normally prevent closer approach than Van der Waals
Covalent bond and antibond interactions of 2 radical electrons
Provides energy bonus for sub-VdW distances
But is not enough by itself
Dispersion interactions
Also favor shorter distances, but also not enough by itself!
Electrostatics
1-electron bond beats the 2-electron one in the phenalenyl dimer!
See [Link]. 126, 13850 (2004).
(Grateful) Acknowledgements.
To all students and postdocs, past and present, in my group.
Rohini Lochan is their representative here!
Without them, little of this research would have been accomplished.
Funding: Dept of Energy, NSF,
NIH (SBIR program), AFOSR
In collaboration with Q-Chem
(MHG is a part-owner of Q-Chem)