GW tutorial. Convergence and approximations (BN)

We chosen hexagonal boron nitride to explain the use of yambopy. Along this tutorial we show how to use yambopy to make efficient convergence tests, to compare different approximations and to analyze the results.

The initial step is the ground state calculation and the non self-consistent calculation using the gs_bn.py file:

<source lang="bash">python gs_bn.py python gs_bn.py -sn</source> We have set 50 bands and the k-grid 12x12x1.

GW convergence

(a) Calculations We converge the main parameters of a GW calculation independently. We make use of the plasmon pole approximation for the dielectric function and the Newton solver to find the GW correction to the LDA eigenvalues. The quantity to converge is the band gap of the BN (conduction and valence band at the K point of the Brillouin zone). We can select this calculation by calling the YamboIn with the right arguments:

<source lang="python">y = YamboIn('yambo -d -g n -p p -V all',folder='gw_conv')</source> The main variables are:

EXXRLvcs: Exchange self-energy cutoff

BndsRnXp: Number of bands in the calculation of the dielectric function (PPA).

NGsBlkXp: Cutoff of the dielectric function.

GbndRnge: Self-energy. Number of bands.

The convergence with the k-grid is done after these variables are converged and in principle is also independent of them. The convergence is set with a dictionary in which we choose the parameter and the values. Be aware of setting the right units and format for each parameter.

<source lang="python">conv = { 'EXXRLvcs': [[1,20,40,60,80,100],'Ry'],

        'NGsBlkXp': [[0,1,2,3,4], 'Ry'],
        'BndsRnXp': [[[1,5],[1,10],[1,20],[1,30],[1,40]]],] ,
        'GbndRnge': [[[1,5],[1,10],[1,20],[1,30],[1,40]],] }</source>

The class YamboIn includes the function optimize, which is called here:

<source lang="python">y.optimize(conv,run=run,ref_run=False)</source> This optimization function just needs the convergence dictionary and the run instructions, given by the function:

<source lang="python">def run(filename):

   """ Function to be called by the optimize function """
   folder = filename.split('.')[0]
   print(filename,folder)
   shell = bash() 
   shell.add_command('cd gw_conv; %s -F %s -J %s -C %s 2> %s.log'%(yambo,filename,folder,folder,folder))
   shell.run()
   shell.clean()</source>

We set an interactive run, in the folder gw_conv. All the calculations will be made there with the corresponding jobname.

(b) Analysis

Once all the calculations are finished it's time to pack all the files in the json format for posterior analysis. For this use the YamboOut() class:

<source lang="python">#pack the files in .json files for dirpath,dirnames,filenames in os.walk('gw_conv'):

 #check if there are some output files in the folder
 if ([ f for f in filenames if 'o-' in f ]):
     y = YamboOut(dirpath,save_folder=dirpath)
     y.pack()</source>

This snippet of code can be called using the function:

<source lang="python">pack_files_in_folder('gw_conv',save_folder='gw_conv')</source> Yambopy provides the function yambopy analysegw to perform the analysis of the json files in an automatic way. By running the script selecting the bands and kpoints, together with the parameter we will obtain the convergence plot.

<source lang="python">yambopy analysegw -bc 5 -kc 19 -bv 4 -kv 19 gw_conv EXXRLvcs yambopy analysegw -bc 5 -kc 19 -bv 4 -kv 19 gw_conv NGsBlkXp yambopy analysegw -bc 5 -kc 19 -bv 4 -kv 19 gw_conv BndsRnXp yambopy analysegw -bc 5 -kc 19 -bv 4 -kv 19 gw_conv GbndRnge</source>

From the convergence plot we can choose now a set of parameters and repeat the calculation for finer k-grids until we reach convergence with the k-points. The convergence criteria are left to the user.

GW calculation in a regular grid and plot in a path in the Brillouin zone

We will work in the PPA for the screening. We have chosen the following parameters:

<source lang="bash">FFTGvecs = 20 Ha BndsRnXp = 24 bands NGsBlkXp = 500 mHa GbndRnge = 20 bands EXXRLvcs = 20 Ha QPkrange = [1,19,2,6]</source> We can just simply run the code to calculate the GW corrections for all the points of the Brillouin zone by setting the convergence parameters in the function gw of the script and doing:

<source lang="bash">python gw_conv_bn.py -g</source> The first image show all the GW energies along all the k-points of the Brillouin zone. A clearer picture can be obtained by plotting the band structure along the symmetry points GMKG by using the analyser:

<source lang="bash">python gw_conv_bn.py -r</source> We first pack the results in a json file and subsequently we use the analyser to create the object which contains all the information.

<source lang="bash">pack_files_in_folder('gw') ya = YamboAnalyser('gw')</source> The object ya contains all the results written in the output. We can plot any output variable. In yambopy we provide a function to plot the band structure along a given path. The BN band structure is shown below. The GW correction opens the LDA bandgap as expected.

Approximations of the dielectric function (COHSEX, PPA, Real axis integration)

We can use yambopy to examine different run levels. For instance, the approximations used to obtain the screening are the: (i) static screening or COHSEX, plasmon-pole approximations (PPA), or real axis integration. We have set the same parameters for each run, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">COHSEX BndsRnXs = 24 bands NGsBlkXs = 500 mHa PPA BndsRnXp = 24 bands NGsBlkXp = 500 mHa RA BndsRnXd = 24 bands NGsBlkXd = 500 mHa</source> We have set the converged parameters and the function works by running:

<source lang="bash">python gw_conv_bn.py -x</source> We plot the band structure using the analyzer explained above.

<source lang="bash">python gw_conv_bn.py -xp</source> The PPA and the RA results are basically on top of each other. On the contrary, the COHSEX (static screening) makes a poor job, overestimating the bandgap correction.

Solvers (Newton, Secant, Green's function)

The solvers to find the QP correction from the self-energy can also be tested. We have included the Newton and the secant method. In the resulting band structures we do not appreciate big differences. In anycase it is worthy to test during the convergence procedure.