Axel Becke

Report
A Density-Functional Theory for
Covalent and Noncovalent
Chemistry
Non-empirical and fast 
• Review of the XDM (exchange-hole dipole moment)
dispersion model of Becke and Johnson
• Combining the model with non-empirical exchange-correlation
GGAs (Kannemann and Becke)
• Tests on standard bio-organic benchmark sets
Dispersion Interactions from the
Exchange-Hole Dipole Moment:
The “XDM” model
Axel D. Becke and Erin R. Johnson*
Email: [email protected]
Department of Chemistry, Dalhousie University,
Halifax, Nova Scotia, Canada
(*now at University of California, Merced)
What is the source of the “instantaneous” multipole moments
that generate the dispersion interaction?
Suggestion:
Becke and Johnson, J. Chem. Phys. 122, 154104 (2005)
> Becke and Johnson, J. Chem. Phys. 127, 154108 (2007) <
The dipole moment of the exchange hole!
(Position, rather than time, dependent)
In Hartree-Fock theory, the total energy of many-electron system
is given by:
EHF
1
1  (r1 )  (r2 ) 3 3
2
3
3
    i   i d r   Vnuc d r  
d r1d r2  E X
2 i ,
2
r12
    
EX  EX  EX
E X
N
   i2
i 1
 i (r1 ) j (r1 ) i (r2 ) j (r2 ) 3 3
1
   
d r1d r2
2 ij
r12
The “exchange” energy
The Exchange Hole
The ( -spin) exchange energy can be rewritten as follows:
E X
where
hX (r1 , r2 ) 3 3
1
   (r1 )
d r2 d r1
2
r12
hX (r1 , r2 )  
1
 i (r1 ) j (r1 ) i (r2 ) j (r2 )

 (r1 ) ij
is called the “exchange hole”.
Physical interpretation: each electron interacts with a “hole” whose shape
(in terms of r2) depends on the electron’s position r1.
When an electron is at r1, the hole measures the depletion of probability,
with respect to the total electron density, of finding another electron of the
same spin at r2. This arises from exchange antisymmetry.
It is simple to prove that…
• The hole is always negative.
• The probability of finding another same-spin electron at r2=r1 (“on
top” of the reference electron) is completely extinguished:
hX (r1 , r1 )   (r1 )
Pauli or exchange “repulsion”!
• The hole always contains exactly (minus) one electron:
3
h
(
r
,
r
)
d
 X 1 2 r2  1
Dispersion Model: the Basic Idea
An electron plus its exchange hole always has zero total charge
but in general a non-zero dipole moment! This r1-dependent
dipole moment is easily obtained by integrating over r2:
 1

d X (r1 )  
rij i (r1 ) j (r1 )  r1

  (r1 ) ij

rij   r2 i (r2 ) j (r2 )d 3 r2
Note that only occupied orbitals are involved.
(This reduces, for the H atom, to the exact dipole moment of the
H atom when the electron is at r1.)
In a spherical atom, consider the following simplified “2-point”
picture:
e- (r , W)
dX
h (r - dX , W)
Notice that this
picture generates
higher multipole
moments as well
(with respect to the
nucleus as origin)
given by
Nucleus
   [r   (r  d X ) ]
and that all these moments depend only
on the magnitude of the exchange-hole
dipole moment.
This is significant because the magnitude dX(r) of the
exchange-hole dipole moment can be approximated using
local densities and the Becke-Roussel exchange-hole model
[Phys. Rev. A 39, 3761 (1989)], a 2nd-order GGA:
d X (r )  b
where b is the displacement from the reference point of the
mean position of the BR model hole [Becke and Johnson,
J. Chem. Phys. 123, 154101 (2005)].
Therefore, the entire van der Waals theory that follows has two
variants:
Orbital (XX) based, or Density-functional (BR) based
XX performs slightly better in rare-gas systems.
BR performs better in intermolecular complexes.
All our current work employs the BR (DFT) variant.
The Dispersion Interaction: Spherical Atoms
rA
Vint(rA,rB)
rB
Vint(rA,rB) = multipole moments of electron+hole at rA interacting with
multipole moments of electron+hole at rB
2nd-Order Ground-State Perturbation Theory in the
Closure (Ünsold) Approximation
If the first-order, ground-state energy correction arising from a
perturbation Vpert is zero:
E (1)  V pert  0
Then the second-order correction is approximately given by
 V pert 
