by Dmitrii Nabok, Keith Gilmore & Sven Lubeck for exciting neon
(Jupyter notebook by Megha Arya and Mara Voiculescu)
Purpose: In this tutorial you will learn how to perform a basic G0W0 calculation. As an example, the electronic band structure of bulk Si is calculated. Notice that self-consistent GW is not yet implemented in exciting.
Read the following paragraphs before starting with the rest of this tutorial!
Before running any Jupyter tutorials, please refer to the 00_before_starting.md
document on how to correctly set up the environment. This only needs to be done once. After which, the venv can be (re)activated from exciting
's root directory:
source tools/excitingjupyter/venv/excitingvenv/bin/activate
Important note: All input parameters are in atomic units!
The quasiparticle energies in the G0W0 approximation are given as a solution of the linearized quasi-particle (QP) equation Hybertsen1986 as
\begin{equation} \tag{1} \epsilon_{n\mathbf k}^{QP} = \epsilon_{n\mathbf k} + Z_{n\mathbf k} \Big[ Re \sum_{n\mathbf k} (\epsilon_{n\mathbf k})- V_{n\mathbf k}^{xc}\Big], \end{equation}where $\epsilon_{n\mathbf k}$ are the Kohn-Sham (KS) eigenvalues, $\sum_{n\mathbf k}$ and $V_{n\mathbf k}^{xc}$ are, respectively, the diagonal matrix elements of the self-energy and the exchange-correlation (xc) potential that is employed in the single-particle KS Hamiltonian. The QP renormalization factor $Z_{n\mathbf k}$ accounts for the energy-dependence of the self-energy. The self-energy $\sum (\mathbf r , \mathbf r^{\prime};\omega)$ is given by Hedin1965
\begin{equation} \tag{2} \sum (\mathbf r , \mathbf r^{\prime};\omega) = \frac{i}{2\pi} \int G_{0}(\mathbf r , \mathbf r^{\prime};\omega + \omega^{\prime}) \; W_{0}(\mathbf r , \mathbf r^{\prime}; \omega^{\prime}) \; e^{i \omega^{\prime}\eta} \; d \omega^{\prime} , \end{equation}where $G_{0}$ is the noninteracting single-particle Green function obtained from the Kohn-Sham states and $W_{0}$ is the dynamically screened Coulomb potential
\begin{equation} \tag{3} W_{0}(\mathbf r , \mathbf r^{\prime}; \omega) = \int \epsilon^{-1} (\mathbf r ,\mathbf r^{\prime\prime}; \omega) \; \nu_{C}(\mathbf r^{\prime\prime}, \mathbf r^{\prime}) \; d\mathbf r^{\prime\prime} . \end{equation}Here, $\nu_{C} ( \mathbf r, \mathbf r^{\prime}) = 1/|\mathbf r - \mathbf r^{\prime} | $ is the Coulomb potential and $\epsilon (\mathbf r ,\mathbf r^{\prime}; \omega) $ is the dielectric function calculated in the random-phase approximation (RPA).
As a first step, create a new directory named run_Si_GW and move inside it.
%%bash
mkdir -p run_Si_GW
Below is the input file (input.xml) for bulk Si, which we consider in this tutorial.
<input>
<title>Silicon: Ground-State Calculation</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal>
<basevect>5.13 5.13 0.00</basevect>
<basevect>5.13 0.00 5.13</basevect>
<basevect>0.00 5.13 5.13</basevect>
</crystal>
<species speciesfile="Si.xml" rmt="2.1">
<atom coord="0.00 0.00 0.00"></atom>
<atom coord="0.25 0.25 0.25"></atom>
</species>
</structure>
<groundstate
do="fromscratch"
rgkmax="7.0"
ngridk="3 3 3"
xctype="LDA_PW"
>
</groundstate>
</input>
The next step is writing the complete input as a string and saving it in your working directory as input.xml.
N.B.: Do not forget to replace in the input.xml the string "$EXCITINGROOT" by the actual value of the environment variable $EXCITINGROOT using the command
%%bash
cd run_Si_GW
python3 -m excitingscripts.setup.excitingroot
cd ..
In order to run exciting
, simply execute the exciting_smp binary in the running directory.
%%bash
cd run_Si_GW
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
After the run, one can inspect the INFO.OUT file for checking the convergence behaviour of the SCF cycle and for other information related to the ground-state calculation. Please note that the quality of the groundstate calculations (the starting point) plays an important role when calculating G0W0 QP corrections. Therefore one should always check the convergence of results with respect to the groundstate
parameters (such as rgkmax
, ngridk
, etc.) in production calculations.
Using the data from the previous run as an input, one performs a G0W0 calculation to obtain the QP energies.
Disable the execution of the ground-state run by changing the attribute do="fromscratch"
to do="skip"
in the groundstate
element of the input.xml file.
To initialize G0W0 calculations, a new block, corresponding to the element gw
, must be added to the input file inside the input
element. Now the file input.xml will look like the following.
<input>
<title>Silicon: GW Single-Shot Run</title>
<structure speciespath="$EXCITINGROOT/species">
<crystal>
<basevect>5.13 5.13 0.00</basevect>
<basevect>5.13 0.00 5.13</basevect>
<basevect>0.00 5.13 5.13</basevect>
</crystal>
<species speciesfile="Si.xml" rmt="2.1">
<atom coord="0.00 0.00 0.00"></atom>
<atom coord="0.25 0.25 0.25"></atom>
</species>
</structure>
<groundstate
do="skip"
rgkmax="7.0"
ngridk="3 3 3"
xctype="LDA_PW"
>
</groundstate>
<gw
taskname="g0w0"
ngridq="3 3 3"
nempty="22"
ibgw="1"
nbgw="10"
coreflag="xal"
>
<mixbasis
lmaxmb="3"
epsmb="1.d-4"
gmb="1.0"
></mixbasis>
<freqgrid
nomeg="12"
></freqgrid>
</gw>
</input>
Let us look closer at the gw
element. Attributes and elements used in this example are :
Parameter | Description |
---|---|
gw /taskname |
Specifies the task name. |
gw /ngridq |
Specifies the number of both the k and the q-points used specifically in GW calculations. Notice that these values are independent from the one specified by groundstate /ngridk . |
gw /nempty |
The number of empty states used as an input in GW calculations. Note that this value is independent from the one specified inside the groundstate element. |
gw /ibgw , gw /nbgw |
Give the range of resulting QP states to be calculated. ibgw (nbgw ) is the lowest (highest) band index for which the GW correction is calculated. Notice that the numbering of the band indexes is such that the index 1 is the one for the lowest occupied valence band. Furthermore, nbgw must be chosen to be larger than the number of occupied bands. |
gw /coreflag |
Specifies the treatment of core electrons. |
gw /mixbasis |
Specifies the parameters for the mixed product basis. |
gw /freqgrid |
Specifies the frequency grid for dynamical quantities. |
For all available parameters and options and their description, one should check the GW Section of Input Reference. It is interesting to notice that the computational setup given for the mixed-product basis and the frequency grid is rather universal and could be applied for most materials without modification.
%%bash
cd run_Si_GW
python3 -m excitingscripts.setup.excitingroot
cd ..
We are now ready to run the calculation with the usual command.
%%bash
cd run_Si_GW
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
Once the calculation has finished, the main output of the G0W0 run will be stored in GW_INFO.OUT and EVALQP.DAT. The former file contains information on the GW program work-flow. In particular, energy band-gap values for KS and QP spectra can be found at the end of this file.
The text file EVALQP.DAT contains all important quantities that are calculated during the G0W0 run. Below, we describe the structure of this file. The line "k-point # 1:" starts the block of data and contains the lattice coordinates of the irreducible k-point and the corresponding Brillouin-zone-integration weight. Then, for each electronic state (within the range determined by the ibgw
and nbgw
input attributes) one can find the following quantities:
It is important to note that the zero of the E_HF and E_GW quasiparticle energies listed in EVALQP.DAT is the KS Fermi energy (that can be found, e.g., in GW_INFO.OUT). Therefore, in order to obtain absolute QP eigenenergies one has to add to the data listed in EVALQP.DAT the corresponding KS Fermi energy.
i) Band Structure Plot
Now, we will visualize the results to compare the Kohn-Sham (KS) and quasiparticle (QP, G0W0) band structures. To calculate them, one should introduce the following changes into the input.xml file.
groundstate
block, add the new element block properties
as indicated below....
<properties>
<bandstructure>
<plot1d>
<path steps="100">
<point coord="1.0 0.0 0.0" label="Gamma"/>
<point coord="0.625 0.375 0.0" label="K"/>
<point coord="0.5 0.5 0.0" label="X"/>
<point coord="0.0 0.0 0.0" label="Gamma"/>
<point coord="0.5 0.0 0.0" label="L"/>
</path>
</plot1d>
</bandstructure>
<dos
inttype="trilin"
winddos="-0.5 0.5"
nwdos="200"
></dos>
</properties>
...
gw
element replace taskname="g0w0"
by taskname="band"
to calculate the QP band structure on the basis of results obtained from the previous G0W0 run. In this case, Fourier interpolation is used to evaluate the band energies for the k-point path specified in the properties
element....
<gw
taskname="band"
... >
</gw>
...
Finally, run exciting
one more time.
%%bash
cd run_Si_GW
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
The output files BAND.OUT and BAND-QP.OUT can now be used to plot the electronic band structure. This can be easily done using the script excitingscripts.plot.band_structure
, which is described in details in The python script "plot.band_structure".
%%bash
cd run_Si_GW
python3 -m excitingscripts.plot.band_structure -a KS GW -e -15 10
cd ..
The script produces a file PLOT.png which you can visualize on the screen with standard tools. The result should look like the following.
An interesting result is obtained if the KS and G0W0 electronic band structures are compared by setting the energy zero at the valence band maximum (VBM):
%%bash
cd run_Si_GW
python3 -m excitingscripts.plot.band_structure -a KS GW -e -15 10 -z vbM
cd ..
The result should look like the following.
This plot shows that the main difference between the two band structures is due to the conduction bands which seem to be uniformly shifted upwards by a constant value. This is the theoretical justification of the so-called scissor approximation.
Now, we will compute and compare the Kohn-Sham (KS) and quasiparticle (QP, G0W0) density of states. To calculate them, one should introduce the following changes into the input.xml file.
gw
element replace taskname="band"
by taskname="dos"
to calculate the QP density of states on the basis of results obtained from the previous G0W0 run....
<gw
taskname="dos"
... >
</gw>
...
%%bash
cd run_Si_GW
time $EXCITINGROOT/bin/exciting_smp input.xml
cd ..
The results for the KS and G0W0 density of states are stored in TDOS.OUT and TDOS-QP.OUT correspondingly. To visualize them one can simply execute the script excitingscripts.plot.dos
(discussed in The python script "plot.dos").
%%bash
cd run_Si_GW
python3 -m excitingscripts.plot.dos -a KS GW
cd ..
The script produces a file PLOT.png which you can visualize on the screen with standard tools. The result should look like the following.
The quasiparticle energies obtained from G0W0 depend on many parameters (some of them in a crucial way). In order to keep the runtimes of the above calculations short, certain parameters were reduced from their converged values. The most important parameters are
Parameter | Description | Info |
---|---|---|
gw /ngridq |
k-mesh | very crucial and should always be converged |
gw /nempty |
number of unoccupied states | very crucial and sometimes difficult to converge |
In the general case, additional attention should be payed to the proper choice of the system-dependent parameters used in the mixbasis
, freqgrid
, etc., subelements of gw
to obtain physically relevant results (see Input Reference).
Once you have successfully run the previous GW calculation, you can play around with computational parameters in order to achieve the desired convergence for the value of the electronic band gap.
groundstate
and the gw
elements.Repeat the complete GW calculation for a different system, e.g., diamond or boron nitride.