Band gap of ZnO

[okuno]

Dear Yambo Users and developers.

I'm new comer of Yambo users.

I want to calculate the correct band gap of the ZnO,
and try the GW calculation with Yambo.

However, the calculated band gap is underestimated , about 0.8 ev (calc) and the experimental
gap is about 3.0eV.

I want to know what is wrong in my calculation.

My input file for PWSCF scf and nscf calculation are like below,

for scf PWSCF calculation

&control
calculation = 'scf'
pseudo_dir = './'
outdir = './'
prefix = 'ZnO'
forc_conv_thr = 1.0D-4
/
&system
ibrav=4
celldm(1)=6.14134490124
celldm(3)=1.60210593687
nat=4
ntyp=2
ecutwfc =30.0
ecutrho = 300.0
degauss=0.01
occupations='smearing'
degauss = 0.01
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
/
ATOMIC_SPECIES
Zn 65.39 zn.ncpp.UPF
O 15.9994 O.pz-mt.UPF
ATOMIC_POSITIONS {crystal}
Zn 0.3333333333333357 0.6666666666666643 0.0000000000000000
Zn 0.6666666666666643 0.3333333333333357 0.5000000000000000
O 0.3333333333330017 0.6666666666669983 0.3825000000000000
O 0.6666666666669983 0.3333333333330017 0.8825000000000000
K_POINTS {automatic}
4 4 4 0 0 0

and the nscf calculation are

&control
calculation = 'nscf'
pseudo_dir = './'
outdir = './'
prefix = 'ZnO'
wf_collect=.true.
/
&system
ibrav=4
celldm(1)=6.14134490124
celldm(3)=1.60210593687
nat=4
ntyp=2
ecutwfc =30.0
ecutrho = 300.0
degauss=0.01
occupations='smearing'
degauss = 0.01
nbnd=140
nosym=.false.
force_symmorphic=.true.
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
/
ATOMIC_SPECIES
Zn 65.39 zn.ncpp.UPF
O 15.9994 O.pz-mt.UPF
ATOMIC_POSITIONS {crystal}
Zn 0.3333333333333357 0.6666666666666643 0.0000000000000000
Zn 0.6666666666666643 0.3333333333333357 0.5000000000000000
O 0.3333333333330017 0.6666666666669983 0.3825000000000000
O 0.6666666666669983 0.3333333333330017 0.8825000000000000
K_POINTS {automatic}
8 8 8 0 0 0

and Yambo GW calculation are
gw0 # [R GW] GoWo Quasiparticle energy levels
xxvxc # [R XX] Hartree-Fock Self-energy and Vxc
ppa # [R Xp] Plasmon Pole Approximation
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
FFTGvecs= 1225 RL # [FFT] Plane-waves
EXXRLvcs= 28159 RL # [XX] Exchange RL components
XfnQPdb= "none" # [EXTQP Xd] Database
XfnQP_N= 1 # [EXTQP Xd] Interpolation neighbours
% XfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP Xd] E parameters (c/v)
%
% QpntsRXp
1 | 65 | # [Xp] Transferred momenta
%
% BndsRnXp
1 | 140 | # [Xp] Polarization function bands
%
NGsBlkXp= 200 RL # [Xp] Response block size
CGrdSpXp= 100.0000 # [Xp] [o/o] Coarse grid controller
% EhEngyXp
-1.000000 |-1.000000 | eV # [Xp] Electron-hole energy range
%
% LongDrXp
1.000000 | 0.000000 | 0.000000 | # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy
GfnQPdb= "none" # [EXTQP G] Database
GfnQP_N= 1 # [EXTQP G] Interpolation neighbours
% GfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v)
%
% GbndRnge
1 | 140 | # [GW] G[W] bands range
%
GDamping= 0.10000 eV # [GW] G[W] damping
QPreport= "kpbne0ees0" # [GW] QP info. Keys: kp/bn/xx/xc/s0/sq/e0/eq/ee/zf/ds/lm/lf
%QPkrange # [GW] QP generalized Kpoint/Band indices
1| 1| 1|24|
%

All input and output files are attached in this posts.

If I can get your suggestions, I'm very glad.

Sincerely

Yukihiro Okuno.

[Conor Hogan]

