Abinit Workshop

1
Abinit Workshop

• Sornthep Vannarat

2
Lesson 1

• Hydrogen Molecule

3
Lesson 1

• cd abinit/Tutorial
• mkdir work
• copy ../t1x.files .
• edit t1x.files
• copy ../t11.in .
• ../../abinis lt t1x.files gt log

../t11.in t1x.out t1xi t1xo t1x ../../Psps_for_tes
ts/01h.pspgth
t11.in t11.out t1xi t1xo t1x ../../Psps_for_tests/
01h.pspgth
4
Lesson 1 Starting abinit

• Program name abinis.exe
• Input files t11.in and t1xi
• Output files t11.out and t1xo
• Temporary files t1x
• Pseudo potential files ../../Psps_for_tests/01h.ps
  pgth

5
Lesson 1 Input parameters

• Input file t11.in

acell 10 10 10 Cell size is 103 ntypat 1
One type of atom znucl 1
Atomic number is one natom 2 There
are two atoms typat 1 1 They both are
of type 1 xcart Location of the
atoms -0.7 0.0 0.0 atom 1 in Bohr 0.7
0.0 0.0 atom 2 in Bohr ecut 10.0
Cut-off energy in Hartree nkpt 1
One k point nstep 10 Maximal number of
SCF cycles toldfe 1.0d-6 Threshold diemac
1.0 Preconditioner diemix 0.5
Using standard preconditioner
for molecules in a big box
6
Lesson 1 Output

• Abinit version
• Input/output files
• Values of input parameters
• Data Set and Pseudo potential file
• Number of plan waves
• Iterations
• Stress tensor
• Eigen values
• Max/Min Electronic density
• Total energy
• Values of input parameters (after calculation)
• Log file interactive input more details of
  iterations

7
Lesson 1 Inter-atomic distance (1)

• 3 approaches
• compute total energy E(d) or force F(d)
• Use relaxation
• Multiple datasets
• t12.in ndtset xcart xcart getwfk nband
• Edit t1x.files and run
• Look at t12.out
• Data sets symmetry and number of plane waves
• Data sets coordinates xangst xcart xred
• Data sets Number of iterations
• etotal and fcart
• Plot data

8
Lesson 1 Inter-atomic distance (2)

• Use relaxation ionmove ntime tolmxf toldff
• Multiple datasets
• t13.in
• Edit t1x.files and run
• Look at t13.out
• BROYDEN STEP
• value of coordinates after relaxation xangst
  xcart xred

9
Lesson 1 Charge density

• prtden 1 t14.in
• move atoms to middle of box
• cut3d convert to OpenDX

10
Lesson 1 Atomization energy

• Eatomization (2EHatom EH2molecule) per
  molecule
• Caution
• Calculations with the same setting
• Spin nsppol spinat
• Degeneracy HOMO and LUMO (see lesson_4)
• Ground-state charge density NON-spherical
  automatic determination of symmetries should be
  disabled (nsym)
• For Hydrogen
• ground state is spherical (1s orbital)
• HOMO and LUMO have a different spin
• t15.in define occupation of each spin occopt
  and occ
• Output file Eigen values Max/Min spin
• Atomization energy

11
Lesson 1 Summary

• System
• H2 molecule in a big box 10x10x10 Bohr3
• Method
• Using cut-off energy 10 Ha
• LDA (LSDA for atom) with ixc1 (default)
• Pseudopotential from Goedecker-Hutter-Teter table
• Results
• Bond length 1.522 Bohr (1.401 Bohr)
• Atomisation energy 0.1656 Ha 4.506 eV (4.747
  eV)

12
Lesson 2

• Convergence Study

13
Lesson 2 Combined calculation
t21.in
ndtset 2 acell 10 10 10 ecut 10 First
dataset natom1 2 ionmov1 3 BD
algorithm ntime1 10 tolmxf1 5.0d-4 xcart1 -0.7
0.0 0.0 0.7 0.0 0.0 toldff1
5.0d-5 nband1 1
Second dataset natom2 1 nsppol2 2
LSDA occopt2 2 nband2 1 1 Spin up down occ2
1.0 0.0 toldfe2 1.0d-6 xcart2 0.0 0.0
0.0 spinat2 0.0 0.0 1.0 Init spin
14
Lesson 2 Convergence

