Ab initio Simulation of Nanostructures Workshop Sornthep Presentation by HC120421094435

VIEWS: 66 PAGES: 84

									Abinit Workshop

Sornthep Vannarat




                    1
   Lesson 1

Hydrogen Molecule




                    2
Lesson 1
•   cd ~abinit/Tutorial
•   mkdir work
•   copy ../t1x.files .
•   edit t1x.files
•   copy ../t11.in .
•   ../../abinis < t1x.files > log


../t11.in                         t11.in
t1x.out                           t11.out
t1xi                              t1xi
t1xo                              t1xo
t1x                               t1x
../../Psps_for_tests/01h.pspgth   ../../Psps_for_tests/01h.pspgth




                                                              3
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.pspgth




                                                           4
Lesson 1: Input parameters
• Input file: t11.in
acell 10 10 10     #   Cell size is 10**3
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




                                                       5
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
                                                              6
Lesson 1: Inter-atomic distance (1)
• 3 approaches:                                             Energy
   – compute total energy E(d) or force F(d)   -1.02
   – Use relaxation                            -1.04
                                               -1.06
• Multiple datasets                            -1.08
                                               -1.10
• t12.in: ndtset, xcart, xcart+, getwfk,       -1.12

  nband                                             0.00   0.20      0.40   0.60


• Edit t1x.files and run                                    Force

• Look at t12.out                              0.40
   – Data sets: symmetry and number of         0.20
     plane waves
   – Data sets: coordinates xangst, xcart,      0.00
                                                    0.00   0.20      0.40   0.60
     xred                                      -0.20
   – Data sets: Number of iterations
   – etotal and fcart
• Plot data                                                                 7
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




                                             8
Lesson 1: Charge density
• prtden = 1, t14.in
• move atoms to middle of box
• cut3d convert to OpenDX




                                9
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 =?

                                                               10
Lesson 1: Summary
• System
  – H2 molecule in a big box 10x10x10 Bohr^3
• Method
  – Using cut-off energy 10 Ha
  – LDA (LSDA for atom) with ixc=1 (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)


                                                        11
   Lesson 2

Convergence Study




                    12
Lesson 2: Combined calculation
t21.in
ndtset 2

acell 10 10 10
ecut 10
                           #Second dataset
#First dataset             natom2 1
natom1 2                   nsppol2 2 # LSDA
ionmov1 3 # BD algorithm   occopt2 2
ntime1 10                  nband2 1 1 # Spin up, down
tolmxf1 5.0d-4             occ2 1.0 0.0
xcart1 -0.7 0.0 0.0        toldfe2 1.0d-6
         0.7 0.0 0.0       xcart2 0.0 0.0 0.0
toldff1 5.0d-5             spinat2
nband1   1                   0.0 0.0 1.0 # Init spin


                                                   13
Lesson 2: Convergence
• Calculation parameters:
  – ecut
  – acell
  – number of k-points
  – convergence of SCF cycle: toldfe, toldff
  – convergence of geometry optimization: tolmxf




                                               14
Lesson 2: ecut
t22.in
ndtset 12 udtset 6 2          #Second dataset: H-atom
 acell 10 10 10                  natom?2 1
 ecut:? 10    ecut+? 5          nsppol?2 2
#First dataset: bond length     occopt?2 2
   natom?1 2                     nband?2 1 1
  ionmov?1 3                       occ?2 1.0 0.0
   ntime?1 10                   toldfe?2 1.0d-6
  tolmxf?1 5.0d-4                xcart?2 0.0 0.0 0.0
   xcart?1 -0.7 0.0 0.0         spinat?2 0.0 0.0 1.0
            0.7 0.0 0.0
  toldff?1 5.0d-5
   nband?1 1




                                                        15
Lesson 2: ecut
               Estim ated Required Mem ory (MB)




 100
                                                             What determines ecut?
 75

 50
                                                             What if H is changed to Cl?
 25

  0
       0       10        20        30      40     50    60



                              E-Atom ize
                                                                          Bond (Bohr)
 0.1800
                                                             1.5300
 0.1760
                                                             1.5000
 0.1720
                                                             1.4700
 0.1680

                                                             1.4400
 0.1640
           0        10        20     30     40     50   60
                                                                      0     20          40   60




                                                                                                  16
