{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Back to the main [Index](../index.ipynb) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# $G_0W_0$ band structure with star-function interpolation\n", "[[back to top](#top)]\n", "\n", "Standard functionals (LDA and GGA), systematically underestimate band gaps, giving values\n", "that are about 30-40% smaller than experimental data.\n", "The inability of standard Kohn-Sham (KS) theory to give band gaps close to experiment is often referred to as the **band-gap problem**.\n", "\n", "From a theoretical point of view this is not surprising since KS eigenvalues are not supposed to give the correct band energies.\n", "The band structure of a crystal is rigorously defined as the energies needed to add or subtract electrons from the many-body system\n", "which, in turn, are related to the difference between total energies of **many-body states** differing by one electron.\n", "\n", "An alternative, more traditional, approach to the study of exchange-correlation effects in\n", "many-body systems is provided by Many-Body Perturbation Theory (MBPT) which defines a rigorous approach to the description of excited-state properties, based on the Green's function formalism.\n", "In this lesson, we discuss how to use the MBPT part of ABINIT to compute the band-structure of silicon\n", "within the so-called $G_0W_0$ approximation.\n", "For a very brief introduction to the many-body formalism, please consult the \n", "[MBPT_notes](https://docs.abinit.org/theory/mbt/)\n", "\n", "### Related ABINIT variables\n", "\n", " * optdriver\n", " * ecuteps\n", " * ecutsigx\n", " * nband\n", " * gwcalctyp\n", " * gw_qprange\n", " * all gw** variables\n", " \n", "## Table of Contents\n", "[[back to top](#top)]\n", "\n", "- [Description of the lesson](#Description-of-the-lesson)\n", "- [Analysis of the results](#Analysis-of-the-results)\n", "- [Interpolating the QP corrections](#Interpolating-the-QP-corrections)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Description of the lesson\n", "[[back to top](#top)]\n", "\n", "In this lesson, we construct an AbiPy flow made of two works.\n", "The first work is a standard KS band-structure calculation consisting of\n", "an initial GS calculation to get the density, followed by two NSCF calculations.\n", "The first NSCF task computes the KS eigenvalues on a high-symmetry path in the BZ,\n", "whereas the second NSCF task employs a homogeneous k-mesh so that one can compute\n", "the DOS from the KS eigenvalues.\n", "\n", "The second work represents the real GW workflow that uses the density computed in the first task of\n", "the previous work to compute the KS bands for *many empty states*.\n", "The WFK file produced in this step is then used to compute the screened interaction $W$.\n", "Finally, we perform a self-energy calculation that uses the $W$ produced\n", "in the previous step and the WFK file to compute the matrix elements of the self-energy and\n", "the $G_0W_0$ corrections for all the k-points in the IBZ and 8 bands (4 occupied + 4 empty).\n", "Once the flow is completed, we can interpolate the $G_0W_0$ corrections.\n", "\n", "Do not worry if there are steps of the entire procedure that are not clear to you.\n", "$GW$ calculations are much more complicated than standard KS band structures and\n", "the main goal of this lesson is to give you an overview of the Abipy capabilities." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's start with the usual import section:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Use this at the beginning of your script so that your code will be compatible with python3\n", "from __future__ import print_function, division, unicode_literals\n", "\n", "import numpy as np\n", "import warnings\n", "warnings.filterwarnings(\"ignore\") # to get rid of deprecation warnings\n", "\n", "from abipy import abilab\n", "abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook\n", "\n", "import abipy.flowtk as flowtk\n", "\n", "# This line configures matplotlib to show figures embedded in the notebook.\n", "# Replace `inline` with `notebook` in classic notebook\n", "%matplotlib inline \n", "\n", "# Option available in jupyterlab. See https://github.com/matplotlib/jupyter-matplotlib\n", "#%matplotlib widget " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "and now we import the `make_inputs` function that generates the inputs files required by our `Flow`.\n", "Let's have a look at the code:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "
\n", "def make_inputs(ngkpt, dos_ngkpt=(6, 6, 6), paral_kgb=0):\n",
" """\n",
" Crystalline silicon: calculation of the G0W0 band structure with the scissors operator.\n",
"\n",
" Args:\n",
" ngkpt: list of 3 integers. Abinit variable defining the k-point sampling.\n",
" dos_ngkpt: list of 3 integers. k-point sampling for DOS.\n",
" paral_kgb: Option used to select the eigensolver in the GS part.\n",
"\n",
" Return:\n",
" Six AbinitInput objects:\n",
"\n",
" [0]: Ground state run to get the density.\n",
" [1]: NSCF run to get the KS band structure on a high-symmetry k-path.\n",
" [2]: NSCF run with a homogeneous sampling of the BZ to compute the KS DOS with `dos_ngkpt`.\n",
" [3]: NSCF run with empty states to prepare the GW steps.\n",
" [4]: Calculation of the screening from the WFK file computed in dataset 4.\n",
" [5]: Use the SCR file computed at step 5 and the WFK file computed in dataset 4 to get the GW corrections.\n",
" """\n",
" multi = abilab.MultiDataset(abidata.cif_file("si.cif"),\n",
" pseudos=abidata.pseudos("14si.pspnc"), ndtset=6)\n",
"\n",
" # Add mnemonics to input file.\n",
" multi.set_mnemonics(True)\n",
"\n",
" # This grid is the most economical and we will use it for the GS.\n",
" # Note that it does not contain Gamma point so we cannot use it to\n",
" # compute the QP corrections at the Gamma point.\n",
" scf_kmesh = dict(\n",
" ngkpt=ngkpt,\n",
" shiftk=[0.5, 0.5, 0.5,\n",
" 0.5, 0.0, 0.0,\n",
" 0.0, 0.5, 0.0,\n",
" 0.0, 0.0, 0.5]\n",
" )\n",
"\n",
" # k-point sampling for DOS (gamma-centered)\n",
" dos_kmesh = dict(\n",
" ngkpt=dos_ngkpt,\n",
" shiftk=[0.0, 0.0, 0.0])\n",
"\n",
" # This grid contains the Gamma point, which is the point at which\n",
" # we will compute the (direct) band gap.\n",
" gw_kmesh = dict(\n",
" ngkpt=ngkpt,\n",
" shiftk=[0.0, 0.0, 0.0,\n",
" 0.0, 0.5, 0.5,\n",
" 0.5, 0.0, 0.5,\n",
" 0.5, 0.5, 0.0]\n",
" )\n",
"\n",
" # Global variables\n",
" ecut = 6\n",
" multi.set_vars(\n",
" ecut=ecut,\n",
" istwfk="*1",\n",
" paral_kgb=paral_kgb,\n",
" gwpara=2,\n",
" )\n",
"\n",
" # Dataset 1 (GS run to get the density)\n",
" multi[0].set_kmesh(**scf_kmesh)\n",
" multi[0].set_vars(\n",
" tolvrs=1e-6,\n",
" nband=4,\n",
" )\n",
" multi[0].set_kmesh(**scf_kmesh)\n",
"\n",
" # Dataset 2 (NSCF run with k-path)\n",
" multi[1].set_vars(iscf=-2,\n",
" tolwfr=1e-12,\n",
" nband=8,\n",
" )\n",
" multi[1].set_kpath(ndivsm=8)\n",
"\n",
" # Dataset 3 (DOS NSCF)\n",
" multi[2].set_vars(iscf=-2,\n",
" tolwfr=1e-12,\n",
" nband=8,\n",
" )\n",
" multi[2].set_kmesh(**dos_kmesh)\n",
"\n",
" # Dataset 4 (NSCF run to produce WFK for GW, note the presence of empty states)\n",
" multi[3].set_vars(iscf=-2,\n",
" tolwfr=1e-12,\n",
" nband=35,\n",
" )\n",
" multi[3].set_kmesh(**gw_kmesh)\n",
"\n",
" # Dataset3: Calculation of the screening.\n",
" multi[4].set_vars(\n",
" optdriver=3,\n",
" nband=25,\n",
" ecutwfn=ecut,\n",
" symchi=1,\n",
" inclvkb=0, # Disable [Vnl, r] contribution for q --> 0\n",
" ecuteps=4.0,\n",
" )\n",
" multi[4].set_kmesh(**gw_kmesh)\n",
"\n",
" multi[5].set_vars(\n",
" optdriver=4,\n",
" nband=10,\n",
" ecutwfn=ecut,\n",
" ecuteps=4.0,\n",
" ecutsigx=6.0,\n",
" symsigma=1,\n",
" gw_qprange=-4, # Compute GW corrections for all kpts in IBZ,\n",
" # all occupied states and 4 empty states,\n",
" )\n",
" multi[5].set_kmesh(**gw_kmesh)\n",
"\n",
" return multi.split_datasets()\n",
"
def build_g0w0_flow(options=None, ngkpt=(2, 2, 2)):\n",
" """\n",
" Build and return a flow with two works.\n",
" The first work is a standard KS band-structure calculation that consists of\n",
" an initial GS calculation to get the density followed by two NSCF calculations.\n",
"\n",
" The first NSCF task computes the KS eigenvalues on a high-symmetry path in the BZ,\n",
" whereas the second NSCF task employs a homogeneous k-mesh so that one can compute\n",
" the DOS from the KS eigenvalues.\n",
"\n",
" The second work represents the real GW workflow that uses the density computed in the first task of\n",
" the previous work to compute the KS bands for many empty states.\n",
" The WFK file produced in this step is then used to compute the screened interaction $W$.\n",
" Finally, we perform a self-energy calculation that uses the $W$ produced\n",
" in the previous step and the WFK file to compute the matrix elements of the self-energy and\n",
" the $G_0W_0$ corrections for all the k-points in the IBZ and 8 bands (4 occupied + 4 empty)\n",
" """\n",
"\n",
" # Call make_input to build our 4 input objects.\n",
" scf, bands_nscf, dos_nscf, gw_nscf, scr, sig = make_inputs(ngkpt=ngkpt)\n",
"\n",
" workdir = options.workdir if (options and options.workdir) else "flow_g0w0"\n",
" flow = flowtk.Flow(workdir=workdir)\n",
"\n",
" # Add KS band structure work (SCF-GS followed by two NSCF runs\n",
" # (the first one is done on a k-path, the second on the IBZ to compute the DOS\n",
" work0 = flowtk.BandStructureWork(scf, bands_nscf, dos_inputs=dos_nscf)\n",
" flow.register_work(work0)\n",
"\n",
" # Create new Work for GW\n",
" work1 = flowtk.Work()\n",
"\n",
" # NSCF run with empty states\n",
" gw_nscf_task = work1.register_nscf_task(gw_nscf, deps={work0[0]: "DEN"})\n",
"\n",
" # SCR run with WFK produced by previous task.\n",
" scr_task = work1.register_scr_task(scr, deps={gw_nscf_task: "WFK"})\n",
"\n",
" # SIGMA task (requires WFK with empty states and SCR file)\n",
" sigma_task = work1.register_sigma_task(sig, deps={gw_nscf_task: "WFK", scr_task: "SCR"})\n",
"\n",
" # Add GW work to flow.\n",
" flow.register_work(work1)\n",
"\n",
" return flow\n",
"
\n", " | e0 | \n", "qpe | \n", "vxcme | \n", "sigxme | \n", "sigcmee0 | \n", "ze0 | \n", "qpeme0 | \n", "
---|---|---|---|---|---|---|---|
0 | \n", "-6.232 | \n", "-5.693 | \n", "-10.362 | \n", "-17.428 | \n", "7.925 | \n", "0.629 | \n", "0.540 | \n", "
1 | \n", "5.633 | \n", "6.264 | \n", "-11.160 | \n", "-13.117 | \n", "2.755 | \n", "0.790 | \n", "0.631 | \n", "
2 | \n", "5.633 | \n", "6.264 | \n", "-11.160 | \n", "-13.117 | \n", "2.755 | \n", "0.790 | \n", "0.631 | \n", "
3 | \n", "5.633 | \n", "6.264 | \n", "-11.160 | \n", "-13.117 | \n", "2.755 | \n", "0.790 | \n", "0.631 | \n", "
4 | \n", "8.148 | \n", "9.444 | \n", "-9.967 | \n", "-5.517 | \n", "-2.817 | \n", "0.793 | \n", "1.296 | \n", "
5 | \n", "8.148 | \n", "9.444 | \n", "-9.967 | \n", "-5.517 | \n", "-2.817 | \n", "0.793 | \n", "1.296 | \n", "
6 | \n", "8.148 | \n", "9.444 | \n", "-9.967 | \n", "-5.517 | \n", "-2.817 | \n", "0.793 | \n", "1.296 | \n", "
7 | \n", "8.526 | \n", "9.980 | \n", "-10.766 | \n", "-5.897 | \n", "-3.037 | \n", "0.793 | \n", "1.454 | \n", "