Simplified Methods for Electronic Structure Calculations

```SIMPLIFIED METHODS FOR ELECTRONIC
STRUCTURE CALCULATIONS
Jorge Kohanoff
Queen’s University Belfast
United Kingdom
[email protected]
CCMS Summer Institute
LLNL 23 July 2012
QUEEN’S UNIVERSITY BELFAST
NORTHERN IRELAND, UK
THE PROBLEM:
NUCLEI AND ELECTRONS INTERACTING VIA

COULOMB FORCES
Hamiltonian of the universe
H  Tn  Vnn  Te  Vee  Vne
2
Tn  
2
P
1 2
I

M
I 1
I
ZI ZJ
e2 P P
Vnn  
2 I 1 J  I | R I  R J |
2 N 2
Te  
i

2m i 1

e2 N N
1
Vee  
2 i 1 j i | ri  r j |
e2 N P
ZI
Vne   
2 i 1 I 1 | ri  R I |
Schrödinger equation:
 (r, R; t )
i
 H (r, R; t )
t
r={ri } i=1,...,N
R={RI } I=1,...,P
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
CHEMISTRY VS PHYSICS OR …
WHEN PHYSICS MEETS CHEMISTRY?
Quantum Chemistry
Postulate an approximate many-electron wave function and
find it by minimizing the expectation value of the energy for
this wave function (Rayleigh-Ritz variational principle).
Wave function (ab initio) methods
Condensed Matter Physics
Express the energy of the system as an approximate functional
of the electronic density. This functional is minimized with
respect to the density, subject to a normalization constraint.
Density functional (first principles) methods
SIMPLIFIED ELECTRONIC STRUCTURE METHODS:
GENERAL PRINCIPLES (POPLE AND BEVERIDGE)


1.
2.
3.
4.
5.
Treating electrons explicitly is expensive.
No electrons at all is generally insufficient.
The method should be sufficiently simple to allow for the
study of large systems
Approximations shouldn’t modify significantly the
physical and chemical properties of the system, e.g.
structure, energy levels, vibrations, etc.
The approximate wave function should be as unbiased as
possible
Method should take into account active (valence) electrons
It should be sufficiently general to allow for systematic
improvements. Ideally, it should be possible to derive it
from ab initio through controlled approximations.
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
EMPIRICAL (CLASSICAL) FORCE FIELDS





Many-body expansion (2-, 3-, 4-body): general purpose FF, have been
tuned mostly for biological systems (AMBER, CHARMM, GROMACS,
LAMMPS). Can be used for other systems, but generally require validation
and re-calibration (e.g. ionic liquids).
Embedded atom models (many-body): useful for systems with
delocalized electrons like metals (Al, Au, Ni). The electronic information is
included in the embedding many-body potential. BOPs.
Polarizable and Shell models: Developed to introduce oxygen
(electronic) polarization in oxides (e.g. ferroelectric perovskites). They
introduce either point dipoles or an electronic shell variable that can
displace with respect to the core, thus generating a dipole. Dipoles/shells
Reactive force fields: Useful to describe systems where there are specific
chemical changes (reactions). Extremely efficient, but generally ad hoc for
a particular system.
Combining systems is not obvious, e.g. metals with water or oxides
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
ORBITAL-FREE OR DENSITY-ONLY
FUNCTIONAL METHODS

Thomas-Fermi (1927): Approximation for the kinetic energy from
the homogeneous electron gas.
1  (r )  (r ' )
ETF [  ]  TTF [  ]   Vext (r )  (r ) dr  
dr dr '
2
| r  r'|
3 2

TTF   tTF [  ]  (r) dr   
(3 ) 2 / 3  2 / 3   (r) dr  CK   5 / 3 (r) dr
 5 2m


Minimizing the functional with respect to (r), under the constraint
that (r) integrates to N, we obtain an integral equation for (r):
5
 (r' )
CK  2 / 3 (r)  Vext (r)  
dr '  
3
| r  r'|
Functional minimization
 = electronic
chemical potential
E X [  ]  C X   4 / 3 (r ) dr
Dirac exchange
EC [  ]    A 4 / 3 /( B   1/ 3 ) dr
Wigner correlation
TVW [  ] 


1
|   | 2 /  dr

8
Von Weiszäcker
ORBITAL-FREE OR DENSITY-ONLY
FUNCTIONAL METHODS

Based on Thomas-Fermi idea, these methods require solving a
single integro-differential equation.
T [  ]
 (r' )
 Vext (r )  
dr'   XC (r )  
 (r)
| r  r'|



Need kinetic energy as explicit functional of the density. It’s
an unsolved problem. Otherwise, we would all be using OF-DFT.
Inspired in the homogeneous electron gas, they tend to work quite
well for simple metals: Na, K (alkaline).
For “less-simple” (sp) metals, e.g. Al, the TF-vW kinetic functional
is not good enough. Non-local versions that impose linear
response improves a lot (Smargiassi and Madden 1994, Wang et
al, 1999).
ORBITAL-FREE OR DENSITY-ONLY
FUNCTIONAL METHODS

TTF does not reproduce the shell structure for atoms. This can be
introduced by an ad hoc modulating function (r) in the kinetic
functional (King and Handy, 2001)
2




1
5/ 3

 dr
TS [  ]  CK   (r) (r) 

8  




Also non-local linear response functionals help with the shell
structure.
OF-DFT admits only local external potentials.
Core-electrons can be removed by introducing pseudopotentials, but
these must be local, which is possible but not easy (Carter).
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
QM/MM

Treat relevant part of the system quantum-mechanically (QM), and the
rest classically (MM).
Chemical reactions in a solvent
 Active site of enzymes


Two main cases



No chemical bonds between QM and MM
QM and MM regions involve chemical bonds
Q,
R,
Part of the system can be treated as a
polarizable continuum, or reaction field
(RF)

Energy:
EQM MM [](R, Q, )  EQM [](R)  EMM (Q, )  Eint [](R, Q, )
QM/MM

QM-MM interaction energy:
S
 (r)
K 1
r  τK
Eint [  ](R, Q, )   QK 

S
P
dr  
K 1 I 1
S
P
QK Z I
 VSR  R I  τ K
R I  τ K K 1 I 1

Short-range replusion
Avoids Q<0 collapsing
Electrostatics: e-n + n-n
onto QM nuclei
Hamiltonian: QM wfn feels the electrostatics from MM region
S
QK
Hˆ QM MM [](R, Q, )  Hˆ QM [](R)  
K 1 r  τ K

Forces:


QM nuclei feel electrostatics and SR repulsion from MM particles
MM particles feel electrostatics from nuclei and electrons + SR repulsion
QM/MM

Cutting chemical bonds: electrons in the QM region must
“believe” that on the MM side there are orbitals to make bonds.
MM

C
N
C
O
QM
The safest strategy is to cut homoatomic bonds, e.g. C-C, and replace the
MM part with a pseudopotential (Laio) or frozen orbitals (Monard)

MM force fields should be compatible with the QM model

Reaction field: mimics the dielectric response of the medium.

Onsager: QM dipole moment polarizes the medium, which in
return produces an electric field back in QM:
2  1 p


Over-polarization (catastrophe)
Improvements: PCM, Oniom, etc.
E RF 
2  1 Rc3
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
COMPUTATIONAL METHODS
100
1000
10000
100000
High
Low
Size (number of atoms)
1000000
SOLVING KS/HF EQUATIONS

Kohn-Sham equations:
 2 2
 KS
 (r' )
KS



V
(
r
)

d
r
'


[

](
r
)

(
r
)



ext
XC
n n (r )

 n

| r  r'|
 2m

N
 (r )   f n  (r )
n 1
KS
n
2
KS *
KS

(
r
)

(r ) dr   nm
n
m


Density constructed with KS orbitals
fn = occupation numbers
KS/HF orbitals must be orthogonal
to be solved self-consistently (Hartree-Fock similar)
--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
Nuclear-electron interaction
Basis sets
All-electron
Basis functions (), Size (M)
Pseudopotentials
M
 (r )   cn  (r)
KS
n
 1
BASIS SETS

Basis set representation of one-electron orbitals:
M
 (r )   cn  (r)
KS
n

(r) = Basis functions
 1
(r) can depend on energy (eigenvalue). In that case one has
to solve a complex non-linear secular equation (e.g. KKR)

If they are energy-independent, the KS equations become a
Generalized Eigenvalue Problem:
Expansion coefficients
H    (r) H  (r) dr
*
S   * (r)  (r) dr

H   S c   0


M
1
Hamiltonian
n
n
Overlap
The same can be done for Hartree-Fock (Roothaan-Hall eq.)
LOCAL ORBITALS

Atom-centered basis functions:
I (r)   (r  R I )

RI = atomic positions
(can use other positions)
Types:
Atomic orbitals (AO)
 Slater-type orbitals (STO): exp(- r)
 Gaussian-type orbitals (GTO): exp(- r2)
 Numerical


 (r  R I )  (r  R J )
Generally these are non-orthogonal
M

Löwdin orthogonalization:
1/ 2
 ' (r)   S
 (r)
Standard Eigenvalue Problem:
H '  I c'   0


 1
M

1
n
n
THE HAMILTONIAN IN A LOCAL BASIS

One-electron orbitals:
M
 1 I 1

M = number of basis functions

P = number of atoms

Eigenvalue equation:
 H
M
P
 1 J 1

Overlap matrix:

Size of basis set (M)
P
 (r)   cnI  (r  R I )
KS
n
IJ


IJ
  n S
cnJ  0
IJ
S
  * (r  R I )  (r  R J ) dr
SZ: One function per occupied orbital (in the neutral atom)
 DZ: Two functions per occupied orbital
 DZP: Two functions per occupied orbital and one per empty orbital


Example: O(1s22s22p4). SZ: 1s+3p / DZ: 2s+6p / DZP: 2s+6p+5d
THE DENSITY MATRIX

IJ
(r)  * (r  RI ) (r  R J )
Differential overlaps: O
M
P
 O (r)



Electronic density:
 (r) 

Density matrix:
   f n cnJ cnI*
IJ
IJ
1IJ 1
N
IJ
n 1
Two-point density (r,r’) expressed in a local basis.

Bond orders: Strength of the bond between two atoms IJ
 
IJ
M
 
IJ
 1
s
ss
s
THE HAMILTONIAN MATRIX


Fock matrix:
KS matrix:
M
H
IJ
1,
IJ
KS ,

P
   LK     
F  h
IJ
 1 KL 1
M
h
IJ
1,

P
IJ
   LK     XC
,
 1 KL 1
2
P

2
ˆ
h1  
   uK (r)
2m
K 1

One-electron operator:

Kinetic integrals:
T

Nuclear attraction:
U    * (r  R I )uK (r  R K ) (r  R J ) dr

I=J=K: One-center integrals

IJ=K: Two-center integrals (and permutations)

IJK: Three-center integrals
IJ
2
*



(r  R I )   (r  R J ) dr

2m 
P
IJ
K 1
TOTAL ENERGY AND DENSITY MATRIX

Two-electron integrals (Coulomb and Exchange):
   
* (r  R I )* (r  R K ) (r  R J ) (r  R L )
| r  r' |
dr dr'
One-, two-, three- and four-center integrals

HF energy:
EHF




KS energy:
EKS


1 M P JI IJ
1
     h1,  FIJ  Vnn  Tr ˆ hˆ1  Fˆ  Vnn
2  1 IJ 1
2




1 M P JI IJ
1
IJ
     h1,  H KS ,  Vnn  Tr ˆ hˆ1  Hˆ KS  Vnn
2  1 IJ 1
2
Also forces on atoms and stresses can be expressed in terms
of the density matrix. If non-orthogonal, also overlap.
BASIC APPROXIMATIONS
1.
Eliminate core electrons, as they do not participate in
bonding. Reduce size of H matrix and number of
necessary eigenstates.
2.
Use minimal basis set (SZ). Reduces size of H matrix.
3.
Simplify the calculation of matrix elements: Makes
computation of H faster.
1.
2.
3.
Eliminate expensive integrals (3- and 4-center)
Parameterize the remaining integrals
Cut off the range of the interactions.
QUANTUM
CHEMISTRY
CONDENSED MATTER
PHYSICS
SEMIEMPIRICAL
TIGHT-BINDING
THE TIGHT-BINDING APPROACH
 Non-interacting electrons approximation: electronelectron interaction is subsumed into the electronnuclear interaction:
2
P

2
Hˆ TB  
  U K (r)
2m
K 1
 Hamiltonian matrix elements:
H
IJ
TB ,
 2 2 P

   (r  R I ) 
  U K (r  R K ) (r  R J ) dr
K 1
 2m

*
 Onsite terms: I=J
II
I
HTB


,

atomic energies
 Hoppings: IJ
IJ
IJ
HTB
,  t
interatomic bonding
TIGHT-BINDING: INTERPRETATION
 (r  R I )
 Two-center hopping integrals:
 (r  R J )
U I (r  R I )
The electron feels the attraction of both nuclei
I
 Three-center hopping integrals:
U J (r  R J )
 (r  R I )
J
 (r  R J )
The hopping between I and J is influenced by
the presence of another atom K
I
U K (r  R K )
K
 Onsite terms:
Also influenced by the presence of neighboring atoms
P
H
II
TB ,
       * (r  R I ) U K (r  R K )  (r  R I ) dr
I ( 0)
K I
J
TIGHT-BINDING: TWO CENTER APPROXIMATION
 Two-center approximation: Neglect environmental
dependence of the TB parameters by ignoring:
 3-center hopping integrals
H
IJ
TB ,
 2 2

   (r  R I ) 
  U I (r  R I )  U J (r  R J ) (r  R J ) dr
 2m

*
 2-center integrals in onsite terms
H
II
TB ,
 
I ( 0)
 2 2

   (r  R I ) 
  U I (r  R I ) (r  R I )dr
 2m

*
 More sophisticated TB models re-introduce this
environmental dependence at a parametric level.
TIGHT-BINDING: PARAMETERIZATION
 Two-center integrals are not calculated explicitly. It
wouldn’t make sense, as neglected 3- and 4-center
integrals have been incorporated into 2-center ones.
 Instead, onsite energies and hoppings are treated as
fitting parameters.
 Hence, basis functions (r) become symbolic objects.
Can’t reconstruct electronic density from eigenvectors.
 Fitting TB parameters is a craft. A set of target properties
to reproduce is selected from experiment, from ab initio
or DFT calculations, or a combination of both.
 Initially, it is good practice to try and fit by hand to get a
feeling of the role of the various parameters. Then, since
parameters interact, it is better to refine them using some
automatic procedure (e.g. non-linear fitting, genetic, etc.)
EXAMPLE: TB BANDS IN A SOLID
 Simple cubic solid. Parameters: Onsite
0 and hopping t.
 Hoppings cut-off at first neighbors
 TB energy bands:
 (k )   0  2tcos(kx a)  cos(k y a)  cos(kz a)
 Onsite  Center of the band. Hopping  Band width
 TB works when atomic picture is a good start: TM, also sp
TIGHT-BINDING: ENERGY
 Being TB non-interacting electrons, the electronic energy
is just the sum of the occupied eigenvalues
N
ETB    n  Tr[ ˆHˆ TB ]
n 1
 This produces purely attractive forces. At short distances
these have to be compensated by Pauli’s repulsion forces
between overlapping electronic shells.
 The inter-nuclear (bare Coulomb) potential is modified into
a generic replusive potential Vrep(|RI-RJ|) to be fitted.
 Tight-binding band model:
1 P
ETB    n  Vrep  R I  R J
2 I J
n 1
N

 Onsite elements are related to charge on atoms, which are
given by eigenvectors of H. Points at introducing a charge
variable and imposing self-consistency (or LCN).
SELF-CONSISTENT TIGHT-BINDING
 Introduce charge inbalance for each atom, defined as:
M
qI  qI (out)  qI (in)    = Mulliken charge
 1
 A functional of second order in the charge inbalance can be
written as (Finnis):
N
ESCCT
1 P
1 P
2
   n  U I qI  U IJ qI qJ
2 I 1
2 I 1
n 1
Onsite Coulomb
(Hubbard U)
Resists charge
accumulation
Screening
UIJ = 1/|RI-RJ|
 Lowest level of self-consistency: self-consistent charge
transfer (SCCT) model.
SELF-CONSISTENT TIGHT-BINDING
 Self-consistency in the entire charge distribution can be
achieved by generalizing the second-order functional to
multipole moments (up to some reasonable order).
N
ESCTB
1 P
1 P
1 P
2
l 'm '
   n  U I qI  U IJ qI qJ   qlm (R I ) Blm
ql 'm' (R J )
2 I 1
2 I J
2 I  J lm l 'm'
n 1
 Terminating the multipolar expansion at l=1 gives a
polarizable tight-binding model (dipole polarizability)
 SCTB Hamiltonian:
IJ
IJ
Hˆ 
(SCTB)  Hˆ TB
,  U I qI   IJ  Vlm (R I )l IJ
lm
HOPPINGS
 Two-center integrals depend on type of orbital (SlaterKoster tables), and on relative orientation
 Example: Si (sp)
2 onsite: s and p
4 hoppings: ss, sp, ps, pp
 Example: Pt(spd)
3 onsite: s , p and d
13 hoppings: ss, sp, ps, pp,
sd, ds, pd, dp, pd, dp,
dd, dd, dd
HOPPINGS AND PAIR POTENTIAL
 Distance dependence: hoppings decay with distance in a
way that depends on the type of orbitals
 Ducastelle: exponential
 Harrison: Power-law
 Cut-off: Decay too slow. Apply cutoff to try and retain only
relevant neighbors (preferably only nearest).