VIEWS: 66 PAGES: 84 POSTED ON: 4/21/2012 Public Domain
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 q0 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 qk1-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 de 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 ei(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