The Yambo Project - User contributions [en]

Yambopy tutorial: band structures

2022-04-04T16:04:05Z

Alejandro: /* Tutorial 1. BN (semiconductor). Band structure */

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful not only to execute Quantum Espresso but also to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The class '''PwXML''' reads the data file generated by Quantum Espresso and post-process the data. The class is initiated doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. Previously to plot the bands, we define the k-points path using the function Path:

npoints = 50
path_kpoints = 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)])

It is worth to note that the path should coincide with the selected path for the QE calculations.

In order to plot we do:

xml.plot_eigen(path_kpoints)

This function will automaticall plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px | center |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively we can use the function plot_eigen_ax. This functions requires as input the figure object and with given axes, as in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic bands can be displayed using the class ProjwfcXML. Before we need to use the QE function projwfc.x to create the file '''atomic_proj.xml'''. The function projects wavefunctions onto orthogonalized atomic wavefunctions, among others functionalities. The orbital order is explained in the input documentation of the projwfc.x function (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). To run projwf.x we can use the class ProjwfcIn:

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class PwXML, we define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the boron and nitrogen orbitals

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class ProjwfcXML runs with this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | center | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). We can also plot the bands varying the size. In this case we pass only the variable selected_orbitals (see next tutorial). The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial the example of using the Classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py to check an example. Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. We have to run the file:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. In the case of representing the size as a function of the weight of the orbitals, we have represented exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script).

==Tutorial 3==

Yambopy is used mainly to run Yambo and QE calculations. Nevertheless, there are classes to analyse the results of Yambo and deals with the database generated by Yambo. In the case of GW quasi-particle calculations we can use the class YamboQPDB to find the scissor operator, plot the GW bands and to interpolated the GW bands in a smoother k-path. The example runs by doing:

python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
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)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, just the interpolation will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle resutls ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')

The first option is to plot and calculate the scissor operator. We need to select the index of the top valence band:

n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not linear would be very rare and you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

valence bands:
slope: 1.0515886598785766
conduction bands:
slope: 1.026524081134514
scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure, in the path. The first choice would be to plot exact GW calculations, without interpolation, to check that the results are meaningfull (or not). This plot is independent of the number of k-points ('''npoints''') that we want to plot in the interpolation. The class YamboQPDB finds the calculated points that belong to the path and plot them. Be aware that if we use coarse grids the class would not find any point and the function will not work.

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')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interplation of the DFT and GW band structures has a similar aspect:

ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
ks_bs.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

The '''lpratio''' can be increase if the interpolation does not work. The interpolation is the one implemented by abipy.

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally, we can compare the calculated GW eigenvalues with the interpolation.

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-04T16:03:42Z

Alejandro: /* Tutorial 3 */

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful not only to execute Quantum Espresso but also to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The class '''PwXML''' reads the data file generated by Quantum Espresso and post-process the data. The class is initiated doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. Previously to plot the bands, we define the k-points path using the function Path:

npoints = 50
path_kpoints = 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)])

It is worth to note that the path should coincide with the selected path for the QE calculations.

In order to plot we do:

xml.plot_eigen(path_kpoints)

This function will automaticall plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively we can use the function plot_eigen_ax. This functions requires as input the figure object and with given axes, as in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic bands can be displayed using the class ProjwfcXML. Before we need to use the QE function projwfc.x to create the file '''atomic_proj.xml'''. The function projects wavefunctions onto orthogonalized atomic wavefunctions, among others functionalities. The orbital order is explained in the input documentation of the projwfc.x function (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). To run projwf.x we can use the class ProjwfcIn:

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class PwXML, we define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the boron and nitrogen orbitals

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class ProjwfcXML runs with this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). We can also plot the bands varying the size. In this case we pass only the variable selected_orbitals (see next tutorial). The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial the example of using the Classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py to check an example. Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. We have to run the file:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. In the case of representing the size as a function of the weight of the orbitals, we have represented exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script).

==Tutorial 3==

Yambopy is used mainly to run Yambo and QE calculations. Nevertheless, there are classes to analyse the results of Yambo and deals with the database generated by Yambo. In the case of GW quasi-particle calculations we can use the class YamboQPDB to find the scissor operator, plot the GW bands and to interpolated the GW bands in a smoother k-path. The example runs by doing:

python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
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)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, just the interpolation will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle resutls ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')

The first option is to plot and calculate the scissor operator. We need to select the index of the top valence band:

n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not linear would be very rare and you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

valence bands:
slope: 1.0515886598785766
conduction bands:
slope: 1.026524081134514
scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure, in the path. The first choice would be to plot exact GW calculations, without interpolation, to check that the results are meaningfull (or not). This plot is independent of the number of k-points ('''npoints''') that we want to plot in the interpolation. The class YamboQPDB finds the calculated points that belong to the path and plot them. Be aware that if we use coarse grids the class would not find any point and the function will not work.

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')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interplation of the DFT and GW band structures has a similar aspect:

ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
ks_bs.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

