Strange BSE results

[DavidPolito93]

Dear Developers,

I am studying 1D systems (GW + BSE), however I found strange results for the excitonic binding energy for which I hope you could point me what kind of errors might be possibly arising in my use of Yambo.

First of all I have performed self consistent GW on 400 kpoints (which are reduced to 201 by symmetry) and I have checked the convergence of the results. The final database is now called g2w2

Next step I copied the results of the self consistent GW inside the SAVE folder, so that the screening is not recomputed (i.e. cp -r g2w2/* ./SAVE/) during the BSE.

Then I performed BSE computation with the following input:

CASE_IN=BSE.in
CASE_FL=BSE

cat > $WORKDIR/$CASE_IN << EOF
rim_cut
optics # [R] Linear Response optical properties
bss # [R] BSE solver
bse # [R][BSE] Bethe Salpeter Equation.
dipoles # [R] Oscillator strenghts (or dipoles)
em1d # [R][X] Dynamically Screened Interaction
ppa # [R][Xp] Plasmon Pole Approximation for the Screened Interaction
WFbuffIO # [IO] Wave-functions buffered I/O
StdoHash= 80 # [IO] Live-timing Hashes
BS_nCPU_LinAlg_INV= $NCPUS # [PARALLEL] CPUs for Linear Algebra
BS_nCPU_LinAlg_DIAGO= $NCPUS # [PARALLEL] CPUs for Linear Algebra
K_Threads=0 # [OPENMP/BSK] Number of threads for response functions
NonPDirs= "none" # [X/BSS] Non periodic chartesian directions (X,Y,Z,XY...)
RandQpts=0 # [RIM] Number of random q-points in the BZ
RandGvec= 1 RL # [RIM] Coulomb interaction RS components
CUTGeo= "ws Z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws X/Y/Z/XY..
CUTwsGvec= 1.1000 # [CUT] WS cutoff: number of G to be modified
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX/BSfxc
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(s)lepc/(i)nversion/(t)ddft
BSENGexx= 30 Ry # [BSK] Exchange components
BSENGBlk= 10 Ry # [BSK] Screened interaction block size [if -1 uses all the G-vectors of W(q,G,Gp)]
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "E < SAVE/ndb.QP" # [EXTQP BSK BSS] Database action
KfnQP_INTERP_NN= 1 # [EXTQP BSK BSS] Interpolation neighbours (NN mode)
KfnQP_INTERP_shells= 20.00000 # [EXTQP BSK BSS] Interpolation shells (BOLTZ mode)
KfnQP_DbGd_INTERP_mode= "NN" # [EXTQP BSK BSS] Interpolation DbGd mode
% 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.000000 | 0.000000 | 0.000000 | # [EXTQP BSK BSS] W parameters (valence) eV| 1|eV^-1
%
KfnQP_Wv_dos= 0.000000 eV # [EXTQP BSK BSS] W dos pre-factor (valence)
KfnQP_Wc_E= 0.000000 eV # [EXTQP BSK BSS] W Energy reference (conduction)
% KfnQP_Wc
0.000000 | 0.000000 | 0.000000 | # [EXTQP BSK BSS] W parameters (conduction) eV| 1 |eV^-1
%
KfnQP_Wc_dos= 0.000000 eV # [EXTQP BSK BSS] W dos pre-factor (conduction)
BSEprop= "abs" # [BSS] abs/kerr/magn/dichr trace
% BSEQptR
1 | 1 | # [BSK] Transferred momenta range
%
% BSEBands
3 | 8 | # [BSK] Bands range
%
% BEnRange
0.00000 | 5.00000 | eV # [BSS] Energy range
%
% BDmRange
0.100000 | 0.100000 | eV # [BSS] Damping range
%
BEnSteps= 500 # [BSS] Energy steps
% BLongDir
0.000000 | 0.000000 | 1.000000 | # [BSS] [cc] Electric Field
%
#WRbsWF # [BSS] Write to disk excitonic the WFs
XfnQPdb= "E < SAVE/ndb.QP" # [EXTQP Xd] Database action
XfnQP_INTERP_NN= 1 # [EXTQP Xd] Interpolation neighbours (NN mode)
XfnQP_INTERP_shells= 20.00000 # [EXTQP Xd] Interpolation shells (BOLTZ mode)
XfnQP_DbGd_INTERP_mode= "NN" # [EXTQP Xd] Interpolation DbGd mode
% 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.000000 | 0.000000 | 0.000000 | # [EXTQP Xd] W parameters (valence) eV| 1|eV^-1
%
XfnQP_Wv_dos= 0.000000 eV # [EXTQP Xd] W dos pre-factor (valence)
XfnQP_Wc_E= 0.000000 eV # [EXTQP Xd] W Energy reference (conduction)
% XfnQP_Wc
0.000000 | 0.000000 | 0.000000 | # [EXTQP Xd] W parameters (conduction) eV| 1 |eV^-1
%
XfnQP_Wc_dos= 0.000000 eV # [EXTQP Xd] W dos pre-factor (conduction)
% BndsRnXp
1 | 40 | # [Xp] Polarization function bands
%
NGsBlkXp= 969 RL # [Xp] Response block size
% LongDrXp
0.00000 | 0.00000 | 0.100000E-4 # [Xp] [cc] Electric Field
%
PPAPntXp= 54.42276 eV # [Xp] PPA imaginary energy
XTermKind= "BG" # [X] X terminator ("none","BG" Bruneval-Gonze)
GTermKind = "BG"
EOF

mpirun -np $NCPUS /home/romanind/yambo-5.0.0/bin/yambo -F $CASE_IN -J $CASE_FL

My problem is that the difference between the first singlet peak and the QP peak in the alpha file that I compute is ~1.3 eV, while other works find a ~0.4 eV binding energy.

Moreover, if I compute the excitonic wf it extends up to ~20 unit cells which is quite strange considering such a high binding energy in my case.

So I was wondering where could be the error. My various ideas are:

1) cutoff: my chain is disposed along the z axis and I have a orthorhombic cell (ibrav = 8 in pw.x of QE). Therefore I set "ws Z" in the cutoff section and no random points for the integration as I have seen in that they are not necessary for this kind of cutoff, as I have done for the GW computation. Is there something wrong with the cutoff configuration in your opinion?

2) screening: there might be something wrong with the screening. I have given the one computed via evGW by copying the results into the SAVE folder and indeed in the LOG file I have seen that Yambo recognizes it and skips the evaluation of the dipoles and the various X0q. Is it better to recompute it in a static approximation?

3) convergence: I tried to increase the number of k-points and the various parameters in the BSE kernel but the first peak is not moving so I am assuming that I have reached convergence.

I am really wondering where such a high difference in the binding energy might be arising.

Davide Romanin
-----------------------------------------------------
Post-doctoral fellow
Institut des NanoSciences de Paris
Sorbonne Université, CNRS
4 place Jussieu,
75252 PARIS Cedex 05

[Daniele Varsano]

Dear Davide,
I can't see nothing wrong in your input.
The problem in my opinion is that you are using for the BSE a screening calculated using QP eigenvalues, and in this way the screening is reduced. There is a paper discussing this issue, but at the moment I can't find it. Maybe someone else can add it to this post. My suggestion is to recompute the screening using KS eigenvalues.

Best,
Daniele

[DavidPolito93]

Dear Daniele,

Thank you for your answer!

So I just recompute the static screening right?

Davide

[Daniele Varsano]

Yes,
avoiding to use the ndb.QP database.

Daniele

[DavidPolito93]

Dear Daniele,

I did again the BSE computation with a static screening.
Indeed the binding energy has been reduced from 1.3 to 1.0 eV (i.e. by ~0.3 eV).

However there is still a huge discrepancy wrt other theoretical works (Eb ~0.5 eV).
I was wondering if it might be related to some errors I might have introduced in defining the low dimensional system:
I have used a orthorhombic cell (ibrav = 8) and therefore I have employed the ws z cutoff (as the chain is along z). As far as I have understood there is no need for the RIM with this kind of cutoff.

The only thing I didn't check is the vacuum size: right now I am using ~10Å along x and y, do you think that increasing vacuum in non-periodic directions could give rise to a reduction in the singlet binding energy even if there is the ws cutoff active?

Best,

Davide Romanin

[Daniele Varsano]

Dear Davide,

yes, you can check the effect of the vacuum, and also the k point sampling as the accuracy of the ws cutoff depends on the sampling.
You can also check the CUTwsGvec increasing it, even if the default value should be ok.

Best,
Daniele

[DavidPolito93]

Ok, I will check better the vacuum.

Concerning the k-sampling I have used meshes of 1x1x100, 1x1x200, 1x1x300 and 1x1x400 but the position of the first excitonic peak in the alpha hasn't changed.
Do you think that increasing the x and y k meshes can be important even if I have a 1D system?

Davide

[Daniele Varsano]

Dear Davide,

Do you think that increasing the x and y k meshes can be important even if I have a 1D system?

No, if you want to isolate the system the sampling should be in one direction only.
Best.
Daniele