2
E ( 2)  
E avg
where the expectation values are in the ground state and
Eavg is the average excitation energy.
To evaluate the expectation value <Vint2>
square the multipole-multipole interaction Vint(rA,rB):
Vint2(rA,rB) = (dipole-dipole + dipole-quadrupole +
dipole-octopole + quadrupole-quadrupole +…)2
Then integrate the squared interaction over all rA and rB.
This is a “semiclassical” calculation of <Vint2>.
The result is
where
E disp
C6 C8 C10
  6  8  10    
R
R
R
2   12  A   12  B
C6 
3
Eavg
  12  A   22  B    22  A   12  B
C8 
Eavg
4   12  A   32  B    32  A   12  B 14   22  A   22  B
C10 

3
Eavg
5
Eavg
with atomic moment integrals given by
  2     (r )[r   (r  d X )  ]2 d 3 r

What about
Eavg ?
Assume that Eavg  E A  EB and that, for each atom,
2   12 
E 
3
where α is its dipole polarizability. This easily follows from
the same “semiclassical” 2nd-order perturbation theory applied
to the polarizability of each atom.
Thus our C6, C8, C10’s depend on atomic polarizabilities and
moment integrations from Hartree-Fock (or KS) calculations!
No fitted parameters or explicitly correlated wavefunctions!
How well does it work?
On the 21 pairs of the atoms H, He, Ne, Ar, Kr, Xe,
the mean absolute percent errors are:
C6
3.4 %
C8
21.5 %
C10
21.5 %
From Free Atoms to Atoms in Molecules
• Partition a molecular system into “atoms” using Hirshfeld weight functions:
 iat (r )
wi (r ) 
  nat (r )
n
where  at is a spherical free atomic density placed at the appropriate nucleus
and the n summation is over all nuclei.
wi (r ) has value close to 1 at points near nucleus i and close to 0 elsewhere.
Also,
 w (r )  1
i
i
• Assume that an intermolecular Cm (m=6,8,10) can be written as a sum of
interatomic Cm,ij :
A
B
i
j
Cm   Cm,ij
i in A
j in B
In the previous expressions for Cm replace A and B with i and j :
C 6,ij
C8,ij 
C10,ij
2
2
2   1 i   1  j

3
(E i  E j )
  12  i   22  j    22  i   12  j
(Ei  E j )
2
2
2
2
2
2
4   1 i   3  j    3 i   1  j 14   2 i   2  j


3
(Ei  E j )
5
(Ei  E j )
with   2  i generalized to
  2  i    wi (r )  (r )[r   (r  d X )  ]2 d 3 r

and,
where
i
2   12  i
Ei 
3 i
is the effective polarizability of atom i in A. We propose that
 r 3 i
i 
 i , free
3
 r  i , free
 r 3  i   r 3 wi ( r )  ( r )d 3 r
 r 3  i , free   r 3  i , free (r )d 3 r
These are effective volume integrations. This is motivated by the well known
qualitative (if not quantitative) general relationship between polarizability
and volume. See Kannemann and Becke, JCP 136, 034109 (2012).
All radii r in the above integrals and in the integrals for   2  i
are with respect to the position of nucleus i.
Atom-Molecule and Molecule-Molecule
Dispersion Coefficients
Test set:
H2 and N2 with He, Ne, Ar, Kr, and Xe.
Cl2 with He, Ne, Ar, Kr, and Xe (except C10).
H2-H2
H2-N2
N2-N2
(only H2-H2 for C10).
Fully-numerical Hartree-Fock calculations on the monomers using
the NUMOL program (Becke and Dickson, 1989).
MAPEs with respect to dispersion coefficients from frequency
dependent MBPT polarizabilities:
12.7% for C6
16.5% for C8
11.9% for C10
(On a much more extensive test set of 178 intermolecular C6’s, the
model has a MAPE of 9.1%)
Everything, so far, has been about the asymptotic
dispersion series between atom pairs,
E disp
C6 C8 C10
  6  8  10    
