Some problem of calculated exciton spectral levels

[Tianshu]

Dear Developers,

I want to calculate the exciton spectrum for black phosphorous, like energy levels of 1s, 2s, 2p, .... I have a experimental result as comparison. My calculation is as follows,

1. Calculate the corrected band gap with GW method (yambo -F gw.in -p p -g n -r -k hartree -V qp)
2. Adding the BSE effect to calculate the exciton peak (1s) and exciton binding energy (yambo -F bse.in -r -o b -p p -y d -k sex -V all)

After these, I have obtained the band gap and exciton binding energy which are consistent with the experimental values.

However,, when I want to get the energy levels of exciton spectrum, I use the command "ypp -J bse -e s", "ypp -F ypp.in -J bse -e a" and run it, I got results that can be plotted as the attached figure. From this figure we can see, the calculated 2s and 3s levels are all higher than 2eV. But the experimental results for 2s and 3s levels are 1.68 and 1.85eV, respectively.
The setting parameters for BSE calculation are as following, hope it can be helpful.

Code: Select all

#                                                           
#  __  __   ________   ___ __ __    _______   ______        
# /_/\/_/\ /_______/\ /__//_//_/\ /_______/\ /_____/\       
# \ \ \ \ \\::: _  \ \\::\| \| \ \\::: _  \ \\:::_ \ \      
#  \:\_\ \ \\::(_)  \ \\:.      \ \\::(_)  \/_\:\ \ \ \     
#   \::::_\/ \:: __  \ \\:.\-/\  \ \\::  _  \ \\:\ \ \ \    
#     \::\ \  \:.\ \  \ \\. \  \  \ \\::(_)  \ \\:\_\ \ \   
#      \__\/   \__\/\__\/ \__\/ \__\/ \_______\/ \_____\/   
#                                                           
#                                                           
#              GPL Version 4.1.3 Revision 76                
#                      Serial Build                         
#                http://www.yambo-code.org                  
#
rim_cut                      # [R RIM CUT] Coulomb potential
optics                       # [R OPT] Optics
ppa                          # [R Xp] Plasmon Pole Approximation
bss                          # [R BSS] Bethe Salpeter Equation solver
bse                          # [R BSE] Bethe Salpeter Equation.
bsk                          # [R BSK] Bethe Salpeter Equation kernel
em1d                         # [R Xd] Dynamical Inverse Dielectric Matrix
StdoHash=  40                # [IO] Live-timing Hashes
Nelectro= 20.00000           # Electrons number
ElecTemp= 0.000000     eV    # Electronic Temperature
BoseTemp=-1.000000     eV    # Bosonic Temperature
OccTresh=0.1000E-4           # Occupation treshold (metallic bands)
DBsIOoff= "none"             # [IO] Space-separated list of DB with NO I/O. DB=(DIP,X,HF,COLLs,J,GF,CARRIERs,W,SC,BS,ALL)
DBsFRAGpm= "none"            # [IO] Space-separated list of +DB to be FRAG and -DB NOT to be FRAG. DB=(DIP,X,W,HF,COLLS,K,BS,QINDX,
FFTGvecs= 21149        RL    # [FFT] Plane-waves
#WFbuffIO                    # [IO] Wave-functions buffered I/O
NonPDirs= "none"             # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=1000000                   # [RIM] Number of random q-points in the BZ
RandGvec= 200            RL    # [RIM] Coulomb interaction RS components
#QpgFull                     # [F RIM] Coulomb interaction: Full matrix
% Em1Anys
 0.00     | 0.00     | 0.00     |        # [RIM] X Y Z Static Inverse dielectric matrix
%
IDEm1Ref=0                   # [RIM] Dielectric matrix reference component 1(x)/2(y)/3(z)
CUTGeo= "none"               # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY..
% CUTBox
 0.00     | 0.00     | 0.00     |        # [CUT] [au] Box sides
%
CUTRadius= 0.000000          # [CUT] [au] Sphere/Cylinder radius
CUTCylLen= 0.000000          # [CUT] [au] Cylinder length
#CUTCol_test                 # [CUT] Perform a cutoff test in R-space
Chimod= "Hartree"            # [X] IP/Hartree/ALDA/LRC/BSfxc
BSEmod= "retarded"           # [BSE] resonant/retarded/coupling
BSKmod= "SEX"                # [BSE] IP/Hartree/HF/ALDA/SEX
BSSmod= "d"                  # [BSS] (h)aydock/(d)iagonalization/(i)nversion/(t)ddft`
BSENGexx= 40       Ry    # [BSK] Exchange components
#ALLGexx                     # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 8            Ry    # [BSK] Screened interaction block size
#WehDiag                     # [BSK] diagonal (G-space) the eh interaction
#WehCpl                      # [BSK] eh interaction included also in coupling
KfnQPdb= "E < gw/ndb.QP"              # [EXTQP BSK BSS] Database
KfnQP_N= 1                   # [EXTQP BSK BSS] Interpolation neighbours
% KfnQP_E
 0.000000 | 1.000000 | 1.000000 |        # [EXTQP BSK BSS] E parameters  (c/v) eV|adim|adim