Dear Yukihiro,
I've no experience with ZnO, but there are several places to start.
First, I recommend you look carefully (if you havent done so already) at:
Dixit et al, J. Phys.: Condens. Matter 22 (2010) 125505, especially the discussion related to the states included in the pseudopotential. I notice they also find a 0.8eV discrepancy in their quasiparticle gap.
Second, note the use of diago* keywords in running the bands/nscf calculation: http://www.yambo-code.org/doc/p2y_direct.php, otherwise your unoccupied states will not be very well converged.
Third, if you are happy with everything so far - well, there are many many parameters to vary (in the yambo.in file) to check convergence. I have no way of knowing if the number of G-vectors or bands are enough for your system, although they don't look completely incorrect at first glance. Your large discrepancy may well be linked to the pseudopotential, so check as well as you can all the LDA stuff before doing tedious convergence checks with yambo.

Hope it helps for now!
Conor

[okuno]

Dear Dr. Hogan.

Thank you for your reply. I have changed the pseudopotential and also, in nscf PWSCF calculation,
I put diag_... option for unoccupied states. But the results are still inconsistent to the experimental gap.
(calculation gap is about 2.0 ev experimental gap is 3.4 eV)
I have chaged the energy cutoff for pwscf scf calculation, and also GbndRnge in yambo, the calculated results
are not changed.) There many papers on ZnO which is used GW calculation, and seems GW calculation should
be more correct ( for example Phys.Rev.B 52 R15977 or Phys.Rev.B 60 10754)
pseudo potential seems good in LDA level. ( Are there any pseudopotential for Zn ? I want to use the pseudo
potential which are used by other people)

GW calculation more sensitive to pseudo potential than in LDA case ? or what is wrong in my calculation.

Sincerely.

Yukihiro Okuno
Here my input file

PWSCF scf

&control
calculation = 'scf'
pseudo_dir = './'
outdir = './'
prefix = 'ZnO'
forc_conv_thr = 1.0D-4
/
&system
ibrav=4
celldm(1)=6.14134490124
celldm(3)=1.60210593687
nat=4
ntyp=2
ecutwfc =50.0
ecutrho = 500.0
degauss=0.01
occupations='smearing'
degauss = 0.01
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
/
ATOMIC_SPECIES
Zn 65.39 zn.ncpp.UPF
O 15.9994 O.pz-mt.UPF
ATOMIC_POSITIONS {crystal}
Zn 0.3333333333333357 0.6666666666666643 0.0000000000000000
Zn 0.6666666666666643 0.3333333333333357 0.5000000000000000
O 0.3333333333330017 0.6666666666669983 0.3825000000000000
O 0.6666666666669983 0.3333333333330017 0.8825000000000000
K_POINTS {automatic}
4 4 4 0 0 0

PWSCF non-scf

&control
calculation = 'nscf'
pseudo_dir = './'
outdir = './'
prefix = 'ZnO'
wf_collect=.true.
/
&system
ibrav=4
celldm(1)=6.14134490124
celldm(3)=1.60210593687
nat=4
ntyp=2
ecutwfc =50.0
ecutrho = 500.0
degauss=0.01
occupations='smearing'
degauss = 0.01
nbnd=200
nosym=.false.
force_symmorphic=.true.
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
diago_thr_init=1.0e-6
diago_full_acc=.true.

/
ATOMIC_SPECIES
Zn 65.39 Zn.cpi.UPF
O 15.9994 O.pz-mt.UPF
ATOMIC_POSITIONS {crystal}
Zn 0.3333333333333357 0.6666666666666643 0.0000000000000000
Zn 0.6666666666666643 0.3333333333333357 0.5000000000000000
O 0.3333333333330017 0.6666666666669983 0.3825000000000000
O 0.6666666666669983 0.3333333333330017 0.8825000000000000
K_POINTS {automatic}
8 8 8 0 0 0

Yambo input

xxvxc # [R XX] Hartree-Fock Self-energy and Vxc
ppa # [R Xp] Plasmon Pole Approximation
gw0 # [R GW] GoWo Quasiparticle energy levels
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
FFTGvecs= 2527 RL # [FFT] Plane-waves
EXXRLvcs= 60707 RL # [XX] Exchange RL components
XfnQPdb= "none" # [EXTQP Xd] Database
XfnQP_N= 1 # [EXTQP Xd] Interpolation neighbours
% XfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP Xd] E parameters (c/v)
%
% QpntsRXp
1 | 65 | # [Xp] Transferred momenta
%
% BndsRnXp
1 | 200 | # [Xp] Polarization function bands
%
NGsBlkXp= 256 RL # [Xp] Response block size
CGrdSpXp= 100.0000 # [Xp] [o/o] Coarse grid controller
% EhEngyXp
-1.000000 |-1.000000 | eV # [Xp] Electron-hole energy range
%
% LongDrXp
1.000000 | 0.000000 | 0.000000 | # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy
GfnQPdb= "none" # [EXTQP G] Database
GfnQP_N= 1 # [EXTQP G] Interpolation neighbours
% GfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v)
%
% GbndRnge
1 | 200 | # [GW] G[W] bands range
%
GDamping= 0.10000 eV # [GW] G[W] damping
QPreport= "kpbne0ees0" # [GW] QP info. Keys: kp/bn/xx/xc/s0/sq/e0/eq/ee/zf/ds/lm/lf
%QPkrange # [GW] QP generalized Kpoint/Band indices
1| 1| 1|25|
%