Lesson 2: acell
t23.in
ndtset 12 udtset 6 2          #Second dataset: H-atom
 acell:? 8 8 8                   natom?2 1
 acell+? 2 2 2                  nsppol?2 2
 ecut 10                        occopt?2 2
#First dataset: bond length      nband?2 1 1
   natom?1 2                       occ?2 1.0 0.0
  ionmov?1 3                    toldfe?2 1.0d-6
   ntime?1 10                    xcart?2 0.0 0.0 0.0
  tolmxf?1 5.0d-4               spinat?2 0.0 0.0 1.0
   xcart?1 -0.7 0.0 0.0
            0.7 0.0 0.0
  toldff?1 5.0d-5
   nband?1 1



                                                        17
Lesson 2: acell
                   Estimated Required Memory (MB)



 60

 50
                                                              What determines acell?
                                                              What if H is changed to Cl?
 40

 30

 20

 10

 0
      8       10           12            14         16   18



                           E-Atom ize                                      Bond (Bohr)

 0.1720
                                                              1.5900
 0.1680
                                                              1.5600
 0.1640

 0.1600                                                       1.5300

 0.1560
                                                              1.5000
          8        10           12         14       16   18
                                                                       8   10   12   14   16   18




                                                                                                    18
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




                                                             19
   Lesson 3

Crystalline Silicon




                      20
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
                                                  21
Lesson 3: Introduction
• Parameters:
  – rprim: premitive vectors
  – xred: reduced coordinates
  – K-point sampling
    • kptopt: 0 read from input, 1,2,3 generates,
      negative for band calculation
    • ngkpt: numbers of k-points in 3 directions
    • nshiftk shiftk kptrlatt
    • alternatively use kptrlen
    • Larger cell  smaller Brillouin zone

                                                    22
Lesson 3: Sample k-points generation
• ngkpt = 2,2,2
• First mesh has 8 points, (0,0,0), (0,0,½), (0,½,0), (0,½,½), (½,0,0),
  (½,0,½), (½,½,0), (½,½,½)
• nshiftk = 4, shiftk = (½,½,½), (½,0,0), (0,½,0), (0,0,½)
• First shift with (½,½,½), 8 k-points become (¼,¼,¼), (¼,¼,¾),
  (¼,¾,¼), (¼,¾,¾), (¾,¼,¼), (¾,¼,¾), (¾,¾,¼), (¾,¾,¾)
• Second shift with (½,0,0), 8 k-points become (¼,0,0), (¼,0,½),
  (¼,½,0), ...
• Third shift with (0,½,0), 8 k-points become (0,¼,0), (0,¼,½),
  (0,¾,0), ...
• Forth shift with (0,0,½), 8 k-points become (0,0,¼), (0,0, ¾),
  (0,½,¼), ...
• Value ki larger than ½ can be translated to k i-1 e.g. ¾  -¼
• Totally 32 points obtained by shifting first mesh, but can be
  reduced by symmetry
                                                                      23
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 3,2,2 etc.
  – nshiftk, shiftk
  – using kptrlen instead, with prtkpt = 1,0
                                                 24
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

                                                          25
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)




                                                26
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?

                                                    27
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                                           20



       •   0.5   0.0   0.0 # L point                        15



       •   0.0   0.0   0.0 # Gamma point                    10



       •   0.0   0.5   0.5 # X point                         5


       •   1.0   1.0   1.0 # Gamma point in another cell.    0
                                                                  0   5   10   15   20   25   30   35   40


• Results: kpt#                                              -5



                                                            -10


                                                                                                   28
    Lesson 4

Crystalline Aluminum
and Surface Energy



                       29
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


                                           30
Lesson 4: Smearing
• occopt = 4,5,6,7
  – 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
                                                          31