%
KfnQP_Z= ( 1.000000 , 0.000000 )       # [EXTQP BSK BSS] Z factor  (c/v)
KfnQP_Wv_E= 0.000000   eV    # [EXTQP BSK BSS] W Energy reference  (valence)
% KfnQP_Wv
 0.00     | 0.00     | 0.00     |        # [EXTQP BSK BSS] W parameters  (valence) eV| eV|eV^-1
%
KfnQP_Wc_E= 0.000000   eV    # [EXTQP BSK BSS] W Energy reference  (conduction)
% KfnQP_Wc
 0.00     | 0.00     | 0.00     |        # [EXTQP BSK BSS] W parameters  (conduction) eV| eV|eV^-1
%
Gauge= "length"              # [BSE] Gauge (length|velocity)
#MetDamp                     # [BSE] Define \w+=sqrt(\w*(\w+i\eta))
DrudeWBS= ( 0.00     , 0.00     )  eV  # [BSE] Drude plasmon
#Reflectivity                # [BSS] Compute reflectivity at normal incidence
BoseCut=  0.10000            # [BOSE] Finite T Bose function cutoff
% BEnRange
  0.00000 | 10.00000 | eV    # [BSS] Energy range
%
% BDmRange
  0.10000 |  0.10000 | eV    # [BSS] Damping range
%
BEnSteps= 300                # [BSS] Energy steps
% BLongDir
 0.000000 | 0.000000 | 1.000000 |        # [BSS] [cc] Electric Field
%
% BSEBands
   7 |  13 |                 # [BSK] Bands range
%
% BSEEhEny
-1.000000 |-1.000000 | eV    # [BSK] Electron-hole energy range
%
WRbsWF                      # [BSS] Write to disk excitonic the WFs
#BSSPertWidth                # [BSS] Include QPs lifetime in a perturbative way
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) eV|adim|adim
%
XfnQP_Z= ( 1.000000 , 0.000000 )       # [EXTQP Xd] Z factor  (c/v)
XfnQP_Wv_E= 0.000000   eV    # [EXTQP Xd] W Energy reference  (valence)
% XfnQP_Wv
 0.00     | 0.00     | 0.00     |        # [EXTQP Xd] W parameters  (valence) eV| eV|eV^-1
%
XfnQP_Wc_E= 0.000000   eV    # [EXTQP Xd] W Energy reference  (conduction)
% XfnQP_Wc
 0.00     | 0.00     | 0.00     |        # [EXTQP Xd] W parameters  (conduction) eV| eV|eV^-1
%
% QpntsRXp
  1 |  9 |                   # [Xp] Transferred momenta
%
% BndsRnXp
   1 |  30 |                 # [Xp] Polarization function bands
%
NGsBlkXp= 8            Ry    # [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
 0.000000 | 0.000000 | 1.000000 |        # [Xp] [cc] Electric Field
%
PPAPntXp= 27.21138     eV    # [Xp] PPA imaginary energy
XTermKind= "none"            # [X] X terminator ("none","BG" Bruneval-Gonze)
XTermEn= 40.81708      eV    # [X] X terminator energy (only for kind="BG")

So I want to know why my results have such large energy difference compaired to experimental results? How can I solve this problem? Thank you.

Best regards,
Tianshu Li
Jilin University

[Daniele Varsano]

Dear Tianshu,
the workflow you use to calculate the spectra I would say it is right.
The fact you do not reproduce experimental excitations can have many causes you need to investigate.
Just some suggestions for a few checks you can do:
1) quality of the ground state DFT calculations in terms of approximations (xc functional)
2) Beside the gap a careful check of the convergences of the QP energies for other states (in particular the ones building up the excitation you want to insepct)
3) Comparison of the QP energies, not only the gap if you have experimental data available
4) Finally, if I have well understood you are dealing with a monolayer right? In this case, you may want to check the convergence with respect the used vacuum and
eventually activate a coulomb cutoff technique in order to avoid sporious interactions among celll replica (CUTGeo= "box (dir)" ). Where "dir" stands for x,y or z , ie the direction
orthogonal to your monolayer (I presume x or y) as you set the electric field along z.

..there could be many other reason, but this is not a problem with the code.

Best,
Daniele

[Tianshu]

Dear Daniele,

Thank you for such detailed suggestions. I still have one question.

In this case, you may want to check the convergence with respect the used vacuum and
eventually activate a coulomb cutoff technique...

How do I know the used vacuum is converged? Is this also checked by QP energies?

Thank you for your help and patience.

Best regards,
Tianshu Li
Jilin University

[Daniele Varsano]

Dear Tianshu,
the convergence of QP energies (gap) with respect the vacuum turns to be quite slow.
But also the binding energy is affected, usually in the opposite direction, compensating the GW gap dependence.
A full GW+BSE calculation would be better. When considering a larger value of the vacuum to be consistent with the previous calculations you would need
a larger number of G vectors (if you assign it by energy cutoff as you already did, this is automatically taken into account by Yambo) and a larger number of bands
in the sum over states (e.g. check the energy of the larger conduction band you include).

In any case, usually for non 3D periodic system the cutoff technique it is used (see e.g. Phys. Rev. B 73, 205119, Phys. Rev. B 73, 233103 and many other recipes on the topic)

Best,

Daniele