Transcript Electronic Structure with DFT: GGA and beyo

ELECTRONIC STRUCTURE WITH DFT:
GGA AND BEYOND
Jorge Kohanoff
Queen’s University Belfast
United Kingdom
[email protected]
DENSITY FUNCTIONAL THEORY (DFT):
KOHN-SHAM EQUATIONS

Kohn-Sham equations:
 2 2
 KS
KS



V
[

](
r
)
KS

 n (r)   n n (r)
 2m

N
 KS (r )   (r )    (r )
n 1
VKS [  ]  Vext (r )  

KS
n
 (r' )
| r  r'|
2
dr '   XC [  ]
This PDE must be solved self-consistently, as the KS
potential depends on the density, which is constructed
with the solutions of the KS equations.
EXCHANGE AND CORRELATION IN DFT:
THE LOCAL DENSITY APPROXIMATION (LDA)

The inhomogeneous electron gas is considered as locally
homogeneous:
LDA
E XC [  ]    (r )  XC [  ](r ) dr

LDA
XC
[ ]  
~XC (r, r ' )
| r  r'|
dr '
XC energy density of the HEG
  (r ) 
LDA
HEG
g~XC
(r, r ' )  g~XC
(| r  r ' |,  (r )) 

  (r ' ) 


LDA XC hole centred at r, interacts with the electron also at r.
The exact XC hole is centred at r’
This is partially compensated by multiplying the pair correlation
function with the density ratio (r)/(r’)
EXCHANGE AND CORRELATION IN DFT:
THE LOCAL DENSITY APPROXIMATION (LDA)

Location of the XC hole (Jones and Gunnarsson, 1982)
LDA-LSDA: TRENDS AND LIMITATIONS

Favors more homogeneous electron densities

Overbinds molecules and solids (Hartree-Fock underbinds)

Geometries, bond lengths and angles, vibrational frequencies
reproduced within 2-3%

Dielectric constants overestimated by about 10%

Bond lengths too short for weakly bound systems (H-bonds, VDW)

Correct chemical trends, e.g. ionization energies

Atoms (core electrons) poorly described (HF is much better)

XC potential decays exponentially into vacuum regions. It should
decay as –e2/r. Hence, it is poor for dissociation and ionization

Poor for metallic surfaces and physisorption

Very poor for negatively charged ions (self-interaction error)

Poor for weakly bound systems: H-bonds (), VDW (non-local)

Band gap in semiconductors too small (~40%)

Poor for strong on-site correlations (d and f systems, oxides, UO2)
BEYOND THE LDA

Inhomogeneities in the density

Self-interaction cancellation

Non-locality in exchange and correlation

Strong local correlations

Gradient expansions

Weighted density approximation

Exact exchange in DFT (OEP local vs HF non-local)

DFT-HF hybrids

Self-interaction correction

Van der Waals and RPA functionals

LSDA+U

Multi-reference Kohn-Sham

GW approximation (Many-body)
GRADIENT EXPANSIONS:
GENERALIZED GRADIENT APPROXIMATION

EXC expanded in gradients of the density
LSDA
EGGA[  ]    (r )  XC
[  ,  ](r ) FXC [  ,  , s](r )dr
where  is the spin polarization
(r+dr)
r+dr
s=||/2kF is the density gradient
And FXC is the enhancement factor


First-order term is fine, but higher-order terms diverge. Only by
some re-summation to ∞-order the expansion converges.
GGA: FXC is designed to fulfil a number of exactly known
properties, e.g. Perdew-Burke-Ernzerhof (PBE)
1.
Exchange: uniform scaling, LSDA limit, spin-scaling
relationship, LSDA linear response, Lieb-Oxford bound
2.
Correlation: second-order expansion, hole sum rule, vanishes
for rapidly varying densities, cancels singularity at high densities
PROPERTIES OF THE GGA

Improves atomization and surface energies

Favors density inhomogeneities

Increases lattice parameters of metals

Favors non-spherical distortions

Improves bond lengths

Improves energies and geometries of H-bonded systems

There is error cancellation between X and C at short range

XC potential still decays exponentially into vacuum regions

Some improvement in band gaps in semiconductors

What was correct in LDA is worsened in GGA

Still incorrect dissociation limit. Fractionally charged fragments

Inter-configurational errors in IP and EA

Error cancellation between X and C is not complete at long-range.
X hole is more long-ranged than XC hole
HYBRID FUNCTIONALS

Combine GGA local exchange with Hartree-Fock non-local
exchange:
EHYBRID[  ]   E XGGA[  ]  (1   ) E XHF [  ]  ECGGA[  ]

Parameter  fitted to experimental data for molecules (~0.75), or
determined from known properties.

PBE0, B3LYP, HSE06

Properties:
1.
Quite accurate in many respects, e.g. energies and geometries
2.
Improve on the self-interaction error, but not fully SI-free
3.
Improve on band gaps
4.
Improve on electron affinities
5.
Better quality than MP2
6.
Fitted hybrids unsatisfactory from the theoretical point of view
SELF-INTERACTION CORRECTION (SIC)

Self-interaction can be removed at the level of classical
electrostatics:
1  (r )  (r ' )
EH  
drdr '
2
| r  r'|
 (r )  n (r ' )
1 N
Perdew-Zunger 1982
ESIC  EH    n
drdr '
Mauri, Sprik, Suraud
2 n 1
| r  r'|
(n)
VSIC
(r )  VH (r )  




 n (r ' )
| r  r'|
dr '
Potential is state-dependent. Hence it is not an eigenvalue
problem anymore, but a system of coupled PDEs
Orthogonality of SIC orbitals not guaranteed, but it can be
imposed (Suraud)
Similar to HF, but the Slater determinant of SIC orbitals is
not invariant against orbital transformations
The result depends on the choice of orbitals (localization)
DYNAMICAL CORRELATION: VDW


Van der Waals (dispersion) interactions: are a dynamical
non-local correlation effect
Dipole-induced dipole interaction due to quantum density
fluctuations in spatially separated fragments
1(r,t)

2(r,t)
Functional (Dion et al 2004):
1(t)
EVDW    (r ) (r, r ' )  (r ' ) drdr '
E1(t).2(t)
 = VDW kernel fully non-local.
Depends on (r) and (r’)

Expensive double integral

Efficient implementations (Roman-Perez and Soler 2009)

Good approximations based on dynamical response theory

Beyond VDW: Random Phase Approximation (Furche)
DYNAMICAL CORRELATION: VDW AND BEYOND

Empirical approaches:
P
EVDW  
IJ




C6IJ ( RIJ )
f IJ ( RIJ )
RIJ6
With f a function that removes smoothly the singularity at
R=0, and interferes very little with GGA (local) correlation.
Grimme (2006): C6 parameters from atomic calculations.
Extensive parameterization: DFT-D.
Tkatchenko and Scheffler (2009): C6 parameters dependent
on the density.
Random Phase Approximation (RPA): captures VDW and
beyond. Can be safely combined wit exact exchange (SIC).
Infinite order perturbation (like Coupled clusters in QC).
Furche (2008); Paier et al (2010); Hesselmann and Görling (2011)
STRONG STATIC CORRELATION: LSDA+U



Strong onsite Coulomb correlations are ot captured by LDA/GGA
These are important for localized (d and f) electronic bands, where
many electrons share the same spatial region: self-interaction
problem
Semi-empirical solution: separate occupied and empty state by an
additional energy U as in Hubbard’s model:
1
1
E LSDAU  ELSDAU  U N ( N  1)  U  f i f j
2
2 i j

fi=orbital occupations
This induces a splitting in the KS eigenvalues:
E
1

 i  LSDAU   iLSDA  U   f 
f i
2

 iocc   iLSDA  U / 2
 iempty   iLSDA  U / 2
SUMMARY OF DFT APPROXIMATIONS
ELECTRONIC STRUCTURE OF
UO2
Using the quantum-espresso package
(http://www.quantum-espresso.org/)
Pseudopotentials
Plane wave basis set
•
•
PROPERTIES










fluorite structure
fcc, 3 atoms un unit cell
Lattice constant = 10.26 Bohr
Electronic insulator. Eg=2.1 eV
Electronic configuration of U: [Rn]7s26d15f3
U4+: f2
5f-band partially occupied (2/7)
UO2: splitted by crystal field:
t1u(3)+t2u(3)+ag(1)
Still partially occupied (2/3)
Jahn-Teller distortion opens gap.
PSEUDOPOTENTIAL
CONVERGENCE WITH ENERGY CUTOFF
ENERGY-VOLUME CURVE
GGA(PBE) DENSITY OF STATES
GGA+U DENSITY OF STATES
GGA+U DENSITY OF STATES: DISTORTED
VAN DER WAALS FOR IMIDAZOLIUM SALTS
 An important family of Room Temperature Ionic
Liquids (Green solvents)
 Competing electrostatic vs dispersion interactions
 Large systems studied with DFT, within LDA or GGA
[Del Popolo, Lynden-Bell and Kohanoff, JPCB 109, 5895 (2005)]
 Force fields fitted to DFT-GGA calculations
[Youngs, Del Popolo and Kohanoff, JPCB 110, 5697 (2006)]
 Electrostatics well described in DFT (LDA or GGA)
 Dispersion (van der Waals) interactions are absent in
both, LDA and GGA
DFT IMIDAZOLIUM SALTS
RESULTS FOR SIMPLE DIMERS
M. Dion et al, PRL 92, 246401 (2004)
Ar2 and Kr2
Bond lengths:
(C6H6)2
5-10% too long
Binding energies: 50-100% too large
RESULTS FOR SOLIDS
Polyethylene
Silicon
• Reasonable results for molecular systems
• Keeps GGA accuracy for covalent systems
 General purpose functional
SOLID FCC ARGON (E. ARTACHO)
Some overbinding, and lattice constant still 5% too large
… but much better than PBE
(massive underbinding and lattice constant 14% too large)
THE DOUBLE INTEGRAL PROBLEM

(q1,q2,r12) decays as r12-6

Ecnl = (1/2)   d3r1 d3r2 (r1) (r2) (q1,q2,r12) can be truncated for
r12 > rc ~ 15Å

In principle O(N) calculation for systems larger than 2rc ~ 30Å

But... with x ~ 0.15Å (Ec=120Ry) there are ~(2106)2 = 41012
integration points

Consequently, direct evaluation of vdW functional is much more
expensive than LDA/GGA
FACTORING (Q1,Q2,R12)
G. Román-Pėrez and J. M. Soler, Phys. Rev. Lett. 2009
 ( q1 , q2 , r12 )   p ( q1 ) p ( q2 )   ( r12 )
 ,
Expand  in a basis set of functions p(q)
1
E   dr1dr2  (r1 )  (r2 )  ( q1 , q2 , r12 )
2
1
   dr1dr2  (r1 )   (r2 )  ( r12 )
2  ,
nl
c
1
   dk  ( k )   ( k )  ( k )
2  ,
 (r )   (r ) p q (r ),  (r ) 
FT
INTERPOLATION AS AN EXPANSION
Basis functions p(q) interpolate between grid points
• Lagrange polynomials: grid given by zeros of orthogonal polynomials
• Cubic splines: grid points defined on a logarithmic mesh
f ( x)   f i pi ( x)
f i  f ( xi )
i
 in ( x  x j )
if xi  n  x  xi  n
 
pi ( x)   j j ii n ( xi  x j )

0 otherwise
20 grid points are sufficient
Smoothening of  required at small q
O(NLOGN) ALGORITHM
do, for each grid point i
find i and i
Implemented into SIESTA, but
not SIESTA-specific:
find qi=q(i ,i )
find i = i p(qi )

i
end do
Input:
Fourier-transform i  k 
Output: Exc , vixc on the grid
do, for each reciprocal vector k
find uk =  (k) k
No need for supercells in solids

No cutoff radius of interaction
end do
Inverse-Fourier-transform uk  ui

do, for each grid point i
find i , i , and qi
find i , i /i , and i / i
find vi
end do
on a regular grid

ALGORITHM EFFICIENCY
Atoms
CPU time in
GGA-XC
CPU time in
vdW-XC
vdW/GGA
overhead
2
0.75 s (44%)
7.5 s (89%)
400%
MMX polymer
124
1.9 s (2%)
10.6 s (16%)
17%
DWCN
168
11.9 s (0.6%)
109 s (5.2%)
4%
System
Ar2
Message
If you can simulate a system with LDA/GGA,
you can also simulate it with vdW-DFT
IMIDAZOLIUM CRYSTALS: VOLUMES
LDA
PBE
VDW
EXP
[emim][PF6]
877.4
1088.3
1059.1
1023.9
[bmim][PF6]
513.4
636.6
620.0
605.0
[ddmim][PF6]
1710.9
2258.0
2095.5
2000.6
[mmim][Cl]
607.9
750.5
728.6
687.6
[bmim][Cl]-o
843.8
1039.9
1019.8
961.1
[bmim][Cl]-m
836.3
1051.9
1024.6
966.7
-13.5%
+8.5%
+5%
Error
IMIDAZOLIUM CRYSTALS: HEXAFLUOROPHOSPATES
emimpf6
bmimpf6
ddmimpf6
IMIDAZOLIUM CRYSTALS: HEXAFLUOROPHOSPATES
[bmim][PF6]
C5-C5
C1”-C1”
C4’-C4’
C2-C4’
C4’-C1”
P-P
C2-P
C1”-P
C4’-P
Error
LDA
4.32
4.24
3.73
4.61
3.64
5.05
3.70
3.88
4.75
7.0%
PBE
4.64
4.85
3.99
5.07
3.97
5.69
4.16
4.21
4.95
2.6%
VDWE
4.60
4.55
4.05
4.80
4.07
5.64
4.21
4.20
4.91
2.6%
VDW
4.54
4.60
4.04
4.80
3.95
5.62
3.90
4.11
4.87
2.0%
EXP
4.56
4.52
4.05
5.08
3.86
5.54
3.95
4.29
4.94
[ddmim][PF6]
C4-C4
C1”C1”
C12C12’
C4’-C4’
C2-C4’
C4’C1”
P-P
C2-P
C1”-P
C4’-P
C12’-P
Error
LDA
5.18
4.28
4.83
4.52
4.23
6.00
6.66
3.55
3.95
4.27
4.35
6.0%
PBE
5.62
4.24
5.44
5.12
5.10
6.56
6.89
4.00
4.75
4.96
4.78
7.0%
VDW
5.47
4.29
5.23
4.90
4.79
6.44
6.80
3.82
4.22
4.68
4.65
3.3%
EXP
5.47
4.36
5.01
4.85
4.61
6.25
6.45
4.23
4.01
4.62
4.63
IMIDAZOLIUM CRYSTALS: CHLORIDES
[mmim][Cl]
C5-C5
C1’-C1’
C2-C1’
C2-Cl
C1’-Cl
LDA
3.07
3.54
3.37
3.34
3.81
PBE
3.94
3.86
3.63
3.48
3.74
VDW
3.34
3.82
3.56
3.56
3.77
[bmim][Cl]-m
C5-C5
C1”-C1” C4’-C4’
C2-C4’
C4’-C1”
C2-Cl
C1”-Cl
C4’-Cl
LDA
4.87
4.50
5.17
3.25
3.83
3.28
3.51
3.70
PBE
5.26
5.14
5.72
3.63
3.81
3.38
3.86
3.89
VDW
5.21
4.97
5.42
3.53
3.84
3.46
3.84
3.94
EXP
5.19
4.96
5.33
3.48
3.62
3.39
3.82
3.80
[BMIM][CL]: POLYMORPHISM
[bmim][Cl]-Monoclinic
[bmim][Cl]-Orthorhombic
Energy difference per neutral ion pair
E(M-O) [eV]
LDA
PBE
VDW
FF(CLaP)
-0.05
-0.06
-0.01
0.08
IMIDAZOLIUM CLUSTERS: TRIFLATES
[bmim][Tf]
[bmim][Tf]2
Geometry
Error
LDA
PBE
VDW
FF(CLaP)
7.2%
2.3%
1%
--
Dimer association energy
E [eV]
LDA
PBE
VDW
FF(CLaP)
-1.411
-1.053
-1.453
-1.457
CONCLUSIONS
1.
The non-empirical van der Waals functional of Dion et al.
(DRSLL) improves significantly the description of the geometry
of imidazolium salts.
2.
Volumes are improved respect to PBE, but still overestimated by
5%.
3.
Energetics is also improved. It is similar to that of empirical
force fields such as CLaP.
4.
The cost of calculating the van der Waals correlation correction
is 10 times that of PBE. However, in a self-consistent calculation
for 100 atoms the overhead is only 20%.
THE THEORETICAL LANDSCAPE
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
QM/MM


Treat relevant part of the system quantum-mechanically, and the
rest classically.
The problem is how to match the
two regions. Easy for non-bonded
interactions, more difficult for
chemical bonds

One can also treat part of the
system as a polarizable
continuum, or reaction field (RF)
QUANTUM MECHANICS IN A LOCAL BASIS
TIGHT BINDING
TIGHT BINDING
TIGHT BINDING MODELS FOR WATER
A. T. Paxton and J. Kohanoff, J. Chem. Phys. 134, 044130 (2011)
Ground-up philosophy
Water molecule
1. Minimal basis. On-site energies to reproduce band structure
a) 1s orbital for H: Hs
b) 2s & 2p orbitals for O: Os, Op
2. O-H hopping integrals:
a) Values at equilibrium length to reproduce HOMO-LUMO gap: tss, tsp
b) GSP functional form. Cut-off between first and second neighbours
3. Charge transfer: Hubbard terms fitted to reproduce dipole moment: UO, UH
4. O-H pair potential: GSP form, fitted to reproduce bond length and
symmetric stretching force constant
5. Crystal field parameter spp selected to reproduce polarizability
TIGHT BINDING MODELS FOR WATER
Ground-up philosophy
Water dimer
6. O-O hopping integrals: tss, tsp , tpp, tpp
GSP form. Cut-off after first neighbours
7. O-O pair potential
Various forms (GSP, quadratic) to
reproduce binding energy curve
8. Fitting procedure
a) By hand (intuitive)
b) Genetic algorithm
This is the end of the fitting
All the rest are predictions
ICE-XI
• DFT ice sinks in water!
• Polarizable model marginal
• Point charge model is fine
TIGHT-BINDING LIQUID WATER
A. T. PAXTON
AND
J. KOHANOFF, J. CHEM. PHYS. 134, 044130 (2011)
TIGHT BINDING MODEL FOR TIO2
Band structure: A. Y. Lozovoi, A. T. Paxton and J. Kohanoff
DFT
TB
Rutile
Anatase
O on-site energies Os, Op and O-O hopping integrals: tss, tsp , tpp, tpp
WATER/TIO2 INTERFACES (SASHA LOZOVOI)
Single water adsorption: dissociation
WATER/TIO2 INTERFACES (SASHA LOZOVOI)
A water layer on TiO2
WATER/TIO2 INTERFACES (SASHA LOZOVOI)
Bulk water on TiO2
• 4 TiO2 layers
(480 atoms)
• 184 water
molecules
(552 atoms)
• 1032 atoms
• TBMD 1 ps
(one week)
WATER/TIO2 INTERFACES (SASHA LOZOVOI)
z-density profile of bulk water between TiO2 surfaces
PROTON TRANSFER IN WATER
Chemical bonds broken and formed
ORGANIC MOLECULES (TERENCE SHEPPARD)
Benzylacetone in water
BENZYLACETONE AT WATER/TIO2
INTERFACE
Sasha Lozovoi
and
Terence Sheppard
SUMMARY

In developing a force field there are 3 main aspects to
consider. In order of relevance:
Designing the model
 Choosing the properties to be reproduced
 Choosing a methodology to fit those properties





Fitting procedures cannot work out their magic if the
model is not good enough
To become model-independent we need to introduce
electrons explicitly: ab initio and DFT methods
These are expensive. Simplifications like semi-empirical,
tight-binding and QM/MM methods bridge the gap
Further simplifications lead to extremely efficient bondorder potentials (BOP)