R
R
R
The asymptotic series needs to be damped in order to
avoid divergences when R is small.
i.e. need information about characteristic R values inside
of which the asymptotic series is no longer valid.
The usual approach is to use empirical vdW radii.
However…
“Critical” Interatomic Separation
Since we can compute C6, C8,
and C10 non-empirically, we can
obtain non-empirical range
information.
There is a “critical” Rc,ij where
the three dispersion terms are
approximately equal:
Take Rc,ij as the average of
,
, and
The asymptotic dispersion series is obviously meaningless inside Rc,ij.
Therefore …
General Dispersion Energy Formula
… we use Rc,ij to damp the dispersion energy at small internuclear
separations as follows:
Edisp

C6,ij
C8,ij
C10,ij

  6
 8
 10
6
8

RvdW ,ij  Rij RvdW ,ij  Rij10
i  j  RvdW ,ij  Rij




RvdW ,ij  a1 Rc,ij  a2
where Rvdw,ij is an effective van der Waals separation with only two
universal fit parameters. Best-fit a1 and a2 values depend on the
exchange-correlation theories with which the above is combined.
… which brings us to the second part of the talk …
Key XDM Publications
A. D. Becke and E. R. Johnson, J. Chem. Phys. 127, 154108 (2007)
Dispersion coefficients
E. R. Johnson and A. D. Becke, J. Chem. Phys. 124, 174104 (2006)
Damping functions
A. D. Becke and E. R. Johnson, J. Chem. Phys. 123, 154101 (2005)
Transform to a DFT
Exchange-Correlation GGAs
+
XDM Dispersion
Axel D. Becke and Felix O. Kannemann
Email: [email protected]
Department of Chemistry, Dalhousie University,
Halifax, Nova Scotia, Canada
It would be nice if we could use
“standard” XC-GGAs (instead of, eg., Hartree-Fock)
plus the XDM dispersion model
to treat vdW interactions.
Let’s look at the exchange part first.
Exchange GGA Functionals B86, B86b, B88
E XB86  E XLDA
4 / 3 2
 0.0036 
1  0.0042 

E XB 86 b  E XLDA  0.0036  

E
B88
X
E
LDA
X
 4 / 3  2
1  0.006  
2 4/5

4 / 3 2
 0.0042 
1  0.0252 arcsin h 

1/ 3
where the local (spin) density part is
E XLDA
3 3 
  
2  4 
and the “reduced” (spin) density gradient is
 
 
4/3


4 / 3
  4 / 3

Exchange GGA Functionals PW86, PW91, PBE(96)
E X   e XLDA (  ) g X ( s)
g
g
PW 91
X
PW 86
X
e
where

LDA
X
3 3 1/ 3 4 / 3
( )   ( ) 
4 
(s)  1  1.296s  14s  0.2s
2
4

6 1 / 15
1  0.19645s sinh 1 (7.7956s)  [0.2743 0.1508exp(100s 2 )]s 2
( s) 
1  0.19645s sinh 1 (7.7956s)  0.004s 4
g
PBE
X
0.804
( s )  1  0.804 
1  0.21951 s 2 / 0.804
The “reduced” (spin) density gradient is
and, for spin polarized systems,
E X (  ,   ) 
s

2(3 2 )1 / 3  4 / 3
1