Lesson 4: Smearing delta functions




                                     32
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

                                         33
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


                                        34
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
 -57.0800                                  7.5700
            0   10   20   30   40          7.5600
 -57.0900
                                           7.5500
 -57.1000                                  7.5400
                                                                            A-05
                                    E-05   7.5300
 -57.1100                                                                   A-10
                                    E-10   7.5200
 -57.1200                                  7.5100
                                           7.5000
 -57.1300
                                           7.4900
 -57.1400                                           0   10   20   30   40




  Larger tsmear converges faster, but ...
    Try t43.in
                                                                            35
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

                                                        36
Lesson 4: Al surface energy
• Vacuum layer thickness
• Defining atomic positions in Cartesian
  coordinates is more convenient
• Preconditioner (dielng) for metal+vacuum
  case
• How many layers of vacuum are needed?
• t46.in


                                         37
Lesson 4: Al surface energy
• Al layer thickness
• Preconditioner (dielng) for metal+vacuum 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 iprcel=45
  – 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?
                                                           38
        Lesson 5

Dynamical Matrix, Dielectric
Tensor and Effective Charge



                               39
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).

                                                                     40
Lesson 5: Response functions
• ABINIT gives
  – phonon frequencies
  – electronic dielectric tensor
  – effective charges
  – Derivative DataBase (DDB)
    • Contains all 2DTEs and 3DTEs
    • MRGDDB
    • Anaddb



                                     41
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 (ipert=natom+1)
   – Homogeneous electric field (q=0), only (ipert=natom+2), idir =
     direction
• Homogeneous Strain
   – Uniaxial strain: ipert = natom+3, idir = 1,2,3 for xx,yy,zz
   – Shear strain ipert = natom+4, idir = 1,2,3 for yz, zx, xy
   – No internal coordinate relaxation

                                                                      42
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)




                                             43
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) = E0+xdE + x2d2E / 2 + ...
• x = change in Al position from its equilibrium

                                                        44
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,
                                                         45
Lesson 5: RF Full Dynamical Matrix
• At Gamma [q=(0,0,0)]
• Perturbation J(m,n):
   – m = Atom number (rfatpol 1 natom)
   – n = Direction (Reduced) (rfdir 1 1 1)
• Dynamical matrix M[j1,j2]
• Run
• Output files: *.out, 1WF, 1WF4, DDB,
   –   Perturbation of each atom is applied in each direction in turns
   –   idir 2,3 is symmetric with previous calculation
   –   ipert 1,2,4 (electric field)
   –   Note the symmetry M[i,j] = M[j,i]?
• Rerun with tolvrs = 10-18
• Phonon Energies
                                                                         46
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
• kptopt=1 for the ground state
• kptopt=2 for response functions at q=0
• kptopt=3 for response functions at non-zero q




                                                47
Lesson 5: Recipe: Steps
• Atomic displacement with q=0,
    1. SC GS IBZ (with kptopt=1)
    2. SC RF Phonon Half Set (with kptopt=2)
• Atomic displacement with q=k1-k2 (k1,k2 are special k-points),
    1. SC GS IBZ (with kptopt=1)
    2. SC RF Phonon Full Set (with kptopt=3)
• Atomic displacement for a general q point,
    1. SC GS IBZ (with kptopt=1)
    2. NSC GS k+q (might be reduced due to symmetries, with kptopt=1)
    3. SC RF Phonon Full Set (with kptopt=3)
• Electric Field (with q=0),
    1. SC GS IBZ (with kptopt=1)
    2. NSC RF DDK Half Set (with kptopt=2, and iscf=-3)
    3. SC RF EF Half Set (with kptopt=2)




                                                                        48
Lesson 5: Recipe: Combinations
• Full dynamical matrix, Dielectric tensor and Born
  effective charges
   1. SC GS IBZ (with kptopt=1)
   2. Three NSC RF DDK (one for each direction) Half Set (with
      kptopt=2, and iscf=-3)
   3. SC RF Phonon+EF Half Set (with kptopt=2)
