How to treat low dimensional systems
In this tutorial you will learn how to treat lowdimensional material. A small ksampling is used to reduce the computational cost/time. Later on you can repeat this tutorial but using denser kgrids to check the convergence.
Contents
 1 Prerequisites
 2 Avoid numerical divergences using the Random Integration Method (RIM)
 3 Generate a truncated coulomb potential
 4 Use the truncated coulomb potential in a G_{0}W_{0}PPA calculation
 5 Use the truncated coulomb potential in a BSE calculation
 6 Analyze the differences with corresponding GW and BSE calculations without the use of a truncated coulomb potential
Prerequisites

yambo
executable 
ypp
executable 
gnuplot
, for plotting spectra  Have Completed the First steps: a walk through from DFT to optical properties tutorial
 Have Completed the How to obtain an optical spectrum  a guide through the workflow of calculations tutorial
Avoid numerical divergences using the Random Integration Method (RIM)
In DFT runs of lowdimensional materials low dimensional kgrids are generally used. (i.e. NxNx1 for a 2D sheet perpendicular to the z direction) This can create numerical problems in the convergence of the manybody (MB) results due to the divergence at small q of the coulomb potential, which appears in all the main equations. To eliminate this problem YAMBO uses the socalled Random Integration Method which means to use a Monte Carlo Integration with Random Qpoints whose number RandQpts will be given in input. For example the exchange selfenergy matrix element which is:
assuming that the integrand is a smooth function of momenta the integral can be approximated as
where the small Brillouin Zones (sBZ) relative to a given qpoint are the Brillouin Zones of the momenta vectors lattice. They are chosen in such a way to cover the whole BZ. In the RIM runlevel yambo calculates the integrals of the symmetrized Coulomb potential
So we are ready to start the tutorial
First go in the proper directory and have a look if all the required databases are there
$ cd YAMBO_TUTORIALS/hBN2D/YAMBO $ ls ./SAVE ns.db1 ns.wf ns.kb_pp_pwscf ... $ yambo (initialization, if you haven't already done so)
Create the input to generate the ndb.RIM
database
$ yambo F yambo_WR.in r
and change the following variables
RandQpts = 1000000 # [RIM] Number of random qpoints in the BZ RandGvec = 100 RL # [RIM] Coulomb interaction RS components
N.B RandGvec=100 means to use the RIM for the first 100 Gcomponents of the coulomb potential (Suggestion : later you can check convergence of the HF gap changing these two values) Close input and Run yambo
$ yambo F yambo_WR.in J 2D_WR
At the end you find a new database ndb.RIM and a new report file r2D_WR_rim_cut . Open it and look inside
Gamma point sphere radius [au]: 0.08028 Points outside the sphere : 799246 [Int_sBZ(q=0) 1/q^2]*(Vol_sBZ)^(1/3) = 7.665674 should be < 7.795600 [WR./2D_WR//ndb.RIM] Brillouin Zone Q/K grids (IBZ/BZ): 7 36 7 36 Coulombian RL components : 111 Coulombian diagonal components :yes RIM random points : 1000000 RIM RL volume [a.u.]: 0.390293 Real RL volume [a.u.]: 0.390112 Eps^1 reference component :0 Eps^1 components : 0.00 0.00 0.00 RIM anysotropy factor : 0.000000  S/N 005962  v.04.01.02 r.00120 
Summary of Coulomb integrals for nonmetallic bands Q[au] RIM/Bare:
Q [1]:0.1000E40.9833 * Q [2]: 0.256404 1.094818 Q [5]: 0.444104 1.032268 * Q [3]: 0.512807 1.024055 Q [6]: 0.678380 1.013997 * Q [4]: 0.769211 1.010774 Q [7]: 0.888208 1.008393
The RIM and Real RL Volumes are quite similar so one million random qpoints seems a reasonable number, but again you are invited, later on, to check the convergence of one of the main observables (i.e. HF gap, GW gap, absorption..) changing this number.
Close the report file and perform a calculation of the exchange selfenergy to estimate the HF gap (very fast calculation) using the RIM. Generate the input
$ yambo F yambo_HF_WR.in J 2D_WR r x
In order to calculate the HF gap only for the last kpoint (which is the highsimmetry point K) change the last line as Note also that 4 (5) is the highest occupied (lowest unoccupied) band.
%[Variables#QPkrangeQPkrange]] # [GW] QP generalized Kpoint/Band indices 7 7 4 5
Close the input and run yambo
$ yambo F yambo_HF_WR.in J 2D_WR
At the end you will find the report file r2D_WR_HF_and_locXC_rim_cut and the output file o2D_WR.hf. Open them and have a look In the output file o2D_WR.hf you will see
# Kpoint Band Eo Eh DFT HF 7.00000 4.00000 0.00000 3.18920 16.21949 19.40869 7.00000 5.00000 4.40109 9.69680 11.10752 5.81181
The HF gap is 12.89 eV obtained subtracting the 2 values in the fourth column Now you can do a similar HF calculation without generating and reading the RIM database
$ yambo F yambo_HF_NR.in x
change the variables QPkrange as before, close the input and run yambo
$ yambo F yambo_HF_NR.in J 2D_NR
in this case the HF gap in the file o2D_NR.hf results to be 12.69 eV.
Indeed the presence of the numerical instability, discussed before, is evident only using denser kgrids with respect to the one used in this Tutorial (6x6x1) Suggestion: later you can generate other SAVE directories with denser kgrids and check the HF gap. The results of the HF gap calculated with different kgrids without (noRIM) and with Random Integration Method (RIM) to show up the problem are reported:
noRIM RIM 6x6x1 12.69eV 12.89eV 12x12x1 12.80eV 12.89eV 15x15x1 12.96eV 12.90eV 45x45x1 15.52eV 12.96eV
So the home message is : use always the RIM in MB simulations of lowdimensional materials.
Generate a truncated coulomb potential
To simulate an isolated nanomaterial a convergence with cell vacuum size is in principle required, like in the DFT runs. The use of a truncated Coulomb potential allows to achieve faster convergence eliminating the interaction between the repeated images along the nonperiodic direction (see i.e. C. A. Rozzi et al Phys. Rev. B 73, 205119) In this tutorial we learn how to generate a boxlike cutoff for a 2D system with the nonperiodic direction along z.
In YAMBO you can use :
 spherical cutoff (for 0D systems)
 cylindrical cutoff (for 1D systems)
 boxlike cutoff (for 0D, 1D and 2D systems)
The Coulomb potential with a boxlike cutoff is defined as
where S is a boxshaped region of the space having sides Lx,Ly and Lz, where they can also assume infinite values (untruncated interaction along desired directions).
Then the FT component is
whereFor a 2Dsystem with non period direction along zaxis we have
Important remarks:
 the Random Integration Method (RIM) is required to perform the Qspace integration
 for sufficiently large supercells a choose L_i slightly smaller than the cell size in the idirection ensures to avoid interaction between replicas
Create the input file:
$ yambo F yambo_WR_WC.in r
Change the RIM Variables as before and activate the cutoff generation
RandQpts = 1000000 # [RIM] Number of random qpoints in the BZ RandGvec = 100 RL # [RIM] Coulomb interaction RS components CUTGeo = "box z" # [CUT] Coulomb Cutoff geometry: box/cylinder/sphere X/Y/Z/XY.. % CUTBox 0.00  0.00  32.0  # [CUT] [au] Box sides
Close the input file and run yambo:
$ yambo F yambo_WR_WC.in J 2D_WR_WC
in the directory 2D_WR_WC you will find a new database ndb.cutoff
Use the truncated coulomb potential in a G_{0}W_{0}PPA calculation
Generate the input reading the databases in the 2D_WR_WC directory
$ yambo J 2D_WR_WC F yambo_G0W0.in p p g n r k hartree V qp In the input change the following variables EXXRLvcs = 40 Ry # [XX] Exchange RL components NGsBlkXp = 4 Ry # [Xp] Response block size QPkrange # [GW] QP generalized Kpoint/Band indices 7 7 4 5
Close the input and run yambo
$ yambo F yambo_G0W0.in J 2D_WR_WC
At the end you will find a new report file r2D_WR_WC_em1d_ppa_HF_and_locXC_gw0_rim_cut, open it and have a look. You will find also a new output file o2D_WR_WC.qp. If you open yoi will find that now the G_{0}W_{0} gap is 4.40 eV (LDA) + 3.68 eV (G0W0 correction) = 8.08 eV Indeed as you have learned in the previous tutorials the correlation part of the selfenergy strongly reduces the HF gap. Moreover you should note that the QP correction is larger that one found in the hBN bulk. Why?
Use the truncated coulomb potential in a BSE calculation
Generate the input for the BSE calculation:
$ yambo J 2D_WR_WC F yambo_BSE.in r o b p p y d k sex V all
Some remarks: the largest verbosity is used V all
and a long input file is generated; with the option p p
the static part of the screening matrix is read from the ndb.pp
databases
Put the number of Gvectors in the exchange and correlation part of the kernel as:
BSENGexx= 40 Ry # [BSK] Exchange components BSENGBlk= 4 Ry # [BSK] Screened interaction block size
use the simple rigid scissor to open the correct the KS energies
% KfnQP_E 3.68000000  1.000000  1.000000  # [EXTQP BSK BSS] E parameters (c/v) eVadimadim
Set the number of bands to:
% BSEBands 2  6  # [BSK] Bands range
This choice means to include in the BSE only three occupied and 2 unoccupied bands. Increase the number of energy steps
BEnSteps= 500 # [BSS] Energy steps
To do the next Tutorial we need to write the excitonic WFs, so uncomment the following line
#WRbsWF # [BSS] Write to disk excitonic the WFs
Close the input file and run yambo
$ yambo J 2D_WR_WC F yambo_BSE.in
Look at the report file r2D_WR_WC_optics_bse_bsk_bss_em1d_ppa_rim_cut and use gnuplot to plot o2D_WR_WC.eps_q1_diago_bse
Analyze the differences with corresponding GW and BSE calculations without the use of a truncated coulomb potential
Open the input file yambo_G0W0.in and set back CUTGeo= "none"
Close the input and run yambo
$ yambo J 2D_WR F yambo_G0W0.in
You will find a new report r2D_WR_optics_bse_bsk_bss_em1d_ppa_rim_cut. You are invited to see the difference with the previous one r2D_WR_WC_optics_bse_bsk_bss_em1d_ppa_rim_cut
and also a new output file o2D_WR.qp is present. Open and check the QP correction at the last kpoint. This time a value of 2.79 eV instead of 3.68 eV is obtained. Are you able to explain this result?
Now redoo the BSE calculation but without reading the cutoff database and applying the QP correction just obtained without using the cutoff.
To do that, open the yambo_BSE.in and put CUTGeo= "none"
and set the QP correction of 2.79 eV, just obtained without uisng the cutoff.
% KfnQP_E 2.790000  1.000000  1.000000  # [EXTQP BSK BSS] E parameters (c/v) eVadimadim
Close the input file and run yambo
$ yambo J 2D_WR F yambo_BSE.in
Plot the dielectric functions
$ gnuplot gnuplot> plot 'o2D_WR.eps_q1_diago_bse' u 1:4 w l title 'GW without cutoff' ,'o2D_WR_WC.eps_q1_diago_bse' u 1:4 w l title ' GW with cutoff' , 'o2D_WR.eps_q1_diago_bse' u 1:2 w l title 'BSE without cutoff' ,'o2D_WR_WC.eps_q1_diago_bse' u 1:2 w l title 'BSE with cutoff'
You can note that the energy of the peaks in the two G0W0 spectra (with and without cutoff) are quite different: you find the first peak around 7 eV without cutoff and around 8 eV with cutoff. Instead the two GW+BSE spectra, without and with cutoff, are much more similar. Are you able to explain why? Remarks: The convergence in kspace can be done separately for GW and BSE calculations, because generally more kpoints are needed to converge the excitonic with respect to the quasiparticle properties.
N.B. If you will use the same inputs created in this tutorial remember to set up the cutoff variable CUTGeo= "box z" to use again the truncated coulomb potential in the GW and BSE calculations
Literature on twodimensional hBN: L. Wirtz et al Phys Rev Lett 96 126104 (2006), F. Rasmussen et al Phys Rev. B 94 155406 (2016),N. Berseneva et al Phys. Rev. B 87 035404 (2013).