• Calculation parameters
• ecut
• acell
• number of k-points
• convergence of SCF cycle toldfe toldff
• convergence of geometry optimization tolmxf

15
Lesson 2 ecut
t22.in
ndtset 12 udtset 6 2 acell 10 10 10 ecut 10
ecut 5 First dataset bond length
natom1 2 ionmov1 3 ntime1 10
tolmxf1 5.0d-4 xcart1 -0.7 0.0 0.0
0.7 0.0 0.0 toldff1 5.0d-5
nband1 1
Second dataset H-atom natom2 1
nsppol2 2 occopt2 2 nband2 1 1
occ2 1.0 0.0 toldfe2 1.0d-6 xcart2
0.0 0.0 0.0 spinat2 0.0 0.0 1.0
16
Lesson 2 ecut
What determines ecut What if H is changed to Cl
17
Lesson 2 acell
t23.in
ndtset 12 udtset 6 2 acell 8 8 8 acell 2
2 2 ecut 10 First dataset bond length
natom1 2 ionmov1 3 ntime1 10
tolmxf1 5.0d-4 xcart1 -0.7 0.0 0.0
0.7 0.0 0.0 toldff1 5.0d-5
nband1 1
Second dataset H-atom natom2 1
nsppol2 2 occopt2 2 nband2 1 1
occ2 1.0 0.0 toldfe2 1.0d-6
xcart2 0.0 0.0 0.0 spinat2 0.0 0.0 1.0
18
Lesson 2 acell
What determines acell What if H is changed to Cl
19
Lesson 2 Optimum parameters and GGA calculation

• Use the optimum ecut and acell to determine H2
  bond length and atomization energy.
• Switch to GGA calculation by changing ixc
• No need to change pseudo-potential for H (small
  core)
• No need to change ecut
• No need to change acell
• Compare results

20
Lesson 3

• Crystalline Silicon

21
Lesson 3 Introduction

• Crystalline silicon (Diamond structure 2 FCC)
• the total energy
• the lattice parameter
• the band structure (actually the Kohn-Sham band
  structure)
• Parameters
• k-points
• smearing of cut-off

22
Lesson 3 Introduction

• Parameters
• rprim premitive vectors
• xred reduced coordinates
• K-point sampling
• kptopt 0 read from input 123 generates
  negative for band calculation
• ngkpt numbers of k-points in 3 directions
• nshiftk shiftk kptrlatt
• alternatively use kptrlen
• Larger cell smaller Brillouin zone

23
Lesson 3 Sample k-points generation

• ngkpt 222
• First mesh has 8 points (000) (00½)
  (0½0) (0½½) (½00) (½0½) (½½0)
  (½½½)
• nshiftk 4 shiftk (½½½) (½00) (0½0)
  (00½)
• First shift with (½½½) 8 k-points become
  (¼¼¼) (¼¼¾) (¼¾¼) (¼¾¾) (¾¼¼)
  (¾¼¾) (¾¾¼) (¾¾¾)
• Second shift with (½00) 8 k-points become
  (¼00) (¼0½) (¼½0) ...
• Third shift with (0½0) 8 k-points become
  (0¼0) (0¼½) (0¾0) ...
• Forth shift with (00½) 8 k-points become
  (00¼) (00 ¾) (0½¼) ...
• Value ki larger than ½ can be translated to ki-1
  e.g. ¾ -¼
• Totally 32 points obtained by shifting first
  mesh but can be reduced by symmetry

24
Lesson 3 Silicon crystal

• Look at input file t31.in meaning of acell
  rprim xred kptopt ngkpt nshiftk shiftk
  diemac
• Try running and check the result
• Try changing
• kptopt to 3
• ngkpt to 322 etc.
• nshiftk shiftk
• using kptrlen instead with prtkpt 10

25
Lesson 3 Silicon k-point convergence

• Look at t32.in and try running it
• Why problem occurs
• Change t32.in to correct the problem and to
  perform a series of calculations to test
  convergence against number of k-points
• ndtset 4
• ngkpt1 2 2 2
• ngkpt2 4 4 4
• ngkpt3 8 8 8
• ngkpt4 16 16 16
• Note number of k-points and energy convergence
• Convergence of wavefunction and charge density
  can also be verified

26
Lesson 3 Silicon k-point convergence

• Test k-point when begin the study new material
• Test (at least) three efficient k-point sets
• CPU time is proportional to number of k-points
• Symmetry reduce number of k-points but need to
  be weighted (wtk)
