# Difference between revisions of "GW tutorial. Convergence and approximations (BN)"

(→Plotting GW calculations. Scissor operator and GW band structure) |
(→Plotting GW calculations. Scissor operator and GW band structure) |
||

Line 196: | Line 196: | ||

=== Plotting GW calculations. Scissor operator and GW band structure === | === Plotting GW calculations. Scissor operator and GW band structure === | ||

− | If everything has worked fine now we can start using the yambopy analysis tools. For this purpose, we have created the script <code>plot-qp.py</code>. Using this script we will be able to read the Yambo databases and plot the results. The file is structured as follows: | + | If everything has worked fine now we can start using the yambopy analysis tools. For this purpose, we have created the script <code>plot-qp.py</code>. Using this |

+ | script we will be able to read the Yambo databases and plot the results. Notice that we use matplotlib to make the plots. For all the plots we define a figure and axis by | ||

+ | doing: | ||

+ | |||

+ | <source lang="bash"> | ||

+ | fig = plt.figure(figsize=(6,4)) | ||

+ | ax = fig.add_axes( [ 0.20, 0.20, 0.70, 0.70 ]) | ||

+ | </source> | ||

+ | |||

+ | The file is structured as follows: | ||

'''A. Define a path'''. Using qepy we can define a path to plot the band structure. | '''A. Define a path'''. Using qepy we can define a path to plot the band structure. | ||

Line 225: | Line 234: | ||

[[File:Scissor-gw.png|GW results (dots) and linear fitting (solid lines) |500px]] | [[File:Scissor-gw.png|GW results (dots) and linear fitting (solid lines) |500px]] | ||

+ | '''D. Plot exact QP-GW eigenvalues in a path.''' We can also plot the band structure of calculated points (not interpolated). Yambopy | ||

+ | will find which k-points belong to a given path. We can add the LDA results for comparison. | ||

+ | <source lang="bash"> | ||

+ | ks_bs_0, qp_bs_0 = ydb.get_bs_path(lat,path) | ||

+ | ks_bs_0.plot_ax(ax,legend=True,color_bands='r',label='KS') | ||

+ | qp_bs_0.plot_ax(ax,legend=True,color_bands='b',label='QP-GW') | ||

+ | </source> | ||

− | + | [[File:bands-not-interpolated.png|GW and LDA band structures |500px]] | |

− | |||

− | |||

− | |||

− | + | '''E. Plot interpolated QP-GW eigenvalues in a path.''' In order to obtain results ready for publication or presentation, we can interpolate the GW calculations. | |

− | |||

− | |||