The '''lpratio''' can be increase if the interpolation does not work. The interpolation is the one implemented by abipy.

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally, we can compare the calculated GW eigenvalues with the interpolation.

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-04T16:01:04Z

Alejandro: /* Tutorial 3 */

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful not only to execute Quantum Espresso but also to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The class '''PwXML''' reads the data file generated by Quantum Espresso and post-process the data. The class is initiated doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. Previously to plot the bands, we define the k-points path using the function Path:

npoints = 50
path_kpoints = 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)])

It is worth to note that the path should coincide with the selected path for the QE calculations.

In order to plot we do:

xml.plot_eigen(path_kpoints)

This function will automaticall plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively we can use the function plot_eigen_ax. This functions requires as input the figure object and with given axes, as in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic bands can be displayed using the class ProjwfcXML. Before we need to use the QE function projwfc.x to create the file '''atomic_proj.xml'''. The function projects wavefunctions onto orthogonalized atomic wavefunctions, among others functionalities. The orbital order is explained in the input documentation of the projwfc.x function (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). To run projwf.x we can use the class ProjwfcIn:

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class PwXML, we define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the boron and nitrogen orbitals

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class ProjwfcXML runs with this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). We can also plot the bands varying the size. In this case we pass only the variable selected_orbitals (see next tutorial). The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial the example of using the Classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py to check an example. Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. We have to run the file:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. In the case of representing the size as a function of the weight of the orbitals, we have represented exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script).

==Tutorial 3==

Yambopy is used mainly to run Yambo and QE calculations. Nevertheless, there are classes to analyse the results of Yambo and deals with the database generated by Yambo. In the case of GW quasi-particle calculations we can use the class YamboQPDB to find the scissor operator, plot the GW bands and to interpolated the GW bands in a smoother k-path. The example runs by doing:

python plot-qp.py

See below for an explanation of the tutorial. As usual, we can import the qepy and yambopy libraries:

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
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)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, just the interpolation will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle resutls ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')

The first option is to plot and calculate the scissor operator. We need to select the index of the top valence band:

n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)

Yambopy displays the fitting and also the data of the slope of each fitting. Notice that this is also a test if the GW calculations are running well. '''If the dependence is not linear would be very rare and you should double-check your results!'''

[[File:Figure 6-slope-scissor.png|400px|center]]

In this case the slope is:

valence bands:
slope: 1.0515886598785766
conduction bands:
slope: 1.026524081134514
scissor list (shift,c,v) [eV,adim,adim]: [1.8985204833551723, 1.026524081134514, 1.0515886598785766]

In addition to the scissor operator, we can plot the GW (and DFT) band structure, in the path. The first choice would be to plot exact GW calculations, without interpolation, to check that the results are meaningfull (or not).

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')

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

The interplation of the DFT and GW band structures has a similar aspect:

ks_bs, qp_bs = ydb.interpolate(lat,path,what='QP+KS',lpratio=20)
ks_bs.plot_ax(ax,legend=True,color_bands='r',label='KS')
qp_bs.plot_ax(ax,legend=True,color_bands='b',label='QP-GW')

The '''lpratio''' can be increase if the interpolation does not work. The interpolation is the one implemented by abipy.

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

Finally

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-04T15:53:30Z

Alejandro: /* Tutorial 3 */

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful not only to execute Quantum Espresso but also to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The class '''PwXML''' reads the data file generated by Quantum Espresso and post-process the data. The class is initiated doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. Previously to plot the bands, we define the k-points path using the function Path:

npoints = 50
path_kpoints = 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)])

It is worth to note that the path should coincide with the selected path for the QE calculations.

In order to plot we do:

xml.plot_eigen(path_kpoints)

This function will automaticall plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively we can use the function plot_eigen_ax. This functions requires as input the figure object and with given axes, as in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic bands can be displayed using the class ProjwfcXML. Before we need to use the QE function projwfc.x to create the file '''atomic_proj.xml'''. The function projects wavefunctions onto orthogonalized atomic wavefunctions, among others functionalities. The orbital order is explained in the input documentation of the projwfc.x function (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). To run projwf.x we can use the class ProjwfcIn:

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class PwXML, we define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the boron and nitrogen orbitals

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class ProjwfcXML runs with this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). We can also plot the bands varying the size. In this case we pass only the variable selected_orbitals (see next tutorial). The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial the example of using the Classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py to check an example. Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. We have to run the file:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. In the case of representing the size as a function of the weight of the orbitals, we have represented exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script).

==Tutorial 3==

Yambopy is used mainly to run Yambo and QE calculations. Nevertheless, there are classes to analyse the results of Yambo and deals with the database generated by Yambo. In the case of GW quasi-particle calculations we can use the class YamboQPDB to find the scissor operator, plot the GW bands and to interpolated the GW bands in a smoother k-path. As usual, we can import the qepy and yambopy libraries.

from qepy import *
from yambopy import *
import matplotlib.pyplot as plt

We define the k-points path.
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)] )

