Yambo Community Forum

Dear all,

I am using Yambo 4.5.1 on a 1D system with a WS cutoff. I am trying to compute the excitonic corrections. However the program stops at the kernel step.

First of all I have computed G0W0 corrections to the eigenvalues in the all_BZ folder:

HF_and_locXC # [R XX] Hartree-Fock Self-energy and Vxc
ppa # [R Xp] Plasmon Pole Approximation
gw0 # [R GW] GoWo Quasiparticle energy levels
rim_cut # [R RIM CUT] Coulomb potential
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
DIP_CPU="1 $ncpu 1"
DIP_ROLEs="k c v"
DIP_Threads=0
X_CPU="1 1 1 $ncpu 1"
X_ROLEs="q q k c v"
X_nCPU_LinAlg_INV= $ncpu
X_Threads=0 # [OPENMP/X] Number of threads for response functions
SE_CPU="1 $ncpu 1"
SE_ROLEs="q qp b"
SE_Threads=0 # [OPENMP/GW] Number of threads for self-energy
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
EXXRLvcs= 50 Ry # [XX] Exchange RL components
VXCRLvcs= 424401 RL # [XC] XCpotential RL components
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
% BndsRnXp
1 | 30 | # [Xp] Polarization function bands
%
NGsBlkXp= 3 Ry # [Xp] Response block size
XTermKind = "BG"
% LongDrXp
1.000000 | 1.000000 | 1.000000 | # [Xp] [cc] Electric Field
%
PPAPntXp= 40.00 eV # [Xp] PPA imaginary energy
% GbndRnge
1 | 100 | # [GW] G[W] bands range
%
GTermKind = "BG"
GDamping= 0.10000 eV # [GW] G[W] damping
dScStep= 0.10000 eV # [GW] Energy step to evaluate Z factors
DysSolver= "n" # [GW] Dyson Equation solver ("n","s","g")
%QPkrange # # [GW] QP generalized Kpoint/Band indices
1|51|3|6|
%

Then I wanted do compute the BSE spectrum:

