Transcript New Castep Functionality

New Castep Functionality
Linear response and non-local
exchange-correlation functionals
Stewart Clark
University of Durham, UK
The Authors of Castep







Stewart Clark, Durham
Matt Probert, York
Chris Pickard, Cambridge
Matt Segal, Cambridge
Phil Hasnip, Cambridge
Keith Refson, RAL
Mike Payne, Cambridge
New Functionality

Linear response
1.
2.
3.
4.

Density functional perturbation theory
Atomic perturbations (phonons)
E-field perturbations (polarisabilities)
k-point perturbations (Born charges)
Non-local exchange-correlation
1.
2.
3.
Summary of XC functionals
Why use non-local functionals?
Some examples
Density Functional Perturbation
Theory


Based on compute how the total energy
responds to a perturbation, usually of the
DFT external potential v
Expand quantities (E, n, y, v)
E  E ( 0 )  E (1)  2 E ( 2 )  ...
• Properties given by the derivatives
E (n)
2
1  nE
1

E
( 2)

and
in
particular
E

n! n
2 2
The Perturbations

Perturb the external potential (from the ionic
cores and any external field):
•
•
•
•

Ionic positions
Cell vectors
Electric fields
Magnetic fields




phonons
elastic constants
dielectric response
NMR
But not only the potential, any perturbation to
the Hamiltonian:
•
•
d/dk
d/d(PSP)


Born effective charges
alchemical perturbation
Phonon Perturbations





For each atom i at a time,
in direction ax, y or z
This becomes perturbation
denoted by 
The potential becomes a
function of perturbation 
Take derivatives of the
potential with respect to 
Hartree, xc: derivatives of
potentials done by chain
rule with respect to n and 
uai  2 cos(qRi )
Da i  j (q) 
a i  j
a i  j (q)
mi m j
ij are atom labels. a are directions
 2E
(R ) 
is the force - constant matrix
Ra R 
i
j
So we need this E(2)
The expression for phonon-E(2)

(1)
(0)
(1)
(1)
(0)
(0) (1)
(1)
E (2)   y k,n
HKS
 k,(0)n y k,n
 y (1)
v
y

y
v
y
k,n
k, n
k, n
k,n
k,n
1
 EHxc (1)
(2)
(0)
 
n (r)n (1) (r)   y (0)
v
y
k,n
k,n
2 n(r)n(r )
k,n
2
•Superscripts denote the order of the perturbation
•E(2) given by 0th and 1st order wave functions and densities
•This is a variational quantity – use conjugate gradients minimiser
•Constraint: 1st order wave functions orthogonal to 0th order wave functions
•This expression gives the electronic contribution

The Variational Calculation

E(2) is variational with respect to |y(1)>
The plane-wave coefficients are varied to find
the minimum E(2) under a perturbation
•
•
•



of a given ion i
in a given direction a
and for a given q
Analogous to standard total energy
calculation
Based on a ground state (E(0)) calculation
Can be used for any q value
Sequence of calculation
a i  j (q)  2 Ea i  j (q)
( 2)
Find electronic force constant matrix
 Add in Ewald part
 Repeat for a mesh of q
And with Fourier interpolation:
 Fourier transform to get F(R)
 Fit and interpolate
 Fourier transform and mass weight to get D at any q

Phonon LR: For and against

For
• Fast, each wavevector component about the
•
•
•

same as a single point energy calculation
No supercells requires
Arbitrary q
General formalism
Against
• Details of implementation considerable
Symmetry Considerations




When perturbing the system the symmetry is broken
No time reversal symmetry
Implication is: k-point number increases
For example, Phonon-Si2 (Diamond):
6x6x6 MP set, 48 symmetry operations leads to SCF 28 k-points
q=(0,0,0), 48 symmetry elements
q=(1/2,0,0), 12 symmetry elements
x-displacement leaves 12 elements
y-displacement leaves 4 elements
72 k-points for E(2)
108 k-points for E(2)
q=arbitrary, leaves only identity element and needs 216 k-points
So what can you calculate?

…and phonon DOS
Phonon dispersion curves…
CASTEP Phonon Dispersion
CASTEP Density
Fr
Frequency (cm-1)
Density of Phonon States (1/cm-1)
600
0.010
600
500
0.009
500
400
0.008
400
0.007
300
0.006
300
Thermodynamics


Phonon density of states
Debye temperature

  

F V , T   E0 V   k BT  2 sinh 
 
 2 k BT  

Phase stability via

Entropic terms – derivatives of free energy
Vibrational specific heats