Importantly, the number of points is a free choice. We can increase the variable npoints as we wish, just the interpolation will take more time. In order to analyse GW results we need to have the file related to the basic data of our Yambo calculation ('''SAVE/ns.db1''') and the netcdf file with the quasi-particle resutls ('''ndb.QP'''). We load the data calling the respective classes:

# Read Lattice information from SAVE
lat = YamboSaveDB.from_db_file(folder='SAVE',filename='ns.db1')
# Read QP database
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='qp-gw')

[[File:Figure 6-slope-scissor.png|400px|center]]

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-04T15:45:03Z

Alejandro: /* Tutorial 3 */

The full tutorial, including the Quantum espresso and Yambo databases that we will read, can be downloaded and extracted from the yambo website:
$wget www.yambo-code.org/educational/tutorials/files/databases_qepy.tar
$tar -xvf databases_qepy.tar
$cd databases_qepy

==Tutorial 1. BN (semiconductor). Band structure==

The qepy classes are useful not only to execute Quantum Espresso but also to analyze the results. The full qepy library is imported by simply doing:

from qepy import *

===Plot Band structure===

The class '''PwXML''' reads the data file generated by Quantum Espresso and post-process the data. The class is initiated doing:

xml = PwXML(prefix='bn', path='bands')

The variable prefix corresponds to the same variable of the QE input. The folder location is indicated by variable path. Previously to plot the bands, we define the k-points path using the function Path:

npoints = 50
path_kpoints = 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)])

It is worth to note that the path should coincide with the selected path for the QE calculations.

In order to plot we do:

xml.plot_eigen(path_kpoints)

This function will automaticall plot the bands as shown below running the script:

python plot-qe-bands.py

[[File:Bands BN 1.png| 400 px |Band structure of BN calculated at the level of DFT-LDA]]

Alternatively we can use the function plot_eigen_ax. This functions requires as input the figure object and with given axes, as in the next example.

===Plot atomic orbital projected Band structure===