~/yambo-4.5.1/bin/yambo -J all_BZ -F test.in -r -o b -p p -y d -k sex -V all

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.
dipoles # [R ] Compute the dipoles
em1d # [R Xd] Dynamical Inverse Dielectric Matrix
StdoHash= 40 # [IO] Live-timing Hashes
Nelectro= 8.000000 # Electrons number
ElecTemp= 0.000000 eV # Electronic Temperature
BoseTemp= 0.000000 eV # Bosonic Temperature
OccTresh=0.1000E-4 # Occupation treshold (metallic bands)
NLogCPUs=0 # [PARALLEL] Live-timing CPUs (0 for all)
DBsIOoff= "none" # [IO] Space-separated list of DB with NO I/O. DB=(DIP,X,HF,COLLs,J,GF,CARRIERs,OBS,W,SC,BS,ALL)
DBsFRAGpm= "none" # [IO] Space-separated list of +DB to FRAG and -DB to NOT FRAG. DB=(DIP,X,W,HF,COLLS,K,BS,QINDX,RT,ELP
FFTGvecs= 18733 RL # [FFT] Plane-waves
#WFbuffIO # [IO] Wave-functions buffered I/O
PAR_def_mode= "balanced" # [PARALLEL] Default distribution mode ("balanced"/"memory"/"workload")
DIP_CPU= "" # [PARALLEL] CPUs for each role
DIP_ROLEs= "" # [PARALLEL] CPUs roles (k,c,v)
DIP_Threads=0 # [OPENMP/X] Number of threads for dipoles
X_CPU= "" # [PARALLEL] CPUs for each role
X_ROLEs= "" # [PARALLEL] CPUs roles (q,g,k,c,v)
X_nCPU_LinAlg_INV= 1 # [PARALLEL] CPUs for Linear Algebra
X_Threads=0 # [OPENMP/X] Number of threads for response functions
BS_CPU= "" # [PARALLEL] CPUs for each role
BS_ROLEs= "" # [PARALLEL] CPUs roles (k,eh,t)
BS_nCPU_LinAlg_INV= 1 # [PARALLEL] CPUs for Linear Algebra
BS_nCPU_LinAlg_DIAGO= 1 # [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
#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= "ws z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere/ws 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
CUTwsGvec= 1.100000 # [CUT] WS cutoff: number of G to be modified
#CUTCol_test # [CUT] Perform a cutoff test in R-space
Chimod= "HARTREE" # [X] IP/Hartree/ALDA/LRC/PF/BSfxc
ChiLinAlgMod= "LIN_SYS" # [X] inversion/lin_sys
BSEmod= "retarded" # [BSE] resonant/retarded/coupling
BSKmod= "SEX" # [BSE] IP/Hartree/HF/ALDA/SEX
Lkind= "Lbar" # [BSE] Lbar (default) / full
BSSmod= "d" # [BSS] (h)aydock/(d)iagonalization/(i)nversion/(t)ddft
% DipBands
1 | 30 | # [DIP] Bands range for dipoles
%
#DipBandsALL # [DIP] Compute all bands range, not only valence and conduction
DipApproach= "G-space v" # [DIP] [G-space v/R-space x/Covariant/Shifted grids]
#DipPDirect # [DIP] Directly compute <v> also when using other approaches for dipoles
ShiftedPaths= "" # [DIP] Shifted grids paths (separated by a space)
DbGdQsize= 1.000000 # [X,DbGd][o/o] Percentual of the total DbGd transitions to be used
BSENGexx= 424401 RL # [BSK] Exchange components
#ALLGexx # [BSS] Force the use use all RL vectors for the exchange part
BSENGBlk= 4 Ry # [BSK] Screened interaction block size
#WehDiag # [BSK] diagonal (G-space) the eh interaction
#WehCpl # [BSK] eh interaction included also in coupling
KfnQPdb= "none" # [EXTQP BSK BSS] Database
KfnQP_N= 1 # [EXTQP BSK BSS] Interpolation neighbours
% KfnQP_E
1.51151 | 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| 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.00 | 0.00 | 0.00 | # [EXTQP BSK BSS] W parameters (conduction) eV| 1 |eV^-1
%
KfnQP_Wc_dos= 0.000000 eV # [EXTQP BSK BSS] W dos pre-factor (conduction)
Gauge= "length" # [BSE] Gauge (length|velocity)
#NoCondSumRule # [BSE] Do not impose the conductivity sum rule in velocity gauge
#MetDamp # [BSE] Define
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
% BSEBands
3 | 6 | # [BSK] Bands range
%
% BSEEhEny
-1.000000 |-1.000000 | eV # [BSK] Electron-hole energy range
%
% BEnRange
0.00000 | 10.00000 | eV # [BSS] Energy range
%
% BDmRange
0.10000 | 0.10000 | eV # [BSS] Damping range
%
BDmERef= 0.000000 eV # [BSS] Damping energy reference
BEnSteps= 500 # [BSS] Energy steps
% BLongDir
1.000000 | 0.000000 | 0.000000 | # [BSS] [cc] Electric Field
%
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| 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.00 | 0.00 | 0.00 | # [EXTQP Xd] W parameters (conduction) eV| 1 |eV^-1
%
XfnQP_Wc_dos= 0.000000 eV # [EXTQP Xd] W dos pre-factor (conduction)
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| 1 |eV^-1
%
XfnQP_Wc_dos= 0.000000 eV # [EXTQP Xd] W dos pre-factor (conduction)
% QpntsRXp
1 | 51 | # [Xp] Transferred momenta
%
% BndsRnXp
1 | 30 | # [Xp] Polarization function bands
%
NGsBlkXp= 179 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
0.5774E-5 |0.5774E-5 |0.5774E-5 | # [Xp] [cc] Electric Field
%
PPAPntXp= 40.00000 eV # [Xp] PPA imaginary energy
XTermKind= "none" # [X] X terminator ("none","BG" Bruneval-Gonze)
XTermEn= 40.00000 eV # [X] X terminator energy (only for kind="BG")

However, the execution of Yambo stops at this point:

[07.03] BSE Kernel @q1 (Resonant CORRRELATION EXCHANGE)
=======================================================

[BSE] Exchange components : 424401

[WARNING] Exchange Kernel FFT size is too big. NG_X reduced: 424401 --> 20825

[WARNING] Bigger FFT discarded to avoid slow computation for corr part of Kernel

[07.03.01] Screened interaction header I/O
==========================================

