Transcript LSR Midterm Report 2007

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
t1x.out
t1xi
t1xo
t1x
../../Psps_for_tests/01h.pspgth
t11.in
t11.out
t1xi
t1xo
t1x
../../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
ntypat 1
znucl 1
natom 2
typat 1 1
xcart
-0.7 0.0 0.0
0.7 0.0 0.0
ecut 10.0
nkpt 1
nstep 10
toldfe 1.0d-6
diemac 1.0
diemix 0.5
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
Cell size is 10**3
One type of atom
Atomic number is one
There are two atoms
They both are of type 1
Location of the atoms
atom 1, in Bohr
atom 2, in Bohr
Cut-off energy, in Hartree
One k point
Maximal number of SCF cycles
Threshold
Preconditioner
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:
– 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
Energy
-1.02
-1.04
-1.06
-1.08
-1.10
-1.12
0.00
0.20
0.40
0.60
0.40
0.60
Force
0.40
0.20
0.00
0.00
-0.20
0.20
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
#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
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
acell 10 10 10
ecut:? 10
ecut+? 5
#First dataset: bond length
natom?1 2
ionmov?1 3
ntime?1 10
tolmxf?1 5.0d-4
xcart?1 -0.7 0.0 0.0
0.7 0.0 0.0
toldff?1 5.0d-5
nband?1 1
#Second dataset: H-atom
natom?2 1
nsppol?2 2
occopt?2 2
nband?2 1 1
occ?2 1.0 0.0
toldfe?2 1.0d-6
xcart?2 0.0 0.0 0.0
spinat?2 0.0 0.0 1.0
15
Lesson 2: ecut
Estim ated Required Mem ory (MB)
What determines ecut?
What if H is changed to Cl?
100
75
50
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
acell:? 8 8 8
acell+? 2 2 2
ecut 10
#First dataset: bond length
natom?1 2
ionmov?1 3
ntime?1 10
tolmxf?1 5.0d-4
xcart?1 -0.7 0.0 0.0
0.7 0.0 0.0
toldff?1 5.0d-5
nband?1 1
#Second dataset: H-atom
natom?2 1
nsppol?2 2
occopt?2 2
nband?2 1 1
occ?2 1.0 0.0
toldfe?2 1.0d-6
xcart?2 0.0 0.0 0.0
spinat?2 0.0 0.0 1.0
17
Lesson 2: acell
Estimated Required Memory (MB)
What determines acell?
What if H is changed to Cl?
60
50
40
30
20
10
0
8
10
12
14
16
18
Bond (Bohr)
E-Atom ize
0.1720
1.5900
0.1680
1.5600
0.1640
1.5300
0.1600
0.1560
8
10
12
14
16
18
1.5000
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 ki-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
1.0
0.0
0.0
0.5
1.0
0.0 # L point
0.0 # Gamma point
0.5 # X point
1.0 # Gamma point in another cell.
• Results: kpt#
15
10
5
0
0
5
10
15
20
25
30
35
-5
-10
28
40
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
-57.1100
-57.1200
-57.1300
7.5400
E-05
7.5300
E-10
7.5200
A-05
A-10
7.5100
7.5000
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.
2.
NSC GS with nqpt=1, qpt, getwfk,
getden, kptopt=3, tolwfr, iscf4=-2
SC RF Phonon with rfphon=1, rfatpol,
rfdir, nqpt=1, qpt, getwfk=1, getwfq=4,
kptopt=3, tolvrs, iscf
Notice: splitting between TO and LO
AlAs Phonon Energy (H)
2.0E-03
1.8E-03
1.6E-03
1.4E-03
1.2E-03
1.0E-03
8.0E-04
6.0E-04
4.0E-04
2.0E-04
0.0E+00
G
0
L
0.5
G
1.5
X
1
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.
2.
3.
4.
•
•
•
name output
description
number of DDB files
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
– Elastic constants
– MD (Harmonic approximation)
0.25
0.20
0.15
0.10
0.05
0.00
-0.05 0
-0.10
2
4
6
8
10
12 14
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
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0
0
20
40
60
80
100
120
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.1384756E+03 8.1463588E+03 3.9416450E-01 1.4169102E+00
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
80.0 7.7974376E+03 8.7915524E+03 1.2426435E+01 1.9325165E+01
100.0 7.5004823E+03 9.2274431E+03 1.7269608E+01 2.4175005E+01
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
160.0 6.0396228E+03 1.1028289E+04 3.1179165E+01 3.4847422E+01
180.0 5.3732225E+03 1.1749439E+04 3.5423425E+01 3.7183863E+01
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
–
–
•
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.
–
–
–
–
•
Frequency dependent linear dielectric tensor
Frequency dependent second order nonlinear susceptibility tensor
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