E X (2   )  E X (2   )
2
revPBE (Yang group)
g XPBE ( s )  1  0.804 
0.804
1  0.21951 s 2 / 0.804
g XrevPBE ( s)  1  1.245 
1.245
1  0.21951 s 2 / 1.245
Now plot the exchange-only interaction energy in Ne2
Note that Hartree-Fock (exact exchange) is repulsive!
Plot is from Kannemann and Becke, JCTC 5, 719 (2009)
Ne2
Can get anything!
from artifactual binding to massive over-repulsion,
depending on the choice of functional!
Which exchange GGA best reproduces exact Hartree-Fock?
Lacks and Gordon, PRA 47, 4681 (1993)
Kannemann and Becke, JCTC 5, 719 (2009)
Murray, Lee, and Langreth, JCTC 5, 2754 (2009)
PW86, followed by B86b
Exchange Enhancement Factor [Zhang, Pan, Yang, JCP 107, 7921 (1997)]
PW86 is a completely non-empirical exchange functional!
Its 4 parameters are fit to a theoretical exchange-hole model.
Perdew and Wang, PRB 33, 8800 (1986)
E X   e XLDA (  ) g X ( s)
g
PW 86
X

(s)  1  1.296s  14s  0.2s
2
s
4

6 1 / 15

2(3 2 )1 / 3  4 / 3
That PW86 accurately reproduces Hartree-Fock repulsion energies is remarkable.
The underlying theoretical model (truncated “GEA” hole) knows nothing about
closed-shell atomic or molecular interactions! Could be a fortuitous accident?!
Nevertheless, no parameters need to be fit to data 
Now, what about (dynamical) correlation?
Use the non-empirical “PBE” correlation functional:
Perdew, Burke, and Ernzerhof, PRL 77, 3865 (1996)
Therefore we have,
EXC  E
E
XDM
disp
PW 86
X
E
PBE
C
E
XDM
disp

C6,ij
C8,ij
C10 ,ij
1 

  6
 8
 10
6
8
10 