Electric Field Response






Bulk polarisability
Born effective charges
Phonon G-point LO/TO splitting
Dielectric permittivities
IR spectra
Raman Spectra
Why Electric Fields




Need Born effective
charges to get LO-TO
splitting: originates from
finite dipole per unit cell
Example given later…
Found from d/dk
calculation and similar
cross-derivative
expressions
Technical note: this means
that all expressions for
perturbed potentials
different at zone centre
than elsewhere
Examples of E-field Linear
Response
•Response of silicon to an electric
field perturbation
•Plot shows first order charge density
Blue: where electrons are removed
Yellow: where electrons go
Acknowledgement: E-field work by
my PhD student Paul Tulip
Phonon Dispersion of an Ionic
Solid: G-point problems

At zone
centre: Finite
dipole (hence
E-field) per
unit cell
caused by
atomic
displacements
E-field on a polar system: NaCl


Without E-field: LO/TO Frequencies: 175 cm-1
Born charge: Na
Cl
0.00 
  1.06 0.00


0
.
00

1
.
06
0.00 

 0.00
0.00  1.06 

 1.06 0.00 0.00 


 0.00 1.06 0.00 
 0.00 0.00 1.06 


Ionic character as expected: Na+ and Cl- ions
LO/TO splitting is 90cm-1: Smooth dispersion curve at the G-point
•Polarisability Tensor
•Electronic
Permittivity Tensor
 33.4 0.0 0.0 


 0.0 33.4 0.0 
 0.0 0.0 33.4 


 2.68 0.0 0.0 


 0.0 2.68 0.0 
 0.0 0.0 2.68 


Quantities given in
atomic units
Summary of Density Functional
Perturbation Theory












Phonon frequencies
Phonon DOS
Debye Temperate
PVT phase diagrams
Vibrational specific heats
Born effective charges
LO/TO Splitting
Bulk polarisabilities
Electric permittivities
First order charge densities (where electrons move from and to)
Etc…
Lots of new physics…and more planned in later releases
Non-local XC functionals





Basic background of XC interaction
Description within LDA and GGAs
Some non-local XC functional
Implementation within Castep
Model test cases
The exchange-correlation
interaction
Many body Hamiltonian – many body wavefunction
eZ j
Zi Z j
1 2
1
e2
1
H    
 
 
2
2 i  j | ri  rj | 2 i  j | Ri  R j |
i , j | ri  R j |
Kohn-Sham Hamiltonian – single particle wavefunction
Zi Z j
Z i n(r )dr 1 n(r )
1 2
1
H    
 
dr   
  xc [n(r )]

2
| Ri  r | 2 | r  r |
2 i  j | Ri  R j |
i
Single particle
KE – not many
body
Hartree requires
self-interaction
correction
All goes
in here
Approximations to XC






Local density approximation – only one
Generalised gradient approximations – lots
GGA’s + Laplacians – no real improvements
Meta-GGA’s (GGA + Laplacian + KE) – many
recent papers – but nothing exciting yet.
Hybrid functionals (GGA/meta-GGA + some exact
exchange from HF calculations) – currently
favoured by the chemistry community.
Exc[n(r,r’)] – has been generally ignored recently
by DFT community (although QMC and GW show
it to have several interesting properties!).
New Functionals: The CPU cost





We aim to go beyond the local (LDA) or semilocal (GGA) approximations
Why? Cannot get any higher accuracy with
these class of functionals
The computational cost is high: scaling will be
O(N2) or O(N3)
This is because all pairs (or more) of electrons
must be considered
That is, beyond the single particle model most
DFT users are familiar with
First Class of New XC Functional:
Exact and Screened Exchange

Exact Exchange (Hartree-Fock):
y i* (r )y *j (r )y j (r )y i (r )
1 N N
E x [{y i }]     dr  dr 
2 i 1 j 1
| r  r |

Screened Exchange:
|r  r |
y i (r )y j (r )e TF y j (r )y i (r )
1 N N
Esx [{y i }]     dr  dr 
2 i 1 j 1
| r  r |
*

*
k
And in terms of plane waves:
2
Esx [{ci ,k }] 

c*j ,q (G )c j ,q (G  G  G)
c (G) ci ,k (G) 


2
2
i , k G 
G
j , q ,G | q  k  G  G  |  kTF
*
i ,k
Why Exact/Screened Exchange?