QP
QP
  Vext r  VH r  nk (r)   dr   r,r ,  nk  nk (r )  nk nk (r)
 2

QP equation
 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',) 
2
GW
GW Self-Energy
 de
i 
G(r,r',  )W (r,r',)
Green
Function
Dynamical
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):
DFT
DFT *

r



(0)
nk
nk r'
G (r,r', ) 
DFT
DFT



i

sgn(

nk
nk  )
nk
(0) is the Kohn-Sham electronic structure:
The
basic
ingredient
of
G

65

Lesson 7: Dynamical Screened Interaction
W is approximated by RPA:
WG,G (q,)  1
G,G (q, )v G (q)
Dielectric Matrix
Coulomb
Interaction
Dynamical Screened Interaction
v G (q) 
4
2
q  G
RPA approximation
G,RPAG (q,)  G,G  vG (q) (0)
G,G (q, )

Adler-Wiser expression
 (0)
G,G (q, )  2  f n,k  f n',k q 
n,n ,k
Independent Particle Polarizability
i(q G)r
DFT
DFT i(q G')r
DFT
 nDFT
e


e

,k q
n,k
n,k
n  ,k q
DFT
DFT
n,k
 n',k
q    i
ingredients: KS wavefunctions and KS energies
66
Lesson 7: GWA correction to LDA
QP equation
 2

QP
QP
  Vext r  VH r  nk (r)   dr   r,r ,  nk  nk (r )  nk nk (r)
 2

KS equation
 2

DFT
  Vext r  VH r  nk (r)  Vxc r  nk (r)  nk  nk (r)
 2


Difference = Vxc is replaced by .
Thus GWA
correction to the DFT KS eigenvalues by 1st order PT:

0-order wavefunctions
QP
DFT
DFT
DFT
nk
 nk
 nk
(r,r',) Vxc (r) nk
0-order
Non Self-Consistent G0WRPA, 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 nonlocal, energy dependent and life-time effects)
• LDA generally provides a qualitative understanding.
69
Lesson 7: GWA Calculation Steps
GW calculation scheme
calculation of the
KS ELECTRONIC STRUCTURE

calculation of the
RPA SCREENING
W
calculation of the
GW CORRECTIONS
(SELF-ENERGY MATRIX ELEMENTS)
<-V>
1. 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
2. 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
3. Compute
– Self-energy sigma at the given kpoint, 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
• Two steps
optdriver?1 3
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
– Dielectric constant
– Normalization (grep sum_g)
– Band energies (grep –A 2 –i “E0 <vxclda”)
getscr?2
ecutwfn?2
ecutsigx
nband?2
-1
5.0
6.0
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
• Two steps
– Data Set ?1: calculate screening (1/ε)
– 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 <vxclda”)
• To achieve band energy error < 0.01 eV, how
many bands must be used?
nband11
nband21
nband31
nband41
nband51
25
50
100
150
200
optdriver?1
ecuteps
ecutwfn?1
3
6.0
4.0
optdriver?2
getscr?2
...
4
-1
80
Lesson 7: ecuteps convergence
• tgw_8.in: ndtset 10, udtset 5 2
• Two steps
– Data Set ?1: calculate screening (1/ε)
– 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 <vxclda”)
optdriver?1
nband?1
ecutwfn?1
ecuteps:?
ecuteps+?
optdriver?2
getscr?2
...
3
100
4.0
3.0
1.0
4
-1
• 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