In addition to the band structure, useful information regarding the atomic orbital nature of the electronic bands can be displayed using the class ProjwfcXML. Before we need to use the QE function projwfc.x to create the file '''atomic_proj.xml'''. The function projects wavefunctions onto orthogonalized atomic wavefunctions, among others functionalities. The orbital order is explained in the input documentation of the projwfc.x function (https://www.quantum-espresso.org/Doc/INPUT_PROJWFC.html#idm94). To run projwf.x we can use the class ProjwfcIn:

proj = ProjwfcIn(prefix='bn')
proj.run(folder='bands')

Be aware that this can take a while in a large system with many k-points. As in the class PwXML, we define the path_kpoints and also the selected atomic orbitals to project our bands. We have chosen to represent the projection weight onto the boron and nitrogen orbitals

atom_1 = list(range(16))
atom_2 = list(range(16,32))

The full list of orbitals is written in the file '''bands/prowfc.log'''. We have also defined the figure box

import matplotlib.pyplot as plt
fig = plt.figure(figsize=(5,7))
ax = fig.add_axes( [ 0.12, 0.10, 0.70, 0.80 ])

The class ProjwfcXML runs with this example:

band = ProjwfcXML(prefix='bn',path='bands',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',selected_orbitals=atom_1,selected_orbitals_2=atom_2)

We can run now the file:

python plot-qe-orbitals.py

[[File:Bands AOP BN 1.png| 400 px | Atomic orbital projected band structure of monolayer BN]]

We have chosen the colormap 'viridis' (variable cmap). We can also plot the bands varying the size. In this case we pass only the variable selected_orbitals (see next tutorial). The colormap can be represented in many ways and qepy does not include a particular function for this. An example is:

import matplotlib as mpl
cmap=plt.get_cmap('viridis')
bx = fig.add_axes( [ 0.88, 0.10, 0.05, 0.80 ])
norm = mpl.colors.Normalize(vmin=0.,vmax=1.)
cb1 = mpl.colorbar.ColorbarBase(bx, cmap=cmap, norm=norm,orientation='vertical',ticks=[0,1])
cb1.set_ticklabels(['B', 'N'])

==Tutorial 2. Iron. Ferromagnetic metalic material==

In the case of spin-polarized calculations we can plot the spin up and down band structures. We have included in the tutorial the example of using the Classes Tasks and Flows developed in yambopy. You can check the file flow-iron.py to check an example. Yambopy includes basic predefined workflows to run the common QE and Yambo calculations. In this case we are using the class '''PwBandsTasks'''.

python flow-iron.py

In order to plot the spin-polarized bands. We have to run the file:

python plot-qe-bands.py

The class PwXML automatically detects the spin polarized case (nspin=2 in QE input file). The spin up channel is displayed with red and the spin down channel with blue. In the case of iron we have selected this path:

npoints = 50
path_kpoints = Path([ [[0.0, 0.0, 0.0 ],'G'],
[[0.0, 0.0, 1.0 ],'H'],
[[1./2,0.0,1./2.],'N'],
[[0.0, 0.0, 0.0 ],'G'],
[[1./2, 1./2, 1./2 ],'P'],
[[1./2,0.0,1./2. ],'N']], [npoints,npoints,npoints,npoints,npoints])

xml = PwXML(prefix='pw',path='bands/t0')
xml.plot_eigen(path_kpoints)

[[File:Figure 3-iron-bands.png| 600 px | center |Spin polarized band structure of iron calculated by DFT]]

The analysis of the projected atomic orbitals is also implemented. In this case the results are more cumbersome because the projection is separated in spin up and down channels. In the case of representing the size as a function of the weight of the orbitals, we have represented exclusively the ''d'' orbitals:

atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,selected_orbitals=atom_d,color='red',color_2='blue')

If we run the file:

plot-qe-orbitals-size.py

[[File:Figure 4-iron-bands-size-d-orbitals.png|400px|center|Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.]]

From the plot is clear that ''d'' orbitals are mainly localized around the Fermi energy.

Another option is to plot the orbital composition as a colormap running the file:

plot-qe-orbitals-colormap.py

[[File:Figure 5-iron-bands-colormap.png|400px|center]]

Here we have added the p and d orbitals in the analysis:

atom_p = [0,1,2]
atom_d = [3,4,5,6,7]
band = ProjwfcXML(prefix='pw',path='bands/t0',qe_version='7.0')
band.plot_eigen(ax,path_kpoints=path_kpoints,cmap='viridis',cmap2='rainbow',selected_orbitals=atom_p,selected_orbitals_2=atom_d)

The colormap bar is added in the same way as in Tutorial 1 (see script).

==Tutorial 3==

Yambopy is used mainly to run Yambo and QE calculations. Nevertheless, there are classes to analyse the results of Yambo and deals with the database generated by Yambo. In the case of GW quasi-particle calculations we can use the class YamboQPDB to find the scissor operator, plot the GW bands and to interpolated the GW bands in a smoother k-path.

[[File:Figure 6-slope-scissor.png|400px|center]]

[[File:Figure 7-GW-band-structure-non-interpolated.png|400px|center]]

[[File:Figure 8-GW-band-structure-interpolated.png|400px|center]]

[[File:Figure 8-GW-band-structure-comparison.png|400px|center]]

==Links==
* Back to [[ICTP 2022#Tutorials]]

Yambopy tutorial: band structures

2022-04-04T15:42:14Z

Alejandro: /* Tutorial 2. Iron. Ferromagnetic metalic material */

Yambopy tutorial: band structures

2022-04-04T15:38:30Z

Alejandro: /* Tutorial 1. BN (semiconductor). Band structure */

Yambopy tutorial: band structures

2022-04-04T15:36:15Z

Alejandro: /* Tutorial 2. Iron. Ferromagnetic metalic material */

Yambopy tutorial: band structures

2022-04-04T15:25:13Z

Alejandro: /* Tutorial 2. Iron. Ferromagnetic metalic material */

Yambopy tutorial: band structures

2022-04-04T15:24:42Z

Alejandro: /* Tutorial 2. Iron. Ferromagnetic metalic material */

Yambopy tutorial: band structures

2022-04-04T15:22:43Z

Alejandro: /* Tutorial 2 */

Yambopy tutorial: band structures

2022-04-04T15:14:54Z

Alejandro: /* Tutorial 2 */

Yambopy tutorial: band structures

2022-04-04T15:14:21Z

Alejandro: /* Tutorial 1. BN (semiconductor). Band structure */

Yambopy tutorial: band structures

2022-04-04T15:02:05Z

Alejandro: /* Tutorial 1. BN (semiconductor). Band structure */

Yambopy tutorial: band structures

2022-04-04T14:45:31Z

Alejandro:

Yambopy tutorial: band structures

2022-04-04T14:41:28Z

Alejandro:

Yambopy tutorial: band structures

2022-04-04T14:31:18Z

Alejandro:

Yambopy tutorial: band structures

2022-04-04T14:21:24Z

Alejandro: /* Tutorial 1 */

Yambopy tutorial: band structures

2022-04-04T10:34:49Z

Alejandro:

Yambopy tutorial: band structures

2022-04-02T20:37:51Z

Alejandro: /* Tutorial 3 */

File:Figure 8-GW-band-structure-interpolated.png

2022-04-02T20:37:01Z

Alejandro:

interpolated

File:Figure 8-GW-band-structure-comparison.png

2022-04-02T20:36:02Z

Alejandro:

GW band interpolated

File:Figure 7-GW-band-structure-non-interpolated.png

2022-04-02T20:35:30Z

Alejandro:

GW band non-interpolated

File:Figure 6-slope-scissor.png

2022-04-02T20:34:39Z

Alejandro:

Slope

Yambopy tutorial: band structures

2022-04-02T20:21:52Z

Alejandro: /* Tutorial 2 */

File:Figure 5-iron-bands-colormap.png

2022-04-02T20:21:02Z

Alejandro:

Projected iron bands

Yambopy tutorial: band structures

2022-04-02T20:16:07Z

Alejandro: /* Tutorial 2 */

File:Figure 4-iron-bands-size-d-orbitals.png

2022-04-02T20:15:29Z

Alejandro:

Iron band structure. Size is proportional to the weights of the projection on atomic d-orbitals. Red (blue) is up (down) spin polarization.

Yambopy tutorial: band structures

2022-04-02T19:59:51Z

Alejandro:

File:Figure 3-iron-bands.png

2022-04-02T19:58:18Z

Alejandro:

Spin polarized band structure of Iron

Yambopy tutorial: band structures

2022-04-01T19:37:19Z

Alejandro:

Yambopy tutorial: band structures

2022-03-31T20:38:18Z

Alejandro:

'''Tutorial 1. BN (semiconductor). Band structure'''

'''Plot Band structure'''

[[File:Bands BN 1.png|Band structure of BN calculated at the level of DFT-LDA]]

'''Plot atomic orbital projected Band structure'''

[[File:Bands AOP BN 1.png|0.5px|frame|center]]

File:Bands AOP BN 1.png

2022-03-31T20:37:59Z

Alejandro:

Atomic-projected band structure

Yambopy tutorial: band structures

2022-03-31T20:37:20Z

Alejandro:

Yambopy tutorial: band structures

2022-03-31T20:37:05Z

Alejandro: /* Tutorial 1. BN (semiconductor). Band structure */

'''Tutorial 1. BN (semiconductor). Band structure'''

'''Plot Band structure'''

<gallery>
[[File:Bands BN 1.png|Band structure of BN calculated at the level of DFT-LDA]]
</gallery>

'''Plot atomic orbital projected Band structure'''

Yambopy tutorial: band structures

2022-03-31T20:36:43Z

Alejandro:

== Tutorial 1. BN (semiconductor). Band structure ==

'''Plot Band structure'''

<gallery>
[[File:Bands BN 1.png|Band structure of BN calculated at the level of DFT-LDA]]
</gallery>

'''Plot atomic orbital projected Band structure'''

Yambopy tutorial: band structures

2022-03-31T20:36:08Z

Alejandro:

== Tutorial 1. BN (semiconductor). Band structure ==

=== '''Plot Band structure''' ===

<gallery>
[[File:Bands BN 1.png|Band structure of BN calculated at the level of DFT-LDA]]
</gallery>

=== '''Plot atomic orbital projected Band structure''' ===

Yambopy tutorial: band structures

2022-03-31T20:35:00Z

Alejandro:

== Tutorial 1. BN (semiconductor). Band structure ==

=== '''Plot Band structure''' ===

[[File:Bands BN 1.png|thumb]]

=== '''Plot atomic orbital projected Band structure''' ===

Yambopy tutorial: band structures

2022-03-31T20:34:12Z

Alejandro:

Tutorial 1. BN (semiconductor). Band structure

'''Plot Band structure'''

[[File:Bands BN 1.png|thumb]]

'''Plot atomic orbital projected Band structure'''

File:Bands BN 1.png

2022-03-31T20:33:24Z

Alejandro:

Band structure of BN

Yambopy tutorial: band structures

2022-03-31T20:31:50Z

Alejandro: Created page with "Tutorial 1. BN (semiconductor). Band structure '''Plot Band structure''' <gallery> Bands_BN_1.png|Caption1 </gallery> '''Plot atomic orbital projected Band structure'''"

Tutorial 1. BN (semiconductor). Band structure

'''Plot Band structure'''

<gallery>
Bands_BN_1.png|Caption1
</gallery>

'''Plot atomic orbital projected Band structure'''

GW tutorial. Convergence and approximations (BN)

2020-01-31T08:45:39Z

Alejandro:

We have selected 2D hexagonal boron nitride to explain the use of yambopy. Throughout this tutorial we show how to use yambopy to make efficient convergence tests, to compare different approximations and to analyse 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>

Here you will be using a predefined lattice parameter for the calculations. You can, however, run a geometry optimisation calculation with

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

which at the end will update the input files for <code>pw.x</code> with the relaxed lattice parameters and atomic coordinates.

A quick note: almost all python scripts in this tutorial can be run without any options, but they will only print out the help message. In the case of <code>gs_bn.py</code> it will print:

<source lang="bash">python gs_bn.py
usage: gs_bn.py [-h] [-r] [-s] [-n] [-n2] [-b] [-l] [-o] [-p] [-d]
[-t NTHREADS]

Test the yambopy script.

optional arguments:
-h, --help show this help message and exit
-r, --relax Structural relaxation
-s, --scf Self-consistent calculation
-n, --nscf Non-self consistent calculation
-n2, --nscf_double Non-self consistent calculation for the double grid
-b, --bands Calculate band-structure
-l, --plot Plot band-structure
-o, --orbitals Plot atomic orbital projected band-structure
-p, --phonon Phonon calculation
-d, --dispersion Phonon dispersion
-t NTHREADS, --nthreads NTHREADS
Number of threads
</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">python gw_conv_gw.py -c </source>

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

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.

<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,folder='gw_conv',run=run,ref_run=False)
</source>

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:

<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 have defined interactive run, in the folder <code>gw_conv</code>. 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, <code>yambopy</code> 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 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
</source>

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!

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

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

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

[[File:GbndRnge2.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 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:

<source lang="bash">
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,19,2,6]
</source>

We can change the <code>gs_bn.py</code> scripts to calculate a non self-consistent run for a larger '''k-grid''' (12x12x1 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 <code>gw_bn.py</code>. If we enter in the script we can check that we define a scheduler (the default is bash):

<source lang="bash">
scheduler = Scheduler.factory
</source>

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

<source lang="bash">
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]
</source>

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

<source lang="bash">
python gw_bn.py
</source>

If everything has worked fine now we will have inside the folder <code>gw/yambo</code> the netCDF file <code>ndb.QP</code> 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 <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.

<source lang="bash">
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)] )
</source>

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