• Reference Monkhorst and Pack paper Phys. Rev. B
  13 5188 (1976)

27
Lesson 3 Silicon Lattice Parameter

• Parameters (t34.in)
• optcell 1
• ionmov 3
• ntime 10
• dilatmx 1.05
• ecutsm 0.5
• Experimental value 5.431 Angstrom at 25 degree
  Celsius see R.W.G. Wyckoff Crystal structures
  Ed. Wiley and sons New-York (1963)
• Calculated value
• Using LDA with the 14si.pspnc pseudopotential
• What are 2 data sets

28
Lesson 3 Silicon Band Structure

• Two steps
• SCF calculation of charge density
• Non-SCF calculation of eigen values (bands)
• Use L-Gamma-X-Gamma circuit
• In eight-atom cell coordinates (1/2 1/2 1/2)-(0
  0 0)-(1 0 0)-(1 1 1)
• In two-atom cell coordinates (1/2 0 0)- (0 0 0)-
  (0 1/2 1/2)-(1 1 1)
• Parameters (t35.in)
• prtden iscf getden nband kptopt ndivk
  enunit tolwfr
• kptbounds to
• 0.5 0.0 0.0 L point
• 0.0 0.0 0.0 Gamma point
• 0.0 0.5 0.5 X point
• 1.0 1.0 1.0 Gamma point in another cell.
• Results kpt

29
Lesson 4

• Crystalline Aluminum
• and Surface Energy

30
Lesson 4 Introduction

• Aluminum the bulk and the surface.
• the total energy
• the lattice parameter
• the relaxation of surface atoms
• the surface energy
• Smearing of the Brillouin zone integration
• Preconditioning the SCF cycle

31
Lesson 4 Smearing

• occopt 4567
• Use Fermi-Dirac when trying to mimic physical
  electronic temperature. It is less convenient to
  use due to long-tailed need more bands.
• In general Gaussian-like smearings are
  preferable.
• If you are interested only in the total energies
  you can just use a Gaussian smearing - but need
  to extract corrected energy by taking the semisum
  of the energy and the free energy.
• Methfessel-Paxton and Marzari-Vanderbilt do this
  automatically for you and also provide forces
  stresses and whatever else corrected for the
  leading term in the temperature.
• tsmear

32
Lesson 4 Smearing delta functions
33
Lesson 4 Bulk Al

• Look at t41.in
• ecut 6 Ha compare to previous cases
• H needed 30 Ha
• Si needed 8 Ha
• Run
• Look at output and note 2 points
• Components of energy
• Occupation of each band
• Test ecut convergence

34
Lesson 4 k-point convergence

• How to check k-point convergence
• Look at t42.in
• Run
• Look at the result
• Total energy
• acell
• Try with a different tsmear

35
Lesson 4 tsmear and k-point covergence
Aluminum Total energy (E) and Lattice parameters
(A) calculated using tsmear 0.05 0.10 as
functions of k-point grid
Larger tsmear converges faster but ...
Try t43.in
36
Lesson 4 Al (001) surface energy

• Slab calculation Al layer vacuum layer
• Thicknesses of Al and Vacuum layers
• Reference energy Bulk calculation with
  equivalent parameters i.e. cell shape k-point
  grid ecut
• Esurface (Eslab/nslab Ebulk/nbulk)/(2
  Asurface)
• Look at t44.in and t45.in what do they
  represent
• Difficulties surface reconstruction and
  different top-buttom surfaces

37
Lesson 4 Al surface energy

• Vacuum layer thickness
• Defining atomic positions in Cartesian
  coordinates is more convenient
• Preconditioner (dielng) for metalvacuum case
• How many layers of vacuum are needed
• t46.in

38
Lesson 4 Al surface energy

• Al layer thickness
• Preconditioner (dielng) for metalvacuum case
• Use an effective dielectric constant of about 3
  or 5
• With a rather small mixing coefficient 0.2
• Alternatively Use an estimation of the
  dielectric matrix governed by iprcel45
• Repeat the 3 aluminum layer case for comparison
• t47.in
• See t47_STATUS to check status of long
  calculation
• How many Al layers are needed

39
Lesson 5

• Dynamical Matrix Dielectric
• Tensor and Effective Charge

40
Lesson 5 Response functions