• Phonon at q=0 and general q points
   – Perturbations at different q wavevectors cannot be mixed.
   1. SC GS IBZ (with kptopt=1)
   2. Three NSC RF DDK (one for each direction) Half Set (with
      kptopt=2, and iscf=-3)
   3. SC RF Phonon+EF Half Set (with kptopt=2)
   4. NSC GS k+q points (might be reduced due to symmetries, with
      kptopt=1)
   5. SC RF Phonon q0 Full Set (with kptopt=3)


                                                                 49
Lesson 5: Full RF calculation of AlAs
• Three Data Sets
   – SC GS IBZ
   – NSC RF DDK Half k-point set
   – SC RF Phonon+EF Half k-point set
• New parameters
   – rfelfd, getwfk, getddk
• trf1_5.in
• Run
• Output
   –   Dynamical matrix
   –   Dielectric tensor
   –   Effective charge
   –   Phonon energies

                                        50
Lesson 5: Multiple q Phonon
When qk1-k2,
   1.   NSC GS with nqpt=1, qpt, getwfk,                        AlAs Phonon Energy (H)

        getden, kptopt=3, tolwfr, iscf4=-2        2.0E-03
                                                  1.8E-03

   2.   SC RF Phonon with rfphon=1, rfatpol,      1.6E-03
                                                  1.4E-03

        rfdir, nqpt=1, qpt, getwfk=1, getwfq=4,   1.2E-03
                                                  1.0E-03

        kptopt=3, tolvrs, iscf                    8.0E-04
                                                  6.0E-04
                                                  4.0E-04
Notice: splitting between TO and LO               2.0E-04
                                                  0.0E+00
                                                            G
                                                            0        L
                                                                    0.5            X
                                                                                   1           G
                                                                                              1.5




                                                                                         51
        Lesson 6

Interatomic Force Constants,
Phonon and Thermodynamic
         Properties

                               52
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 occopt=2, 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
                                                        53
Lesson 6: Build DDB
• trf2_1.in has 10 Data Sets
  –   1st SC GS IBZ
  –   2nd NSC RF DDK Half Set
  –   3rd SC RF Phonon+EF Half Set (q=0)
  –   4th-10th SC RF Phonon Full Set (q=Δk)
• 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
                                              54
Lesson 6: Merging DDB Files
•    Steps to use Mrgddb
    1.   name output
    2.   description
    3.   number of DDB files
    4.   file name list