<source lang="bash">
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
</source>

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

<source lang="bash">
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
</source>

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

'''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-LDA-GW.png|GW and LDA band structures |300px]]

'''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:Bands-LDA-GW-interpolated.png| Interpolated GW and LDA band structures |300px]]

'''F. Comparison non-interpolated and interpolated results.''' It is convenient to check how good is the interpolation of GW results.

[[File:Lda-gw-comparison.png| Exact and interpolated GW and LDA band structures |300px]]

=== GW convergence for all k-points ===

If we want to perform GW convergence in all the band states, we can check the
script <code>plot-gw-conv.py</code>. Here we plot the band structure for all the k-points. For this, we pack the results in json files and then use the analyzer:

<source lang="bash">
pack_files_in_folder('gw_conv')
ya = YamboAnalyser('gw_conv')
</source>

For the case of the variable NGsBlkXs, we have set:

<source lang="bash">
ya.plot_gw_all_kpoints_convergence(tag='NGsBlkXs')
</source>

[[File:Chi-cutoff-kpoints.png||500px]]

=== Approximations of the dielectric function (COHSEX, PPA) ===

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, (ii) plasmon-pole approximations (PPA), or (iii) real axis integration. We have set the same parameters for the first two options, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">
COHSEX
EXXRLvcs = 60 Ry
BndsRnXs = 40 bands
NGsBlkXs = 3 Ry