[claudio]

Dear Yukihiro
the band gap of ZnO is still an open question. Have a look to the paper PRB 81 125207(2010)
in table 2 they report the results from different "GW" calculations and they differ by more than 1 ev!
First of all it seems very important to include semi-core states in the calculation (see table I), and second
even in this way it seems that the simple GoWo is not able to reproduce the experimental results, see table I of that paper.

Different strategies have been used to get a result close to the experimental value, partial self-consistency, hybrid-functional+GW, etc...
It is up to you to judge the reliability of these approaches, and also their degree of convergence.

Claudio

[okuno]

Dear Dr. Claudio.

Thank you for you suggestion.

I see your suggested paper. It seems some discrepancy on the results of calculation
and experiments.

ZnO is basically transparent oxyde. Then,
the band gap must be larger than about 3.0 eV.
So, I think my result (2.0eV) is too small.
I also calculate another oxide, MgO which has
band gap 7.8eV as another lesson for calculate the GW method by Yambo.

My results of the calculation by Yambo and PWSCF is 6.6 eV (exp 7.8 eV)
and small compare to the experiment as in the case of ZnO.

Past MgO calculation seems very good results (Phys.Rev.B vol52 8788 (1995) )
(calculated value of this paper is about 7.7eV, but they use LMTO )

I think some systematic errors has occurs on the calculation of
oxides by Yambo. But I cant find what is wrong in my calculation.
(I have increased the unoccupied state, or energy cutoff for the LDA.
but the results does not improved.)

Plasmon pole approximation is not appropriate some oxides or
plasmon pole value effect band gap so much ? (or G0W0 calculation
is not suitable for calculating the correct band gap for transparent oxides ?)

My PWSCF and Yambo input for MgO is below.

Sincerely.

PWSCF scf
&control
calculation = 'scf'
restart_mode = 'from_scratch'
pseudo_dir = './'
outdir = './'
prefix = 'MgO'
forc_conv_thr = 1.0D-4
/
&system
ibrav=2
celldm(1)=7.96708476785
nat=2
ntyp=2
ecutwfc =60.0
ecutrho = 600.0
degauss=0.01
degauss = 0.01
nbnd = 10
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
/
ATOMIC_SPECIES
Mg 65.39 Mg.pz-bhs.UPF
O 15.9994 O.pz-mt.UPF
ATOMIC_POSITIONS {crystal}
Mg 0.0000000000000000 0.0000000000000000 0.0000000000000000
O 0.5000000000000000 0.5000000000000000 0.5000000000000000
K_POINTS {automatic}
6 6 6 0 0 0

# for Non-scf
&control
calculation = 'nscf'
pseudo_dir = './'
outdir = './'
prefix = 'MgO'
wf_collect=.true.
/
&system
ibrav=2
celldm(1)=7.96708476785
nat=2
ntyp=2
ecutwfc =60.0
ecutrho = 600.0
degauss=0.01
degauss = 0.01
nbnd = 150
nosym=.false.
/
&electrons
conv_thr = 1e-9,
mixing_mode = 'plain'
mixing_beta=0.01,
diago_thr_init=1.0e-6
diago_full_acc=.true.
/
ATOMIC_SPECIES
Mg 65.39 Mg.pz-bhs.UPF
O 15.9994 O.pz.cpi.UPF
ATOMIC_POSITIONS {crystal}
Mg 0.0000000000000000 0.0000000000000000 0.0000000000000000
O 0.5000000000000000 0.5000000000000000 0.5000000000000000
K_POINTS {automatic}
8 8 8 0 0 0

Yambo input