[RD./all_BZ//ndb.pp]----------------------------------------
Brillouin Zone Q/K grids (IBZ/BZ): 51 100 51 100
RL vectors (WF): 18733
Coulomb cutoff potential :ws z 1.100
Fragmentation :yes
Electronic Temperature [K]: 0.000000
Bosonic Temperature [K]: 0.000000
PPA diel. fun. energies :Perdew, Burke & Ernzerhof(X)+Perdew, Burke & Ernzerhof(C)
wavefunctions :Perdew, Burke & Ernzerhof(X)+Perdew, Burke & Ernzerhof(C)
Global Gauge :length
*ERR* X matrix size : 179
X band range : 1 30
X e/h energy range [ev]:-1.000000 -1.000000
X Time ordering :R
X xc-Kernel :none
X Drude frequency : 0.00 0.00
X poles [o/o]: 100.0000
RL vectors in the sum : 18733
[r,Vnl] included :yes
Field direction :0.5774E-5 0.5774E-5 0.5774E-5
BZ energy Double Grid :no
BZ energy DbGd points :0
BZ Q point size factor : 1.000000

[WARNING]Variable not compatible in PP/Em1s DB

[WARNING]Bethe Salpter section skipped. Impossible to build the kernel.

Can you help me identifying the problem?

Sincerely,
Davide Romanin
-----------------------------------------------------
PhD student in Physics XXXIII cycle
Representative of the PhD students in Physics
Applied Science and Technology department (DiSAT)
Politecnico di Torino
Corso Duca degli Abruzzi, 24
10129 Torino ITALY
------------------------------------------------------

Dear Davide,

you are asking in input 4y for the screening:

BSENGBlk= 4 Ry

but you calculated the screening using 179 Gblk that most probably correspond t 3Ry according with your GW input file: so you are asking for more blocks than the ones that have been calculated. So either you as for 3Ry in BSENGBlk or you calculate the screening using more blocks NGsBlkXp= 4 Ry

Not related with the error message:
note you calculated the GW correction for 4 bands in he whole BZ, but than you are using a scissor operator. Probably you already checked this is a good approximation anyway you can read the calculated QP energies by setting in input: if the database is in the SAVE directory, or substitute it with the correct path.

Best,
Daniele

Dear Daniele,

Thank you a lot for your reply!

I see, my mistake! I was using the converged values on the screening and I pushed too much on the BSE kernel.

For what concerns the scissor, yes I used it just for a test but I am going to read the QP correction for the final result. Thank you

I just have a final question, I also tried using the ypp executable for analyzing the excitons and I wanted to plot their spatial distribution like in the tutorial. So I did:

~/yambo-4.5.1/bin/ypp -F ypp_WF.in -J all_BZ -e w 1

and in the input I specified:

excitons # [R] Excitons
wavefunction # [R] Wavefunction
Format= "x" # Output format [(c)ube/(g)nuplot/(x)crysden]
Direction= "3" # [rlu] [1/2/3] for 1d or [12/13/23] for 2d [123] for 3D
FFTGvecs= 18733 RL # [FFT] Plane-waves
States= "4 - 4" # Index of the BS state(s)
Qpts= 1 # Index q of the BS state(s)
Degen_Step= 0.0100 eV # Maximum energy separation of two degenerate states
% Cells
1 | 1 | 2 | # Number of cell repetitions in each direction (odd or 1)
%
% Hole
0.00 | 0.00 | 0.00 | # [cc] Hole position in unit cell (positive)
%

However, when I run it the output format is switched from xcrysden to gnu plot:

<---> [01] CPU structure, Files & I/O Directories
<---> CPU-Threads:1(CPU)-6(threads)
<---> [02] Y(ambo) P(ost)/(re) P(rocessor)
<---> [03] Core DB
<---> :: Electrons : 8.000000
<---> :: Temperature [ev]: 0.000000
<---> :: Lattice factors [a.u.]: 18.89727 18.89727 4.78857
<---> :: K points : 51
<---> :: Bands : 100
<---> :: Symmetries : 16
<---> :: RL vectors : 424401
<---> [04] K-point grid
<---> :: Q-points (IBZ): 51
<---> :: X K-points (IBZ): 51
<---> [05] CORE Variables Setup
<---> [05.01] Unit cells
<---> [05.02] Symmetries
<---> [05.03] RL shells
<---> [05.04] K-grid lattice
<---> Grid dimensions : 100
<---> [05.05] Energies [ev] & Occupations
<---> [06] Excitonic Properties
<---> :: Sorting energies
<---> 1 excitonic states selected
<---> [06.01] Excitonic Wave Function
<---> [WARNING] Switching to GnuPlot
<---> [06.01.01] Real-Space grid setup
<01s> [WF-EXCWF] Performing Wave-Functions I/O from ./SAVE
<01s> [FFT-EXCWF] Mesh size: 54 54 14
<05s> :: Extended grid : 54 54 42
<05s> :: Hole position in the DL cell [cc]: 0.00 0.00 0.00
<05s> :: position in the FFT grid [cc]: 0.00 0.00 0.00
<05s> :: translated position [cc]: 0.000000 0.000000 4.788568
<05s> :: [A]: 0.000000 0.000000 2.534001
<05s> Processing 1 states
<06s> ExcWF@4 |########################################| [100%] --(E) --(X)
<06s> 3D Merge |########################################| [100%] --(E) --(X)
<06s> [07] Game Over & Game summary

Am I making a mistake in the input or is it supposed to happen in a 1D system?

Sincerely,

Davide Romanin
-----------------------------------------------------
PhD student in Physics XXXIII cycle
Representative of the PhD students in Physics
Applied Science and Technology department (DiSAT)
Politecnico di Torino
Corso Duca degli Abruzzi, 24
10129 Torino ITALY
------------------------------------------------------

Dear Davide,
that's correct, you indicate in input to have a plot in 1D (i.e. integrated in the x/y directions), xcrysden it is meant for isosurfaces or contour plot (3D and 2D).
Even if you have 1D system, it lives in a 3D space and you can plot an isosurface, so you can indicate:
Direction= "123"

Best,
Daniele

Dear all,

I am still having some troubles with the post processing of excitons.

More precisely, I have sorted the excitons according to:

ypp -J 2D_WR_WC -e s 1

And the output of the first 12 excitons is the following (o-all_BZ.exc_qpt1_E_sorted):

# E [ev] Strength Index
#
0.820348442 0.228709985E-10 1.00000000
0.9082934260 0.3068218346E-8 2.000000000
0.9107534885 0.1027362551E-7 3.000000000
1.206864833832 1.000000000000 4.000000000000
1.512270331 0.2654584996E-8 5.000000000
1.51745868 0.215334001E-09 6.00000000
1.51866841 0.671150982E-10 7.00000000
1.57387865 0.967018132E-09 8.00000000
1.80781424 0.634213645E-11 9.00000000
1.82887924 0.370891790E-09 10.0000000
1.82920873 0.596988153E-10 11.0000000
1.90853464603 0.01866552792 12.00000000000

However, when I try to compute the oscillator strength of the 4th exciton (which has the strongest strength) with:

excitons # [R] Excitons
amplitude # [R] Amplitude
States= "4 - 4" # Index of the BS state(s)
Qpts= 1 # Index q of the BS state(s)
#DipWeight # Weight the contribution to the exciton WFs with the dipoles
Degen_Step= 0.0100 eV # Maximum energy separation of two degenerate states

I obtain an empty output:

# Electron-Hole pairs that contributes to Excitonic State 4 for iq=1 more than 5.000000%
#
# K-point [iku] Weight
# 0.00000 0.00000 0.47250 0.05436
# 0.00000 0.00000 0.47500 0.08326
# 0.000000 0.000000 0.477500 0.124532
# 0.000000 0.000000 0.480000 0.181605
# 0.000000 0.000000 0.482500 0.258419
# 0.000000 0.000000 0.485000 0.357712
# 0.000000 0.000000 0.487500 0.480477
# 0.000000 0.000000 0.490000 0.622146
# 0.000000 0.000000 0.492500 0.770633
# 0.000000 0.000000 0.495000 0.905365
# 0.000000 0.000000 0.497500 1.000000
# 0.000000 0.000000 -0.500000 0.517285
#
# Band_V Band_C Kv-q ibz Symm_kv Kc ibz Symm_kc Weight Energy
#
#
# 12/29/2020 at 16:52 YPP @ r033c01s04 [start]
# 12/29/2020 at 16:52 [end]

I also tried with a range States= "1 - 12" # Index of the BS state(s). Some of them are not empty while others are empty.

Am I making some mistakes?

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

Web site: https://dromanincm.github.io
------------------------------------------------------

Dear Davide,
as reported in the output : by default, yambo prints only the transition that participates with a weight higher than 5%.
In cases where none of the transitions has a weight higher that 5% the output is empty.
You can change the threshold in the input file by using the variable Weight_treshold, if you set: all the transitions with a weight higher than 1% will be reported in the output, and if you set it to 0, all the transitions will be included.

The variable appears in the input if you add verbosity when building the input file: Best,
Daniele