PPA
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 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>

We observe a difference between COHSEX and PPA.

[[File:Gw-chi-s.png|GW bands calculated with the COHSEX and PPA|300px]]

You are welcome to refine the calculations using better parameters and check how this difference changes. In the next tutorial, [[Bethe-Salpeter equation tutorial. Optical absorption (BN)]], you can use the scissor operator computed in section [[#GW_calculation_in_regular_k-grid|GW calculation in regular k-grid]].

GW tutorial. Convergence and approximations (BN)

2020-01-23T16:46:47Z

Alejandro: /* Plotting GW calculations. Scissor operator and GW band structure */

We have selected 2D hexagonal boron nitride to explain the use of yambopy. Throughout this tutorial we show how to use yambopy to make efficient convergence tests, to compare different approximations and to analyse 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>

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.

<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,folder='gw_conv',run=run,ref_run=False)
</source>

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:

<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 have defined interactive run, in the folder <code>gw_conv</code>. 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, <code>yambopy</code> 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 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
</source>

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!

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

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

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

[[File:GbndRnge2.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 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:

<source lang="bash">
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,19,2,6]
</source>

We can change the <code>gs_bn.py</code> scripts to calculate a non self-consistent run for a larger '''k-grid''' (12x12x1 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 <code>gw_bn.py</code>. If we enter in the script we can check that we define a scheduler (the default is bash):

<source lang="bash">
scheduler = Scheduler.factory
</source>

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

<source lang="bash">
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]
</source>

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

<source lang="bash">
python gw_bn.py
</source>

If everything has worked fine now we will have inside the folder <code>gw/yambo</code> the netCDF file <code>ndb.QP</code> 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 <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.

<source lang="bash">
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)] )
</source>

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

<source lang="bash">
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
</source>

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

<source lang="bash">
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
</source>

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

'''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-LDA-GW.png|GW and LDA band structures |300px]]

'''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:Bands-LDA-GW-interpolated.png| Interpolated GW and LDA band structures |300px]]

'''F. Comparison non-interpolated and interpolated results.''' It is convenient to check how good is the interpolation of GW results.

[[File:Lda-gw-comparison.png| Exact and interpolated GW and LDA band structures |300px]]

=== GW convergence for all k-points ===

If we want to perform GW convergence in all the band states, we can check the
script <code>plot-gw-conv.py</code>. Here we plot the band structure for all the k-points. For this, we pack the results in json files and then use the analyzer:

<source lang="bash">
pack_files_in_folder('gw_conv')
ya = YamboAnalyser('gw_conv')
</source>

For the case of the variable NGsBlkXs, we have set:

<source lang="bash">
ya.plot_gw_all_kpoints_convergence(tag='NGsBlkXs')
</source>

[[File:Chi-cutoff-kpoints.png||500px]]

=== Approximations of the dielectric function (COHSEX, PPA) ===

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, (ii) plasmon-pole approximations (PPA), or (iii) real axis integration. We have set the same parameters for the first two options, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">
COHSEX
EXXRLvcs = 60 Ry
BndsRnXs = 40 bands
NGsBlkXs = 3 Ry

PPA
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 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>

We observe a large difference between COHSEX and PPA.

[[File:Gw-chi-s.png|GW bands calculated with the COHSEX and PPA|300px]]

File:Bands-LDA-GW.png

2020-01-23T16:43:56Z

Alejandro: User created page with UploadWizard

=={{int:filedesc}}==
{{Information
|description={{en|1=Yambo tutorial image}}
|date=2020-01-23
|source={{own}}
|author=[[User:Alejandro|Alejandro]]
|permission=
|other versions=
}}

=={{int:license-header}}==
{{self|cc-by-sa-4.0}}

This file was uploaded with the UploadWizard extension.

[[Category:Uploaded with UploadWizard]]

File:Bands-LDA-GW-interpolated.png

2020-01-23T16:43:56Z

Alejandro: User created page with UploadWizard

File:Lda-gw-comparison.png

2020-01-23T16:43:56Z

Alejandro: User created page with UploadWizard

GW tutorial. Convergence and approximations (BN)

2020-01-23T15:42:05Z

Alejandro: /* Approximations of the dielectric function (COHSEX, PPA) */

We have selected 2D hexagonal boron nitride to explain the use of yambopy. Throughout this tutorial we show how to use yambopy to make efficient convergence tests, to compare different approximations and to analyse 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>

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.

<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,folder='gw_conv',run=run,ref_run=False)
</source>

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:

<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 have defined interactive run, in the folder <code>gw_conv</code>. 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, <code>yambopy</code> 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 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
</source>

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!

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

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

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

[[File:GbndRnge2.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 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:

<source lang="bash">
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,19,2,6]
</source>

We can change the <code>gs_bn.py</code> scripts to calculate a non self-consistent run for a larger '''k-grid''' (12x12x1 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 <code>gw_bn.py</code>. If we enter in the script we can check that we define a scheduler (the default is bash):

<source lang="bash">
scheduler = Scheduler.factory
</source>

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

<source lang="bash">
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]
</source>

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

<source lang="bash">
python gw_bn.py
</source>

If everything has worked fine now we will have inside the folder <code>gw/yambo</code> the netCDF file <code>ndb.QP</code> 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 <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.

<source lang="bash">
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)] )
</source>

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

<source lang="bash">
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
</source>

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

<source lang="bash">
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
</source>

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

'''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 |300px]]

'''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:bands-interpolated.png| Interpolated GW and LDA band structures |300px]]

'''F. Comparison non-interpolated and interpolated results.''' It is convenient to check how good is the interpolation of GW results.

[[File:check-interpolation.png| Exact and interpolated GW and LDA band structures |300px]]

=== GW convergence for all k-points ===

If we want to perform GW convergence in all the band states, we can check the
script <code>plot-gw-conv.py</code>. Here we plot the band structure for all the k-points. For this, we pack the results in json files and then use the analyzer:

<source lang="bash">
pack_files_in_folder('gw_conv')
ya = YamboAnalyser('gw_conv')
</source>

For the case of the variable NGsBlkXs, we have set:

<source lang="bash">
ya.plot_gw_all_kpoints_convergence(tag='NGsBlkXs')
</source>

[[File:Chi-cutoff-kpoints.png||500px]]

=== Approximations of the dielectric function (COHSEX, PPA) ===

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, (ii) plasmon-pole approximations (PPA), or (iii) real axis integration. We have set the same parameters for the first two options, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">
COHSEX
EXXRLvcs = 60 Ry
BndsRnXs = 40 bands
NGsBlkXs = 3 Ry

PPA
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 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>

We observe a large difference between COHSEX and PPA.

[[File:Gw-chi-s.png|GW bands calculated with the COHSEX and PPA|300px]]

File:Gw-chi-s.png

2020-01-23T15:41:21Z

Alejandro: User created page with UploadWizard

GW tutorial. Convergence and approximations (BN)

2020-01-16T16:35:32Z

Alejandro: /* Approximations of the dielectric function (COHSEX, PPA) */

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>

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.

<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,folder='gw_conv',run=run,ref_run=False)
</source>

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:

<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 have defined interactive run, in the folder <code>gw_conv</code>. 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, <code>yambopy</code> 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 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
</source>

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!

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

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

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