• Response functions are the second derivatives of
  total energy (2DTE) with respect to different
  perturbations e.g.
• phonons (Dynamical metrix)
• static homogeneous electric field (Dielectric
  tensor Born effective charges)
• strain (Elastic constants internal strain
  piezoelectricity)
• ABINIT computes FIRST-order derivatives of the
  wavefunctions (1WF)
• 2DTE is calculated from 1WF
• References
• X. Gonze Phys. Rev. B55 10337 (1997)
• X. Gonze and C. Lee Phys. Rev. B55 10355
  (1997).

41
Lesson 5 Response functions

• ABINIT gives
• phonon frequencies
• electronic dielectric tensor
• effective charges
• Derivative DataBase (DDB)
• Contains all 2DTEs and 3DTEs
• MRGDDB
• Anaddb

42
Lesson 5 Perturbations

• Phonon
• displacement of one atom (ipert) along one of the
  axis (idir) of the unit cell by a unit of length
  (in reduced coordinates
• characterized by two integer numbers and one
  wavevector
• rfatpol defines the set of atoms to be moved
• rfdir defines the set of directions to be
  considered
• nqpt qpt and qptnrm define the wavevectors to
  be considered
• Electric field
• DDK dH/dk auxiliary for RF-EF (ipertnatom1)
• Homogeneous electric field (q0) only
  (ipertnatom2) idir direction
• Homogeneous Strain
• Uniaxial strain ipert natom3 idir 123
  for xxyyzz
• Shear strain ipert natom4 idir 123 for
  yz zx xy
• No internal coordinate relaxation

43
Lesson 5 Ground State of AlAs

• trf1_1.in trf1_x.files (2 potentials)
• Note tolvrs 1.0d-18
• Run
• Is tolvrs reached (18)
• What is the total energy (15 digits)

44
Lesson 5 Frozen-phonon E

• tr1_2.in
• Read in previous wavefunction file (irdwfk 1)
• Al is moved (xred)
• Need to rename tr1_xo_WFK to tr1_xi_WFK
• Edit tr1_x.files to run tr1_2.in
• Run
• Compare tr1_1.out and tr1_2.out
• Symmetry K-points
• Cartesian forces
• RMS dE/dt
• Estimate E from E(x) E0xdE x2d2E / 2 ...
• x change in Al position from its equilibrium

45
Lesson 5 Response-Function E

• d2E/dx2 x change in Al position from its
  equilibrium
• tr1_3.in
• Read in previous wavefunction file (irdwfk 1)
• kptopt 2
• Atomic positions not changed (xred)
• Phonon perturbation (rfphon)
• Perturbation on Al atom (rfatpol)
• Direction (rfdir) and wave-vector (nqpt qpt)
• Run
• Output files .out (2DEtotal) 1WF DDB

46
Lesson 5 RF Full Dynamical Matrix

• At Gamma q(000)
• Perturbation J(mn)
• m Atom number (rfatpol 1 natom)
• n Direction (Reduced) (rfdir 1 1 1)
• Dynamical matrix Mj1j2
• Run
• Output files .out 1WF 1WF4 DDB
• Perturbation of each atom is applied in each
  direction in turns
• idir 23 is symmetric with previous calculation
• ipert 124 (electric field)
• Note the symmetry Mij Mji
• Rerun with tolvrs 10-18
• Phonon Energies

47
Lesson 5 Recipe K-Points

• Input k-point set for RF should NOT have been
  decreased by using spatial symmetries prior to
  the loop over perturbations
• ABINIT will automatically reduce k-points
• kptopt1 for the ground state
• kptopt2 for response functions at q0
• kptopt3 for response functions at non-zero q

48
Lesson 5 Recipe Steps

• Atomic displacement with q0
• SC GS IBZ (with kptopt1)
• SC RF Phonon Half Set (with kptopt2)
• Atomic displacement with qk1-k2 (k1k2 are
  special k-points)
• SC GS IBZ (with kptopt1)
• SC RF Phonon Full Set (with kptopt3)
• Atomic displacement for a general q point
• SC GS IBZ (with kptopt1)
• NSC GS kq (might be reduced due to symmetries
  with kptopt1)
• SC RF Phonon Full Set (with kptopt3)
• Electric Field (with q0)
• SC GS IBZ (with kptopt1)
• NSC RF DDK Half Set (with kptopt2 and iscf-3)
• SC RF EF Half Set (with kptopt2)

49
Lesson 5 Recipe Combinations

• Full dynamical matrix Dielectric tensor and Born
  effective charges
• SC GS IBZ (with kptopt1)
• Three NSC RF DDK (one for each direction) Half
  Set (with kptopt2 and iscf-3)
• SC RF PhononEF Half Set (with kptopt2)
• Phonon at q0 and general q points
• Perturbations at different q wavevectors cannot
  be mixed.
• SC GS IBZ (with kptopt1)
• Three NSC RF DDK (one for each direction) Half
  Set (with kptopt2 and iscf-3)
• SC RF PhononEF Half Set (with kptopt2)
• NSC GS kq points (might be reduced due to
  symmetries with kptopt1)
• SC RF Phonon q0 Full Set (with kptopt3)

50
Lesson 5 Full RF calculation of AlAs

• Three Data Sets
• SC GS IBZ
• NSC RF DDK Half k-point set
• SC RF PhononEF Half k-point set
• New parameters
• rfelfd getwfk getddk
• trf1_5.in
• Run
• Output
• Dynamical matrix
• Dielectric tensor
• Effective charge
• Phonon energies

51
Lesson 5 Multiple q Phonon

• When qk1-k2
• NSC GS with nqpt1 qpt getwfk getden
  kptopt3 tolwfr iscf4-2
• SC RF Phonon with rfphon1 rfatpol rfdir
  nqpt1 qpt getwfk1 getwfq4 kptopt3
  tolvrs iscf
• Notice splitting between TO and LO

52
Lesson 6

• Interatomic Force Constants
• Phonon and Thermodynamic Properties

53
Lesson 6 DDB File

• DDB contains dE with respect to 3 perturbations
  phonons electric field and stresses
• Header DDB version number natom nkpt nsppol
  nsym ntypat occopt and nband or array nband
  (nkpt nsppol) if occopt2 acell amu ecut
  iscf...
• Data
• Number of data blocks
• For each block Type of the block Number of
  Elements List of Elements
• In most cases each element consists of 4 integer
  and 2 real numbers idir1 ipert1 idir2 ipert2
  Re(2DTE) Im(2DTE)
• Symmetries may reduce number of elements
• DDB files can be merged by Mrgddb

54
Lesson 6 Build DDB

• trf2_1.in has 10 Data Sets
• 1st SC GS IBZ
• 2nd NSC RF DDK Half Set
• 3rd SC RF PhononEF Half Set (q0)
• 4th-10th SC RF Phonon Full Set (qk)
• Note need to overwrite default parameters in
  some Data Sets
• Data Sets 4-10 use selected K-Points which may be
  generated by trf2_2.in
• Run
• See trf2_1o_DS_DDB

55
Lesson 6 Merging DDB Files

• Steps to use Mrgddb
• name output
• description
• number of DDB files
• file name list
• trf2_3.in
• Run
• Look at trf2_3.ddb.out
• How many datasets

56
Lesson 6 ANADDB

• ANADDB analyses DDB for properties e.g. phonon
  spectrum frequency-dependent dielectric tensor
  thermal properties
• Files input output DDB other files
• To run anaddb lt anaddb.files gt log
• Common parameters dieflag elaflag elphflag
  ifcflag instrflag nlflag piezoflag polflag
  thmflag

57
Lesson 6 Interatomic Force Constants

• Dynamical Matrix and Interatomic Force Constants
  are Fourier Transforms of each other
• Calculated Dynamical Matrix on a grid of
  wavevectors IFC IFC vanishes rapidly with
  interatomic distance
• ifcflag1 (ifcflag0 is for checking or when
  there is not enough information in DDB)
• Q-point grid brav nqgpt nqshft q1shft
• Energy conservation and charge neutrality asr
  chneut
• Others dipdip ifcana ifcout natifc atifc
• Run ..\..\anaddb lt trf2_4.files gt trf2_4.log

58
Lesson 6 IFC Results

• On site term of Al trace 0.28080
• First NN 4 As atoms at 4.6 trace -0.06911
• Second NN 12 Al atoms at 7.5 trace 0.00062
• Third NN 12 As atoms at 8.8 trace -0.00037
• Fourth NN 6 Al atoms at 10.6 trace -0.00016
• Fifth NN 12 As atoms at 11.6 trace -0.00056
• Sixth NN Al atoms at 13.0 trace -0.00059
• Applications
• Phonon dispersion curve
• Elastic constants
• MD (Harmonic approximation)

59
Lesson 6 Phonon Band Structure

• ifcflag1
• Q-point grid brav nqgpt nqshft q1shft
• Energy conservation and charge neutrality asr
  chneut
• Others dipdip
• Band eivec nph1l qph1l nph2l qph2l
• Run ..\..\anaddb lt trf2_5.files gt trf2_5.log
• trf2_5_band2eps.freq .dspl are obtained
• Run ..\..\band2eps lt trf2_6.files gt trf2_6.log
• trf2_6.out.eps is obtained
• view with ghostview note discontinuity of
  Optical Phonon at Gamma point
• Edit trf2_5_band2eps.freq lines 1 and 31
  correct LO freq. (trf2_5.out)
• Run band2eps again

60
Lesson 6 Thermodynamic Properties

• Normalized phonon DOS
• Phonon internal energy free energy entropy
  constant volume heat capacity as a function of
  the temperature
• Debye-Waller factors (tensors) for each atom as
  a function of the temperature (DISABLED SORRY)
• Parameters thmflag ng2qpt ngrids q2shft
  nchan nwchan thmtol ntemper temperinc
  tempermin
• ..\..\anaddb lt trf2_7.files gt trf2_7.log

At T F(J/mol-c) E(J/mol-c)
S(J/(mol-c.K)) C(J/(mol-c.K)) (A mol-c is the
abbreviation of a mole-cell that is the
number of Avogadro times the atoms in a unit
cell) 20.0 8.1384756E03 8.1463588E03
3.9416450E-01 1.4169102E00 40.0
8.1061319E03 8.2368069E03 3.2668767E00
7.8985027E00 60.0 7.9980215E03
8.4575659E03 7.6590737E00 1.3992227E01
80.0 7.7974376E03 8.7915524E03 1.2426435E01
1.9325165E01 100.0 7.5004823E03
9.2274431E03 1.7269608E01 2.4175005E01
120.0 7.1069991E03 9.7544363E03
2.2061977E01 2.8411187E01 140.0
6.6189292E03 1.0359248E04 2.6716563E01
3.1955266E01 160.0 6.0396228E03
1.1028289E04 3.1179165E01 3.4847422E01
180.0 5.3732225E03 1.1749439E04
3.5423425E01 3.7183863E01 200.0
4.6241912E03 1.2512641E04 3.9442249E01
3.9069447E01
61
Other Response-Function Tutorials

• Optic Frequency-dependent linear and second
  order nonlinear optical response
• Frequency dependent linear dielectric tensor
• Frequency dependent second order nonlinear
  susceptibility tensor
• Electron-Phonon interaction and superconducting
  properties of Al.
• Phonon linewidths (lifetimes) due to the
  electron-phonon interaction
• Eliashberg spectral function
• Coupling strength
• McMillan critical temperature
• Elastic and piezoelectric properties.
• Rigid-atom elastic tensor
• Rigid-atom piezoelectric tensor (insulators only)
  
• Internal strain tensor
• Atomic relaxation corrections to the elastic and
  piezoelectric tensor
• Static non-linear properties
• Born effective charges
• Dielectric constant
• Proper piezoelectric tensor (clamped and relaxed
  ions)
• Non-linear optical susceptibilities
• Raman tensor of TO and LO modes

62
Lesson 7

• Quasi Particle Band Structure

63
Lesson 7 Introduction

• System
• Nucleus Electrons
• Approach
• Electron wave function
• Electron density DFT
• Quasiparticles
• Quasiparticle Bare particle Decorations
• Modify
• Equation of motion
• Energy Mass
• Life time

64
Lesson 7 GW Approximation
In the quasiparticle (QP) formalism the energies
and wavefunctions areobtained by the Dyson
equation
QP equation
S self-energy (a non-local and energy dependent
operator) is the difference between the energies
of bare particle and quasiparticle.
Within the GW approximationS is given by
GW Self-Energy
Dynamical Screened Interaction
Green Function
65
Lesson 7 Green function
Green function G corresponding to QP equation is
Green function G may be approximated by the
independent particle G(0)
The basic ingredient of G(0) is the Kohn-Sham
electronic structure
66
Lesson 7 Dynamical Screened Interaction
W is approximated by RPA
Dynamical Screened Interaction
Coulomb Interaction
Dielectric Matrix
RPA approximation
Independent Particle Polarizability
Adler-Wiser expression
ingredients KS wavefunctions and KS energies
67
Lesson 7 GWA correction to LDA
QP equation
KS equation
Difference Vxc is replaced by S. Thus GWA
correction to the DFT KS eigenvalues by 1st order
PT
0-order wavefunctions
0-order
Non Self-Consistent G0WRPA Plasmon Pole model
68
Lesson 7 GWA Performance
LDA GWA and experimental energy gaps for
semiconductors and insulators. GWA corrects
most of the LDA band gap underestimation. The
discrepancy for LiO2 results from the neglect of
excitonic effects. The experimental value for
BAS is tentative.
69
Lesson 7 Discrepancy of LDA

• In Kohn-Sham theory eigenvalues ei are Lagrange
  multipliers to ensure the orthogonality of KS
  orbitals
• So both KS eigenvalues and orbitals are not
  physical
• ei are not energy levels eN (highest level) is
  chemical potential for metal or negative
  ionization energy for semiconductor and insulator
• In absence of quasiparticle calculations. LDA
  energy are routinely used to interpreted
  experimental spectra
• LDA energy dispersions are often in fair
  agreement with experiment LDA band gaps are
  sometimes empirically adjusted to fit
  experimental values
• LDA VXC approximate self-energy (neglecting
  non-local energy dependent and life-time
  effects)
• LDA generally provides a qualitative
  understanding.

70
Lesson 7 GWA Calculation Steps

• SC GS (fixed lattice parameters and atomic
  positions)
• self-consistent density potential and Kohn-Sham
  eigenvalues and eigenfunctions at relevant
  k-points and on a regular grid of k-points
• Compute
• susceptibility matrix chi0 and chi on a regular
  grid of q-points for at least two frequencies
  (zero and a pure imaginary frequency a dozen of
  eV)
• Dielectric matrix epsilon and 1/epsilon
• Compute
• Self-energy sigma at the given k-point and
  derive the GW eigenvalues for the target states
  at this k-point

71
Lesson 7 Generation of KSS File

• tgw_1.in 3 Data Sets
• First
• nbandkss1 -1 Number of bands in KSS file
• -1 is full diagonalization see out file for
  number of plane waves and number of bands
• nband1 9 Number of bands to be computed
• istwfk1 101 Do not use time reversal
  symmetry for storing wavefunction
• npwkss 0 for same as ecut
• kssform 1 for full diag 3 for conjugated
  gradient
• symmorphi 0 symmorphic symmetry operations only

72
Lesson 7 Generation of SCR File

• Second
• optdriver2 3 Screening calculation
• getkss2 -1 Obtain KSS file from
  previous dataset
• nband2 17 Bands to be used in the
  screening calculation
• ecutwfn2 2.1 Cut-off energy of the
  planewave set to represent the wavefunctions
• ecuteps2 3.6 Cut-off energy of the
  planewave set to represent the dielectric matrix
• ppmfrq2 16.7 eV Imaginary frequency where
  to calculate the screening

73
Lesson 7 Calculation of Sigma

• Third
• optdriver3 4 Self-Energy calculation
• getkss3 -2 Obtain KSS file from
  dataset 1
• getscr3 -1 Obtain SCR file from
  previous dataset
• nband3 30 Bands for Self-Energy
  calculation
• ecutwfn3 5.0 Planewaves to represents
  wavefunctions
• ecutsigx3 6.0 Dimension of the G sum in
  Sigma_x
• Dimension of Sigma_c size of screening matrix
  (SCR file) or size of Sigma_x whichever is
  smaller
• nkptgw3 1 num of k-point for GW
  correction
• kptgw3 k-points which
  must present in KSS file
• bdgw3 4 5 calculate GW corrections
  for bands
• zcut for avoiding divergence in integration

74
Lesson 7 GWA Output File

• Data Set 1
• Kohn-Sham electronic Structure file
• Note number of plane waves number of bands
• Check Test on the normalization of the
  wavefunctions
• Data Set 2
• Check test on the normalization of the
  wavefunctions
• Is it the same as Data Set 1 Effect of ecutwfk
• total number of electrons per unit cell
• Electron density and plasma frequency
• calculating at frequencies omega eV
• dielectric constant
• Data Set 3
• Band energy E0 ltvxcldagt

75
Lesson 7 Convergence Study

• Simplify Gamma Point only
• tgw_2.in Generate KSS and SCR files
• Check Data Sets KSS is separated from GS
• Note values of ecut
• Run Check normalization number of electrons
  dielectric constant

76
Lesson 7 Sigma ecutwfn convergence

• tgw_3.in ndtset 5 ecutwfn 3.0 ecutwfn 1.0
• Note input KSS and SCR file names
• Rename
• tgw_2o_DS2_KSS to tgw_3o_DS1_KSS
• tgw_2o_DS3_SCR to tgw_3o_DS1_SCR
• Run
• Output
• Num. plane-waves for wave function in Sigma and
  Epsilon calculations
• Number of electrons per unit cell
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• If ecutwfn 5.0 is used what is the error in band
  energy

77
Lesson 7 ecutsigmax convergence

• tgw_4.in ndtset 7 ecutsigmax 3.0 ecutwfn 1.0
  
• Note input KSS and SCR file names
• Rename
• tgw_3o_DS1_KSS to tgw_4o_DS1_KSS
• tgw_3o_DS1_SCR to tgw_4o_DS1_SCR
• Run
• Output
• Num. plane-waves for Sigma_x calculations
• Number of electrons per unit cell
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• What is the appropriate ecutsigmax and the
  associated error in band energy

78
Lesson 7 Sigma nband convergence

• tgw_5.in ndtset 5 nband 50 nband 50
• Note input KSS and SCR file names
• Rename
• tgw_4o_DS1_KSS to tgw_5o_DS1_KSS
• tgw_4o_DS1_SCR to tgw_5o_DS1_SCR
• Run
• Output
• Num. plane-waves are fixed now
• Numbers of bands for KSS Sigma and Epsilon
• Number of electrons per unit cell
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• What is the appropriate nband and the associated
  error in band energy

79
Lesson 7 1/e ecutwfn convergence

• tgw_6.in ndtset 10 udtset 5 2
• Two steps
• Data Set 1 calculate screening (1/e)
• Data Set 2 calculate GW correction (Sigma)
• How to prepare KSS and SCR files
• Run
• Output
• Num. plane-waves for wave function in Sigma and
  Epsilon calculations
• Dielectric constant
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• If ecutwfn 4.0 is used what is the error in band
  energy

80
Lesson 7 1/e nband convergence

• tgw_7.in ndtset 10 udtset 5 2
• Two steps
• Data Set 1 calculate screening (1/e)
• Data Set 2 calculate GW correction (Sigma)
• How to prepare KSS and SCR files
• Run
• Output
• Num. plane-waves for wave function in Sigma and
  Epsilon calculations
• Dielectric constant
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• To achieve band energy error lt 0.01 eV how many
  bands must be used

81
Lesson 7 ecuteps convergence

• tgw_8.in ndtset 10 udtset 5 2
• Two steps
• Data Set 1 calculate screening (1/e)
• Data Set 2 calculate GW correction (Sigma)
• How to prepare KSS and SCR files
• Run
• Output
• Num. plane-waves for wave function in Sigma and
  Epsilon calculations
• Dimension of epsilon
• Dielectric constant
• Normalization (grep sum_g)
• Band energies (grep A 2 i E0 ltvxclda)
• To achieve band energy error lt 0.01 eV how large
  ecuteps must be used

82
Lesson 7 (direct) Egap of Silicon

• Data Set 1 SC GS
• print out the density
• 10 k-points in IBZ 4x4x4 FCC grid (Shifted no
  Gamma point)
• Data Set 2 NSC GS
• Kohn-Sham structure 19 k-points in IBZ but not
  shifted Gamma point included
• Data Set 3 Calculate Screening
• Very time-consuming
• ecutwfn 3.6
• nband 25 CPU time nband (Conduction)
• Accuracy of GW energy 0.2 eV
• Accuracy of energy difference 0.02 eV
• There is no zero of energy defined for bulk
  system
• Data Set 4 Self-energy matrix at Gamma

83
Lesson 7 (direct) Egap of Silicon 2

• What are direct Egap of Silicon by LDA and GWA
• Choice of pseudopotential can contribute to Egap
  variation
• Egap GWA accuracy 0.1 eV
• Full band calculation is possible by shifting
  the k-point
• Simplification GW corrections are quite linear
  with the energy

84
Thank you