2 i  j  RvdW ,ij  Rij RvdW ,ij  Rij RvdW ,ij  Rij 
RvdW ,ij  a1 Rc,ij  a2
with a1 and a2 to be determined.
How to determine a1 and a2?
Fit to the binding energies of the prototypical
dispersion-bound rare-gas systems
He2 HeNe
HeAr
Ne2 NeAr
Ar2
(reference data from Tang and Toennies, JCP 118, 4976 (2003))
At the CBS limit, we find the best-fit values
a1=0.65
a2=1.68
(Kannemann and Becke, to be published)
Gaussian 09, aug-cc-pV5Z, counterpoise, ultrafine grid, BEs in microHartree
a1 = 0.65
a2 = 1.68Å
RMS%E = 4.2%
Note that in all subsequent benchmarking there is
no (re)fitting of parameters!
Our functional is, from here on, essentially
nonempirical in all its parts.
GGA + XDM Publications
F.O. Kannemann and A.D. Becke, JCTC 5, 719 (2009)
Rare-gas diatomics (numerical post-LDA)
F.O. Kannemann and A. D. Becke, JCTC 6, 1081 (2010)
Intermolecular complexes (numerical post-LDA)
A.D. Becke, A.A. Arabi and F.O.Kannemann, Can. J. Chem.
88, 1057 (2010)
Dunning aDZ and aTZ basis-set calculations (post-G09)
Tests on Standard
BioOrganic Benchmark Sets
Axel D. Becke and Felix O. Kannemann
Email: [email protected]
Department of Chemistry, Dalhousie University,
Halifax, Nova Scotia, Canada
The “S22” and “S66” vdW Benchmark Sets of Hobza et al
S22
S66
hydrogen bonding
dispersion
other noncovalent interactions
References
S22: Jurecka, Sponer, Cerny, Hobza, PCCP 8, 1985 (2006)
S66(x8): Rezac, Riley, Hobza JCTC 7, 2427 (2011)
The next slides contain Mean Absolute Deviations (MADs)
for the S22 and the S66 benchmark sets in comparison with
other popular DFT methods. Data for all methods other than
ours are from
Goerigk, Kruse & Grimme, CPC 12, 3421 (2011)
All computations employ the def2-QZVP basis set
PW86PBE-XDM
no CP (0.30)
CP (0.28)
S66
PW86PBE-XDM
no CP (0.27)
CP (0.23)
Aside
Can we do better by combining XDM with
hybrid functionals?
(rather than the pure GGA, PW86+PBE)
Burns, Vazquez-Mayagoitia, Sumpter, Sherrill, JCP 134, 084107 (2011)
S22 MAD (kcal/mol) for B3LYP-XDM and other DFT methods
However
• The pure GGA, PW86+PBE, is completely nonempirical
(B3LYP, with 3 fitted parameters, is not)
• Density-fit basis sets can speed up pure GGA calculations,
with no loss of accuracy, by an order of magnitude!
Tests on the S66x8 Set
(Nonequilibrium Geometries)
Axel D. Becke and Felix O. Kannemann
Email: [email protected]
Department of Chemistry, Dalhousie University,
Halifax, Nova Scotia, Canada
S66x8
A benchmark set of 66 vdW complexes of bio-organic interest,
at 8 intermonomer separations:
0.90, 0.95, 1.00, 1.05, 1.10, 1.25, 1.50, 2.00
(relative to the equilibrium intermonomer separation)
Rezac, Riley, Hobza JCTC 7, 2427 (2011)
Important because complex systems and materials may contain
many vdW interactions between groups at nonequilibrium
(especially stretched) geometries!
Mean Percent Errors (MPEs)
versus geometry
10
0
MPE [%]
-10
-20
PW86PBE+XDM
-30
M06-2X
M05-2X
-40
M06L
-50
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
25
20
15
MPE [%]
10
5
0
PW86PBE+XDM
-5
M06-2X-D3
-10
M05-2X-D3
-15
M06L-D3
-20
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
25
20
15
MPE [%]
10
5
0
PW86PBE+XDM
-5
PBE-D3
-10
BLYP-D3
-15
B97-D3
-20
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
10
5
MPE [%]
0
-5
PW86PBE+XDM
LC-ωPBE-D3
-10
B3LYP-D3
-15
B2-PLYP-D3
-20
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
Mean Absolute Percent Errors (MAPEs)
versus geometry
80
PW86PBE+XDM
70
M06-2X
60
M05-2X
MAPE [%]
50
M06L
40
30
20
10
0
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
35
PW86PBE+XDM
30
M06-2X-D3
MAPE [%]
25
M05-2X-D3
M06L-D3
20
15
10
5
0
0.90
1.00
1.10
1.25
distance multiplier
1.50
2.00
35
PW86PBE+XDM
30
PBE-D3
25
BLYP-D3
MAPE [%]
B97-D3
20
15
10
5
0
0.90
1.00
1.10
1.25
distance multiplier
1.50
2.00
35
PW86PBE+XDM
30
LC-ωPBE-D3
MAPE [%]
25
B3LYP-D3
B2-PLYP-D3
20
15
10
5
0
0.90
1.00
1.10
1.25
distance multiplier 1.50
2.00
pdf file of all S66x8 curves
(Notice the importance of the dispersion term!)
Interaction types in S66x8:
• H-bonding
23
• Dispersion
23
• Mixed (“other”)
20
Dispersion is further divided into:
pi-pi (10), aliphatic-aliphatic (5), and pi-aliphatic (8)
How good is PW86+PBE+XDM for ordinary thermochemistry?
Consider the functional
XDM
EXC  EXGGA  ECPBE  Edisp
with a variety of standard exchange GGAs in the first term.
On the “G3/99” benchmark set of 222 atomization energies of organic/inorganic
molecules (Curtiss, Raghavachari, Pople) we obtain the following error statistics,
in kcal/mol:
For standard hybrid functionals (eg., B3LYP) the MAE is of order 5-6 kcal/mol.
The very best DFTs have MAE as small as 2-3 kcal/mol, but with fitted params!
Availability of XDM code
• Has been implemented (B3LYP-XDM) in Q-Chem by
Kong, Gan, Proynov, Freindorf, and Furlani, PRA 79, 042510 (2009)
• A “post-Gaussian09” code will be available from us by
the end of 2012. Uses G09 to perform the PW86+PBE part, then
adds XDM perturbatively. Can do Berny geometry optimizations using
the EXTERNAL keyword! Very fast with density-fit basis sets!
Many thanks to:
Natural Sciences and Engineering Research Council of Canada
the Killam Trust of Dalhousie University (Killam Chair)
ACEnet (the Atlantic Computational Excellence Network)

similar documents