[[File:GbndRnge2.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 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:

<source lang="bash">
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,19,2,6]
</source>

We can change the <code>gs_bn.py</code> scripts to calculate a non self-consistent run for a larger '''k-grid''' (12x12x1 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 <code>gw_bn.py</code>. If we enter in the script we can check that we define a scheduler (the default is bash):

<source lang="bash">
scheduler = Scheduler.factory
</source>

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

<source lang="bash">
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]
</source>

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

<source lang="bash">
python gw_bn.py
</source>

If everything has worked fine now we will have inside the folder <code>gw/yambo</code> the netCDF file <code>ndb.QP</code> 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 <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.

<source lang="bash">
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)] )
</source>

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

<source lang="bash">
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
</source>

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

<source lang="bash">
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
</source>

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

'''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 |300px]]

'''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:bands-interpolated.png| Interpolated GW and LDA band structures |300px]]

'''F. Comparison non-interpolated and interpolated results.''' It is convenient to check how good is the interpolation of GW results.

[[File:check-interpolation.png| Exact and interpolated GW and LDA band structures |300px]]

=== GW convergence for all k-points ===

If we want to perform GW convergence in all the band states, we can check the
script <code>plot-gw-conv.py</code>. Here we plot the band structure for all the k-points. For this, we pack the results in json files and then use the analyzer:

<source lang="bash">
pack_files_in_folder('gw_conv')
ya = YamboAnalyser('gw_conv')
</source>

For the case of the variable NGsBlkXs, we have set:

<source lang="bash">
ya.plot_gw_all_kpoints_convergence(tag='NGsBlkXs')
</source>

[[File:Chi-cutoff-kpoints.png||500px]]

=== Approximations of the dielectric function (COHSEX, PPA) ===

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, (ii) plasmon-pole approximations (PPA), or (iii) real axis integration. We have set the same parameters for the first two options, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">
COHSEX
EXXRLvcs = 60 Ry
BndsRnXs = 40 bands
NGsBlkXs = 3 Ry

PPA
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 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>

We observe a large difference between COHSEX and PPA.

[[File:GW-COHSEX-PPA.png|GW bands calculated with the COHSEX and PPA|300px]]

GW tutorial. Convergence and approximations (BN)

2020-01-16T16:35:16Z

Alejandro: /* Approximations of the dielectric function (COHSEX, PPA) */

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>

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.

<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,folder='gw_conv',run=run,ref_run=False)
</source>

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:

<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 have defined interactive run, in the folder <code>gw_conv</code>. 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, <code>yambopy</code> 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 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
</source>

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!

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

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

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

[[File:GbndRnge2.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 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:

<source lang="bash">
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 3 Ry
GbndRnge = 30 bands
QPkrange = [1,19,2,6]
</source>

We can change the <code>gs_bn.py</code> scripts to calculate a non self-consistent run for a larger '''k-grid''' (12x12x1 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 <code>gw_bn.py</code>. If we enter in the script we can check that we define a scheduler (the default is bash):

<source lang="bash">
scheduler = Scheduler.factory
</source>

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

<source lang="bash">
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]
</source>

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

<source lang="bash">
python gw_bn.py
</source>

If everything has worked fine now we will have inside the folder <code>gw/yambo</code> the netCDF file <code>ndb.QP</code> 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 <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.

<source lang="bash">
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)] )
</source>

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

<source lang="bash">
lat = YamboSaveDB.from_db_file(folder='gw/SAVE',filename='ns.db1')
ydb = YamboQPDB.from_db(filename='ndb.QP',folder='gw/yambo')
</source>

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

<source lang="bash">
n_top_vb = 4
ydb.plot_scissor_ax(ax,n_top_vb)
</source>

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

'''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 |300px]]

'''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:bands-interpolated.png| Interpolated GW and LDA band structures |300px]]

'''F. Comparison non-interpolated and interpolated results.''' It is convenient to check how good is the interpolation of GW results.

[[File:check-interpolation.png| Exact and interpolated GW and LDA band structures |300px]]

=== GW convergence for all k-points ===

If we want to perform GW convergence in all the band states, we can check the
script <code>plot-gw-conv.py</code>. Here we plot the band structure for all the k-points. For this, we pack the results in json files and then use the analyzer:

<source lang="bash">
pack_files_in_folder('gw_conv')
ya = YamboAnalyser('gw_conv')
</source>

For the case of the variable NGsBlkXs, we have set:

<source lang="bash">
ya.plot_gw_all_kpoints_convergence(tag='NGsBlkXs')
</source>

[[File:Chi-cutoff-kpoints.png||500px]]

=== Approximations of the dielectric function (COHSEX, PPA) ===

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, (ii) plasmon-pole approximations (PPA), or (iii) real axis integration. We have set the same parameters for the first two options, just changing the variable name for the number of bands and the cut-off of the screening.

<source lang="bash">
COHSEX
EXXRLvcs = 60 Ry
BndsRnXs = 40 bands
NGsBlkXs = 3 Ry

PPA
EXXRLvcs = 60 Ry
BndsRnXp = 40 bands
NGsBlkXp = 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>

We observe a large difference between COHSEX and PPA.

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

[[File:GW-COHSEX-PPA.png|GW bands calculated with the COHSEX and PPA|300px]]