Exact exchange gives us access to the
common empirical “chemistry” GGA functionals
such as B3LYP
Screened exchange can lead to accurate band
gaps in semiconductors and insulators, so
improved excitation energies and optical
properties
The cost: scales O(N3), N=number of pl. waves
LDA/GGA is 0.1% of calculation, EXX is 99.9%
Silicon Band Structure
Extra cost can be worthwhile
Some recent results (by my PhD student, Michael Gibson)
•
•
•
•
LDA band gap: 0.54eV
Exact Exchange band gap: 2.15eV
Screened Exchange band gap: 1.11eV
Experimental band gap: 1.12eV
Another approach: Some theory
on exchange-correlation holes…

The exact(!) XC energy within DFT can be written as
nxc (r , r )
1
E xc [n]   n(r )dr 
dr 
2
| r  r |
•Relationship defined as the Coulomb energy between an
electron and the XC hole nxc(r,r’)
•XC hole is described in terms of the electron pair-correlation
function
nxc (r , r )  n(r )[ g xc (r , r )  1]
Determines probability of finding an electron at position r’
given one exists at position r
Properties of the XC hole

Pauli exclusion principle: gxc obeys the
sum rule


n
(
r
,
r
)
d
r
 1
xc

•The size of the XC hole is exactly one electron –
mathematically, this is the Pauli exclusion principle
•The LDA and some GGAs obey this rule
•For a universal XC function (applicable without bias), it must
obey as many (all!) exact conditions as possible
The Weighted Density Method

In the WDA the pair-correlation function is approximated
by
g xc (r , r )  1  GW DA[| r  r  |; n~ (r )]
•The weighted density is fixed at each point by enforcing the
sum rule
W DA
~ (r )]dr   1


n
(
r
)
G
[|
r

r
|];
n

•This retains the non-locality of the function along with
the Coulomb-like integral for Exc[n]
The XC potential


XC potential is determined in the usual
manner (density derivative of XC energy)
We get 3 terms v(r )  v1 (r )  v2 (r )  v3 (r )
where
1
GW DA[| r  r  |; n~ (r )]
v1 (r )   n(r )
dr 
2
| r  r |
1
GW DA[| r  r  |; n~ (r )]
v2 (r )   n(r )
dr 
2
| r  r |
1
n(r ) GW DA[| r  r  |; n~ (r )]
v3 (r )   n(r )dr 
dr 
2
| r  r |
n ( r )
Example of WDA in action: An
inhomogeneous electron gas

Investigate inhomogeneous systems by applying an external
potential of the form
vext (r )  vq cos( q.r )
•Very accurate quantum Monte Carlo results which to
compare
•Will have no pseudopotential effects
•It’s inhomogeneous!
•Given a converged plane wave basis set, we are testing the
XC functional only – nothing else to consider
Investigate XC-hole shapes

XC holes for reference electron at a density maximum
WDA results: acknowledgement to my PhD student Phil Rushton
But at the low density points…

Exchange-correlation hole at a density
minimum
How accurate are the XC
holes?

Compare to
VMC data – M.
Nekovee, W. M.
C. Foulkes and
R. J. Needs PRL
87, 036401
(2001)
Self Interaction Correction



H2+ molecule
contains only
one electron
HF describes it
correctly, DFT
with LDA fails
Self-interaction
correction is the
problem
Realistic Systems - Silicon


Generate self-consistent silicon charge density
Examine XC holes at various points
Si [110] plane
XC holes in silicon
Interstitial Region
Bond Centre Region
XC hole of electron
moving along [100]
direction
Some electronic properties

Band structure of Silicon – band gap opens
III-V semiconductor band gaps
Band gaps in eV
Si
Experiment
LDA
PW91
WDA

Ge
1.17
0.52
0.56
1.11
0.74
0.02
0.55
0.81
GaAs
1.57
1.23
1.32
1.65
InAs
0.42
0.34
0.35
0.45
The non-local potential opens up the band gap
of some simple semiconductors
Timescales
The following is my list of developments - other CGD members
have other new physics and new improvements







First implementation of phonon linear response already in Castep v2.2
Much faster (3-4 times faster) phonon linear response due to
algorithmic improvements in Castep v3.0
Phonon linear response calculations for metals under development.
Aimed to be in v3.0 or v3.1
Electric field response (Born charges, bulk polarisabilities, permittivity,
LO/TO splitting) in Castep v3.0
Exact and screened exchange: currently being tested and developed
aimed at the v3.0 release
Weighted density method: currently working and tested – release
schedule under discussion
Raman and IR intensities currently being implemented. Scientific
evaluation still be performed before a possible release can be
discussed