View source for GW tutorial. Convergence and approximations (BN)


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

=== Initial steps with Quantum Espresso ===

Yambopy includes qepy, a module to handle QE-DFT calculations necessary to run Yambo.

The initial step consists of a self-consistent (scf) ground state calculation together with a non-self-consistent (nscf) calculation with QE using the <code>gs_bn.py</code> file:

<source lang="bash">python gs_bn.py
python gs_bn.py -sn</source>

We have set the non-self-consistent run with a wave-function cutoff of '''60 Ry, 70 bands and a k-grid of <code>12x12x1</code>''' for Yambo calculations.

In addition, we can check if we get a reasonable band structure for our many-body calculations. We can define a path, for instance, in an hexagonal Brillouin zone:

<source lang="bash">
p = Path([ [[0.0, 0.0, 0.0],'$\Gamma$'],
           [[0.5, 0.0, 0.0],'M'],
           [[1./3,1./3,0.0],'K'],
           [[0.0, 0.0, 0.0],'$\Gamma$']], [int(10*2),int(10),int(sqrt(5)*10)])
</source>

We obtain the band structure by doing:

<source lang="bash">python gs_bn.py -b</source>

[[File:Plot-bn-band.png|BN Band Structure|500px]]

The horizontal line marks the top of the valence band. The electronic bandgap has a value of 4.73 eV. In the next sections, we will show how to calculate the GW correction and the excitonic effects with BSE using Yambo in an automatic way.

=== GW convergence of the bandgap ===

'''(a) Calculations'''

We converge the main parameters of a GW calculation independently. In addition, 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. We converge the band gap of BN (difference in energy of the bottom of the conduction and top of the valence band at the K point of the Brillouin zone). We have designed the script '''gw_conv_bn.py''' (folder ~/tutorial/bn) for this purpose.

We can select the GW calculation by calling the <code>YamboIn</code> with the corresponding arguments:

<source lang="python">y = YamboIn.from_runlevel('yambo -d -p p -g n -V all',folder='gw_conv')</source>

The main variables of a GW calculation are:

<blockquote>

<code>EXXRLvcs</code>: Exchange self-energy cutoff

<code>NGsBlkXp</code>: Cutoff of the dielectric function.

<code>BndsRnXp</code>: Number of bands in the calculation of the dielectric function (PPA).

<code>GbndRnge</code>: Self-energy. The number of bands.

</blockquote>

The convergence in these  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': [[10,10,20,30,40,50,60,70,80,90,100],'Ry'],
         'NGsBlkXp': [[0,0,1,2,3,4,5,6], 'Ry'],
         'BndsRnXp': [[[1,10],[1,10],[1,20],[1,30],[1,40],[1,50],[1,60],[1,70]],''] ,
         'GbndRnge': [[[1,10],[1,10],[1,20],[1,30],[1,40],[1,50],[1,60],[1,70]],''] }</source>
The class <code>YamboIn</code> includes the function <code>optimize</code>, 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&gt; %s.log'%(yambo,filename,folder,folder,folder))
    shell.run()
    shell.clean()</source>

<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')
        shell.add_command('rm -f *.json %s/o-*'%folder) #cleanup
        shell.add_command('%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 <code>gw_conv</code>. All the calculations will be made there with the corresponding jobname.

y.optimize(conv,folder='gw_conv',run=run,ref_run=False)

'''(b) Analysis'''

Once all the calculations are finished it's time to analyse them. Before we do that, let's look at the tools yambopy offers.

Yambopy uses the <code>json</code> format for posterior analysis of the results. To pack all the files in this format,
there is a recipe in yambopy to automatically perform this task on a folder:

<source lang="python">
pack_files_in_folder('gw_conv',save_folder='gw_conv')
</source>

Besides the python module, yambopy can also be called in the terminal to perform some post-analysis tasks:

<source lang="python">
     $ yambopy
     analysebse ->     Using ypp, you can study the convergence of BSE calculations in 2 ways:
       plotem1s ->     Plot em1s calculation
      analysegw ->     Study the convergence of GW calculations by looking at the change in bandgap value.
        mergeqp ->     Merge QP databases
           test ->     Run yambopy tests
   plotexcitons ->     Plot excitons calculation
</source>

Calling <code>yambopy analysegw</code> will display the help of the function:

<source lang="python">

    Study the convergence of GW calculations by looking at the change in bandgap value.

    The script reads from <folder> all results from <variable> calculations and display them.

    Use the band and k-point options according to the size of your k-grid
    and the location of the band extrema.

        Mandatory arguments are:
            folder   -> Folder containing SAVE and convergence runs.
            var      -> Variable tested (e.g. FFTGvecs)

        Optional variables are:
            -bc, --bandc   (int)  -> Lowest conduction band number
            -kc, --kpointc (int)  -> k-point index for conduction band
            -bv, --bandv   (int)  -> Highest valence band number
            -kv, --kpointv (int)  -> k-point index for valence band
            -np, --nopack  (flag) -> Do not call 'pack_files_in_folder'
            -nt, --notext  (flag) -> Do not print a text file
            -nd, --nodraw  (flag) -> Do not draw (plot) the result

</source>

Running the function selecting the bands and k-points, together with the parameter of convergence 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>

By calling <code>python gw_conv_bn.py -p</code> in the terminal, these steps will be performed automatically for this tutorial and you will see the plots below.

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. We have
intentionally used non-converged parameters. Nevertheless, along this week
you should have gotten enough expertise to push the convergence of the parameters
and determine the correct convergence set of parameters.
We invite you to enter in the python script, increase the parameters and check
again the convergence for larger values!


[[File:EXXRLvcs.png|FFTGvecs|300px]]

[[File:BndsRnXp.png|BndsRnXp|300px]]

[[File:NGsBlkXp.png|NGsBlkXp|300px]]

[[File:GbndRnge.png|GbndRnge|300px]]

In general, the convergence with the '''k-grid''' is done after these variables are converged and, in principle, it is also independent of them. We invite you to change the number of k-points in the file '''gs_bn.py''' using the variable '''kpoints_nscf'''.

=== 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">
EXXRLvcs = 80 Ry
BndsRnXp = 25 bands
NGsBlkXp =  3  Ry
GbndRnge = 25 bands
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 shows 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 into 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 <code>ya</code> 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 band gap as expected.

[[File:GW-LDA-BN-bands.png|LDA (dashed lines) and GW (solid lines) band structures|500px]]

=== 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 = 3 Ry
PPA
BndsRnXp = 24 bands
NGsBlkXp = 3 Ry
RA
BndsRnXd = 24 bands
NGsBlkXd = 3 Ry</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.

[[File:GW-cohsex-ppa-ra.png|GW bands calculated with the COHSEX, PPA and RA|500px]]

=== 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 any case it is worthy to test during the convergence procedure.

[[File:GW-newton-secant.png|GW bands calculated with the methods Newton and Secant|500px]]