xxvxc # [R XX] Hartree-Fock Self-energy and Vxc
ppa # [R Xp] Plasmon Pole Approximation
gw0 # [R GW] GoWo Quasiparticle energy levels
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
FFTGvecs= 1411 RL # [FFT] Plane-waves
EXXRLvcs= 31283 RL # [XX] Exchange RL components
XfnQPdb= "none" # [EXTQP Xd] Database
XfnQP_N= 1 # [EXTQP Xd] Interpolation neighbours
% XfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP Xd] E parameters (c/v)
%
% QpntsRXp
1 | 29 | # [Xp] Transferred momenta
%
% BndsRnXp
1 | 150 | # [Xp] Polarization function bands
%
NGsBlkXp= 512 RL # [Xp] Response block size
CGrdSpXp= 100.0000 # [Xp] [o/o] Coarse grid controller
% EhEngyXp
-1.000000 |-1.000000 | eV # [Xp] Electron-hole energy range
%
% LongDrXp
1.000000 | 0.000000 | 0.000000 | # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138 eV # [Xp] PPA imaginary energy
GfnQPdb= "none" # [EXTQP G] Database
GfnQP_N= 1 # [EXTQP G] Interpolation neighbours
% GfnQP_E
0.000000 | 1.000000 | 1.000000 | # [EXTQP G] E parameters (c/v)
%
% GbndRnge
1 | 150 | # [GW] G[W] bands range
%
GDamping= 0.10000 eV # [GW] G[W] damping
QPreport= "kpbne0ees0" # [GW] QP info. Keys: kp/bn/xx/xc/s0/sq/e0/eq/ee/zf/ds/lm/lf
%QPkrange # [GW] QP generalized Kpoint/Band indices
1| 1| 1|30|
%

[claudio]

Dear Yukihiro

I think your input is correct, just a small comment,
in the PWSCF file for the non-self-consistent calculation
use

nosym=.false.
force_symmorphic=.true.

you will use more symmetries and save time.

About the result, I'm not an expert in calculation on large-gap oxides but as far as I know
there are some problem with the GoWo approximation.
So if you want make a fair comparison try to find some paper that performed calculations with pseudopotentials in plane-waves.

Moreover calculations with LMTO usually include effects due to core-valence interaction that are not present
if you use a pseudopotential with a large core, try to perform the same calculation but with a psedopotential
with less electrons in the core and check if this changes the result.

Cla

[Daniele Varsano]

Dear Yukihiro,

FFTGvecs= 1411 RL # [FFT] Plane-waves
EXXRLvcs= 31283 RL # [XX] Exchange RL components

may be it's right, but the FFTGvecs value looks to me a little bit small in comparison with th EXXRLvcs.
You can easily check it, looking at the orthonormalization of your wfs. The <n|n'> value is written in
the report for a few bands.

Cheers,

Daniele

[okuno]

Dear Claudio and Daniele.

Thank you for your reply and comment.

I read the J.Phys.Condens Matt vol(22) 125505 (2010) and Phys.Rev.B vol(81) 125207 (2010),
and it depends on the pseudo potential due to an inadequate treatment of exchange correlation
and seems Zn(+12) pseudo potential is not suitable.

Then I calculate with Zn(+20) for WZ-ZnO, and the results are improved a littile.

The results are (Zn(+20) case the energy cutoff for scf are 200 Ry, and K-points are
6 * 6 * 6) (WZ-ZnO exp gap are 3.4eV)
pseudo LDA gap: GW gap :
Zn(+12) 1.12 2.06
Zn(+20) 0.82 2.33

By the way, How the default value of FFTGvec and EXXRLvcs are determined by Yambo ?
Especially, are FFTGvec determined by the plane wave cutoff, and can I set other value
from the determined value of energy cutoff by SCF calculation ?

[claudio]

Dear Yukihiro

it is good that you are getting a result close to the PRB (81) 125207 (2010) and Phys. Rev. B 66, 125101 (2002)
respectively 2.6 eV and 2.45 eV.

About the cutoff on planewaves, FFTGvec is the total number of plane-waves you will use in you calculations,
while EXXRLvcs is the one used in the exchange term.

After you read the wave-function, you can choose the FFTGvec number in the setup "yambo -i -V RL".
Usually a number of plane-waves equal to 1/2 or 1/3 of the total wave-function(WF) plane waves
is sufficient to have converged and fast calculations, but you can make some tests to check the converge.

You can get the total number of planewaves in the WF doing "yambo -D".
If you do not change FFTGvec Yambo automatically will use all the plane waves present in the PWSCF wavefunction,
the same for the EXXRLvcs.

Cla