− | [[File: | + | [[File:bands-interpolated.png| Interpolated GW and LDA band structures |500px]] |

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

## Revision as of 16:06, 16 January 2020

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.

## Contents

### 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 `gs_bn.py`

file:

```
python gs_bn.py
python gs_bn.py -sn
```

We have set the non-self-consistent run with a wave-function cutoff of **60 Ry, 70 bands and a k-grid of 12x12x1** 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:

```
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)])
```

We obtain the band structure by doing:

`python gs_bn.py -b`

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 `YamboIn`

with the corresponding arguments:

`y = YamboIn.from_runlevel('yambo -d -p p -g n -V all',folder='gw_conv')`

The main variables of a GW calculation are:

`EXXRLvcs`

: Exchange self-energy cutoff

`NGsBlkXp`

: Cutoff of the dielectric function.

`BndsRnXp`

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

`GbndRnge`

: Self-energy. The number of bands.

We define a dictionary with all the parameters that we want to converge. Be aware of setting the right units and format for each parameter.

```
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]],''] }
```

The class `YamboIn`

includes the function `optimize`

, which is called here:

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

The function optimize has two main arguments: the convergence dictionary, and a function **run** with the instructions to run
the calculation, defined in our case like:

```
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()
```

We have defined interactive run, in the folder `gw_conv`

. We have also defined the name for each job, associated with the variable and its value.

**(b) Analysis**

Once all the calculations are finished it's time to analyse them. At this point yambopy will facilitate the analysis. Besides
the python module, `yambopy`

can also be called in the terminal to perform some post-analysis tasks:

```
$ 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
```

Calling `yambopy analysegw`

will display the help of the function:

```
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
```

Running the function selecting the bands and k-points, together with the parameter of convergence we will obtain the convergence plot.

```
yambopy analysegw -bc 5 -kc 7 -bv 4 -kv 7 gw_conv EXXRLvcs
yambopy analysegw -bc 5 -kc 7 -bv 4 -kv 7 gw_conv NGsBlkXp
yambopy analysegw -bc 5 -kc 7 -bv 4 -kv 7 gw_conv BndsRnXp
yambopy analysegw -bc 5 -kc 7 -bv 4 -kv 7 gw_conv GbndRnge
```

From the convergence plots, we can choose the set of converged 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!

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 regular k-grid

From the bandgap convergence study made above for a **k-grid** we can decide reasonable parameters. Another option is to decide a convergence threshold to establish
the accuracy of the convergence. In this tutorial, in order to make calculations lighter, we have chosen the following parameters:

```
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,12,2,6]
```

We can change the `gs_bn.py`

scripts to calculate a non self-consistent run for a larger **k-grid** (9x9x1 will do the job). We can also change the number of bands to 40 in order to speed up a bit the QE calculation.

We can just simply run the code to calculate the GW corrections for all the points of the Brillouin zone by setting the converged parameters in the script `gw_bn.py`

. If we enter in the script we can check that we define a scheduler (the default is bash):

`scheduler = Scheduler.factory`

We just need to define the run level and the variables we are going to chang:

```
y = YamboIn.from_runlevel('%s -p p -g n -V all'%yambo,folder='gw')
y['EXXRLvcs'] = [60,'Ry'] # Self-energy. Exchange
y['BndsRnXp'] = [1,40] # Screening. Number of bands
y['NGsBlkXp'] = [3,'Ry'] # Cutoff Screening
y['GbndRnge'] = [1,30] # Self-energy. Number of bands
y['QPkrange'] = [kpoint_start,kpoint_end,2,6]
```

We can now run the script and obtain the GW in the full **k-grid**.

`python gw_bn.py`

If everything has worked fine now we will have inside the folder `gw/yambo`

the netCDF file `ndb.QP`

with the results of the GW calculation. We can analyze the results and/or use them in BSE calculations.

### Plotting GW calculations. Scissor operator and GW band structure

If everything has worked fine now we can start using the yambopy analysis tools. For this purpose, we have created the script `plot-qp.py`

. Using this
script we will be able to read the Yambo databases and plot the results. Notice that we use matplotlib to make the plots. For all the plots we define a figure and axis by
doing:

```
fig = plt.figure(figsize=(6,4))
ax = fig.add_axes( [ 0.20, 0.20, 0.70, 0.70 ])
```

The file is structured as follows:

**A. Define a path**. Using qepy we can define a path to plot the band structure.

```
npoints = 10
path = 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(npoints*2),int(npoints),int(sqrt(5)*npoints)] )
```

**B. Read Yambo lattice and QP databases.**

```
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
```

**C. Plot all QP eigenvalues.** A very typical plot when analyzing GW calculations is the plot of the GW eigenvalues versus LDA eigenvalues. You just need to indicate which band index corresponds to the
top of the valence band.

```
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
```

**D. Plot exact QP-GW eigenvalues in a path.** We can also plot the band structure of calculated points (not interpolated). Yambopy
will find which k-points belong to a given path. We can add the LDA results for comparison.

```
ks_bs_0, qp_bs_0 = ydb.get_bs_path(lat,path)
ks_bs_0.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs_0.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')
```

**E. Plot interpolated QP-GW eigenvalues in a path.** In order to obtain results ready for publication or presentation, we can interpolate the GW calculations.

### 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.

```
COHSEX
BndsRnXs = 24 bands
NGsBlkXs = 3 Ry
PPA
BndsRnXp = 24 bands
NGsBlkXp = 3 Ry
RA
BndsRnXd = 24 bands
NGsBlkXd = 3 Ry
```

We have set the converged parameters and the function works by running:

`python gw_conv_bn.py -x`

We plot the band structure using the analyzer explained above.

`python gw_conv_bn.py -xp`

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