•    trf2_3.in
•    Run
•    Look at trf2_3.ddb.out
    – How many datasets?

                               55
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 < anaddb.files > log &
• Common parameters: dieflag, elaflag,
  elphflag, ifcflag, instrflag, nlflag, piezoflag,
  polflag, thmflag

                                                56
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
• ifcflag=1; (ifcflag=0 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 < trf2_4.files > trf2_4.log



                                                                   57
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:                      0.30

    – Phonon dispersion curve        0.25
                                     0.20
    – Elastic constants              0.15
    – MD (Harmonic approximation)     0.10
                                      0.05
                                      0.00
                                     -0.05 0   2   4   6   8   10   12 14
                                     -0.10
                                                   Distance

                                                                            58
Lesson 6: Phonon Band Structure
• ifcflag=1
• Q-point grid: brav, nqgpt, nqshft, q1shft
• Energy conservation and charge neutrality: asr,
  chneut
• Others: dipdip
• Band: eivec, nph1l, qph1l, nph2l, qph2l
• Run ..\..\anaddb < trf2_5.files > trf2_5.log
• trf2_5_band2eps.freq .dspl are obtained
• Run ..\..\band2eps < trf2_6.files > 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

                                                       59
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 < trf2_7.files > trf2_7.log
    0.08                                       At T F(J/mol-c)      E(J/mol-c)     S(J/(mol-c.K)) C(J/(mol-c.K))
    0.07
                                                (A mol-c is the abbreviation of a mole-cell, that is, the
                                                 number of Avogadro times the atoms in a unit cell )
    0.06                                         20.0 8.1384756E+03 8.1463588E+03 3.9416450E-01 1.4169102E+00
    0.05                                         40.0 8.1061319E+03 8.2368069E+03 3.2668767E+00 7.8985027E+00
                                                 60.0 7.9980215E+03 8.4575659E+03 7.6590737E+00 1.3992227E+01
    0.04
                                                 80.0 7.7974376E+03 8.7915524E+03 1.2426435E+01 1.9325165E+01
    0.03                                        100.0 7.5004823E+03 9.2274431E+03 1.7269608E+01 2.4175005E+01
    0.02                                        120.0 7.1069991E+03 9.7544363E+03 2.2061977E+01 2.8411187E+01
                                                140.0 6.6189292E+03 1.0359248E+04 2.6716563E+01 3.1955266E+01
    0.01
                                                160.0 6.0396228E+03 1.1028289E+04 3.1179165E+01 3.4847422E+01
      0                                         180.0 5.3732225E+03 1.1749439E+04 3.5423425E+01 3.7183863E+01
           0   20   40   60   80   100   120    200.0 4.6241912E+03 1.2512641E+04 3.9442249E+01 3.9069447E+01

                                                                                                            60
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
     –   Electro-optic coefficients



                                                                                    61
        Lesson 7

Quasi Particle Band Structure




                                62
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

                                                   63
    Lesson 7: GW Approximation
    In the quasiparticle (QP) formalism, the energies and wavefunctions are
    obtained by the Dyson equation:
       2                   
        Vext r  VH r  nk (r)   dr   r,r ,  nk  nk (r )  nk  nk (r)
                                                              QP                QP
                                                                                              QP equation
       2                    
     self-energy (a non-local and energy dependent operator) is the difference
    between the energies of bare particle and quasiparticle.
   Within the GW approximation, is given by:
                      i
         (r,r',) 
          GW
                                        de     i 
                                                           G(r,r',  )W (r,r', )
                     2
    GW Self-Energy                                      Green                    Dynamical
                                                       Function                   Screened
                                                                                 Interaction
                                                                                                    64
Lesson 7: Green function
Green function G corresponding to QP equation is




Green function G may be approximated by the independent particle G(0):

                                    nk r  nk r'
                                     DFT      DFT*
           G(0) (r,r', ) 
                           nk  - nk  i sgn(nk  )
                                   DFT             DFT




The basic ingredient of G(0) is the Kohn-Sham electronic structure:
 


                                                                         65
    Lesson 7: Dynamical Screened Interaction
    W is approximated by RPA:
     W G,G (q, )  1G (q, )v G (q)
                      G,                                  Dynamical Screened Interaction

                                                                            4
     Dielectric Matrix                 Coulomb               v G (q)             2   RPA approximation
                                      Interaction                         q  G

          G,G (q, )  G,G  v G (q) (0)G (q, )
           RPA
                                          G,
                                                 
                                                            Independent Particle Polarizability
    Adler-Wiser expression
                                                 n,k q ei(q G)r  n,k  n,k e i(q G')r  n,k q
                                                  DFT                   DFT   DFT                DFT

      (0)G (q,)  2  f n,k  f n',k q 
      G,
                      n,n ,k                                 n,k  n',k q    i
                                                               DFT    DFT



                                ingredients: KS wavefunctions and KS energies

                                                                                                           66
    Lesson 7: GWA correction to LDA
                     2                   
QP equation           Vext r  VH r  nk (r)   dr   r,r ,  nk  nk (r )  nk  nk (r)
                                                                            QP                QP

                     2                    
                     2                   
KS equation           Vext r  VH r  nk (r)  Vxc r  nk (r)  nk  nk (r)
                                                                           DFT

                     2                    
         
    Difference = Vxc is replaced by .
           
    Thus GWA correction to the DFT KS eigenvalues by 1st order PT:
                                                                0-order wavefunctions

             nk  nk  nk (r,r',) Vxc (r) nk
              QP    DFT   DFT                    DFT


                 0-order

     Non Self-Consistent G0W RPA, Plasmon Pole model

                                                                                                        67
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.




                                                  68
Lesson 7: Discrepancy of LDA
• In Kohn-Sham theory, eigenvalues εi are Lagrange
  multipliers to ensure the orthogonality of KS orbitals
• So both KS eigenvalues and orbitals are not physical
• εi are not energy levels; εN (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.

                                                        69
Lesson 7: GWA Calculation Steps
     GW calculation scheme      1. SC GS (fixed lattice parameters and
                                   atomic positions)
        calculation of the         – self-consistent density, potential and
  KS ELECTRONIC STRUCTURE            Kohn-Sham eigenvalues and
                                  eigenfunctions at relevant k-points
                                     and on a regular grid of k-points
                                2. Compute
       calculation of the
       RPA SCREENING               – susceptibility matrix chi0 and chi, on
               W                     a regular grid of q-points, for at least
                                     two frequencies (zero and a pure
                                     imaginary frequency ~ a dozen of
        calculation of the           eV)
       GW CORRECTIONS              – Dielectric matrix epsilon and
(SELF-ENERGY MATRIX ELEMENTS)        1/epsilon
             < -V>
                                3. Compute
                                   – Self-energy sigma at the given k-
                                     point, and derive the GW
                                     eigenvalues for the target states at
                                     this k-point

                                                                          70
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 10*1            #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

                                                       71
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

                                                        72
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


                                                                     73
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 <vxclda>”


                                                               74
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




                                                                 75
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 <vxclda”)
• If ecutwfn 5.0 is used, what is the error in band energy?

                                                               76
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 <vxclda”)
• What is the appropriate ecutsigmax and the associated
  error in band energy?
                                                      77
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 <vxclda”)
• What is the appropriate nband and the associated error
  in band energy?

                                                       78
Lesson 7: 1/ε ecutwfn convergence
• tgw_6.in: ndtset 10, udtset 5 2                 optdriver?1 3
• Two steps                                       ecuteps      6.0
   – Data Set ?1: calculate screening (1/ε)       nband?1      25
   – Data Set *2: calculate GW correction (Sigma) ecutwfn:?   3.0
• How to prepare KSS and SCR files?               ecutwfn+?   1.0
• Run
• Output                                          optdriver?2 4
    – Num. plane-waves for wave function in Sigma
      and Epsilon calculations                      getscr?2    -1
    – Dielectric constant                           ecutwfn?2   5.0
    – Normalization (grep sum_g)                    ecutsigx     6.0
    – Band energies (grep –A 2 –i “E0 <vxclda”)     nband?2     100
• If ecutwfn 4.0 is used, what is the error in
  band energy?




                                                                 79
Lesson 7: 1/ε nband convergence
• tgw_7.in: ndtset 10, udtset 5 2                   nband11   25
• Two steps                                         nband21   50
   – Data Set ?1: calculate screening (1/ε)
                                                    nband31   100
   – Data Set *2: calculate GW correction (Sigma)
• How to prepare KSS and SCR files?                 nband41   150
• Run                                               nband51   200
• Output
   – Num. plane-waves for wave function in Sigma    optdriver?1     3
     and Epsilon calculations                       ecuteps         6.0
   – Dielectric constant                            ecutwfn?1       4.0
   – Normalization (grep sum_g)
   – Band energies (grep –A 2 –i “E0 <vxclda”)      optdriver?2     4
• To achieve band energy error < 0.01 eV, how       getscr?2        -1
  many bands must be used?
                                                    ...



                                                                    80
Lesson 7: ecuteps convergence
• tgw_8.in: ndtset 10, udtset 5 2
• Two steps                                         optdriver?1   3
   – Data Set ?1: calculate screening (1/ε)
                                                    nband?1       100
   – Data Set *2: calculate GW correction (Sigma)
                                                    ecutwfn?1     4.0
• How to prepare KSS and SCR files?
• Run
                                                    ecuteps:?      3.0
• Output
   – Num. plane-waves for wave function in Sigma
                                                    ecuteps+?      1.0
     and Epsilon calculations
   – Dimension of epsilon                           optdriver?2   4
   – Dielectric constant                            getscr?2      -1
   – Normalization (grep sum_g)                     ...
   – Band energies (grep –A 2 –i “E0 <vxclda”)
• To achieve band energy error < 0.01 eV, how
  large ecuteps must be used?



                                                                  81
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
                                                                 82
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




                                                               83
Thank you




            84

								
To top