{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Back to the main [Index](../index.ipynb) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Bethe-Salpeter calculations with AbiPy and Abinit"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This lesson discusses how to calculate the macroscopic dielectric function, $\\epsilon_\\infty(\\omega)$, including excitonic effects within the Bethe-Salpeter equation (BSE). \n",
"Crystalline silicon is used as test case. \n",
"\n",
"For a more detailed description of the Abinit implementation, see the official Abinit \n",
"[BSE tutorial](https://docs.abinit.org/tutorial/bse/).\n",
"A brief description of the formalism can be found in the [BSE_notes](https://docs.abinit.org/theory/bse/).\n",
"\n",
"## Table of Contents\n",
"[[back to top](#top)]\n",
"\n",
"- [Typical BSE flowchart](#Typical-BSE-flowchart)\n",
"- [AbiPy flow for BSE with the model dielectric function](#AbiPy-flow-for-BSE-with-the-model-dielectric-function)\n",
"- [Analyzing the results](#Analyzing-the-results)\n",
"- [Converge study with respect to the $k$ point sampling](#Converge-study-with-respect-to-the-$k$-point-sampling)\n",
"- [Exercises](#Exercises)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Typical BSE flowchart \n",
"[[back to top](#top)]\n",
"\n",
"The flowchart of a typical Bethe-Salpeter run is schematically depicted in the diagram below: "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The WFK file (KSS file in old versions of Abinit) contains the Kohn-Sham (KS) wavefunctions and energies and is represented with an ellipsis. \n",
"The path on the left indicated with blue arrows represents the RPA calculation (`optdriver=3`) that produces the SCR file (see also the first lesson of the $GW$ tutorial). \n",
"Once the WFK (KSS) and the SCR file are available, we can finally contruct the BSE Hamiltonian \n",
"and solve the Bethe-Salpeter problem (the green rectangle at the bottom of the flowchart).\n",
"The construction of the Bethe-Salpeter Hamiltonian represents a significant portion of the overall CPU time due to the large number of transitions (bands and in particular $k$-points) needed for an accurate description of the frequency-dependence of the polarizability.\n",
"\n",
"For BSE computations, it is common practice to simulate the self-energy corrections by employing the scissors operator whose value can be obtained either from experiments or from *ab-initio* calculations. The scissors operator allows one to avoid a costly $GW$ calculation that should performed for all the $k$-points and bands included in the transition space (the optional path on the right indicated with yellow arrows that corresponds to `optdriver=4`).\n",
"\n",
"For this reason, in this lesson, we will employ two commonly used approximations that will\n",
"reduce considerably the computational cost of the BSE flowchart while giving reasonably accurate results:\n",
"\n",
" * The *ab-initio* $W$ is replaced by a model dielectric function \n",
" that is constructed from the GS density $n(r)$ and the additional variable `mdf_epsinf`\n",
" that gives the value of the static limit $\\epsilon_\\infty(\\omega=0)$.\n",
" This approximation allows us to bypass the blue boxes in the diagram above (`optdriver=3`)\n",
" \n",
" * The modifications introduced by the $GW$ self-energy on the initial KS band structure\n",
" are approximated with a scissor operator (`mbpt_sciss`).\n",
" This approximation allows us to bypass the yellow boxes in the diagram above (`optdriver=4`).\n",
" \n",
"Under these assumptions, the BSE flowchart reduces to a simple GS-SCF run to get $n(r)$ plus\n",
"a NSCF calculation of the band structure on a dense $k$-mesh and, finally, the solution\n",
"of the BSE problem (the green box)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## AbiPy flow for BSE with the model dielectric function\n",
"[[back to top](#top)]\n",
"\n",
"After the standard imports:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import warnings\n",
"warnings.filterwarnings(\"ignore\") # to get rid of deprecation warnings\n",
"\n",
"import abipy.abilab as abilab\n",
"abilab.enable_notebook() # This line tells AbiPy we are running inside a notebook\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": [
"We import from `lesson_bse` the function that builds the 3 input objects we are going to use to build the `Flow`:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
def make_scf_nscf_bse_inputs(ngkpt=(6, 6, 6), ecut=6, ecuteps=3,\n",
" mdf_epsinf=12.0, mbpt_sciss="0.8 eV"):\n",
" """\n",
" Build and returns three `AbinitInput` objects to perform a\n",
" GS-SCF + GS-NSCF + BSE calculation with model dielectric function.\n",
"\n",
" Args:\n",
" ngkpt: Three integers giving the number of divisions for the k-mesh.\n",
" ecut: Cutoff energy for the wavefunctions.\n",
" ecuteps: Cutoff energy for the screened interation W_{GG'}.\n",
" mdf_epsinf: Static limit of the macroscopic dielectric functions.\n",
" Used to build the model dielectric function.\n",
" mbpt_sciss: Scissors operator energy (used to open the initial KS gap).\n",
" """\n",
" multi = abilab.MultiDataset(structure=abidata.structure_from_ucell("Si"),\n",
" pseudos=abidata.pseudos("14si.pspnc"), ndtset=3)\n",
" multi.set_mnemonics(True)\n",
"\n",
" # Variables common to the three datasets.\n",
" multi.set_vars(\n",
" ecut=ecut,\n",
" nband=8,\n",
" istwfk="*1",\n",
" diemac=12.0,\n",
" #iomode=3,\n",
" )\n",
"\n",
" # SCF run to get the density.\n",
" multi[0].set_vars(tolvrs=1e-8)\n",
" multi[0].set_kmesh(ngkpt=ngkpt, shiftk=(0, 0, 0))\n",
"\n",
" # NSCF run on a randomly shifted k-mesh (improve the convergence of optical properties)\n",
" multi[1].set_vars(\n",
" iscf=-2,\n",
" nband=15,\n",
" tolwfr=1e-8,\n",
" chksymbreak=0, # Skip the check on the k-mesh.\n",
" )\n",
"\n",
" # This shift breaks the symmetry of the k-mesh.\n",
" multi[1].set_kmesh(ngkpt=ngkpt, shiftk=(0.11, 0.21, 0.31))\n",
"\n",
" # BSE run with Haydock iterative method (only resonant + W + v)\n",
" multi[2].set_vars(\n",
" optdriver=99, # BS calculation\n",
" chksymbreak=0, # To skip the check on the k-mesh.\n",
" bs_calctype=1, # L0 is constructed with KS orbitals and energies.\n",
" mbpt_sciss=mbpt_sciss, # Scissors operator used to correct the KS band structure.\n",
" bs_exchange_term=1, # Exchange term included.\n",
" bs_coulomb_term=21, # Coulomb term with model dielectric function.\n",
" mdf_epsinf=mdf_epsinf, # Parameter for the model dielectric function.\n",
" bs_coupling=0, # Tamm-Dancoff approximation.\n",
" bs_loband=2, # Lowest band included in the calculation\n",
" nband=6, # Highest band included in the calculation\n",
" bs_freq_mesh="0 6 0.02 eV", # Frequency mesh for the dielectric function\n",
" bs_algorithm=2, # Use Haydock method.\n",
" zcut="0.15 eV", # Complex shift to avoid divergences in the continued fraction.\n",
" ecutwfn=ecut, # Cutoff for the wavefunction.\n",
" ecuteps=ecuteps, # Cutoff for W and /bare v.\n",
" inclvkb=2, # The commutator for the optical limit is correctly evaluated.\n",
" )\n",
"\n",
" # Same shift as the one used in the previous dataset.\n",
" multi[2].set_kmesh(ngkpt=ngkpt, shiftk=(0.11, 0.21, 0.31))\n",
"\n",
" scf_input, nscf_input, bse_input = multi.split_datasets()\n",
"\n",
" return scf_input, nscf_input, bse_input\n",
"def build_bse_flow(options):\n",
" """\n",
" Build a flow to solve the BSE with default parameters.\n",
"\n",
" Args:\n",
" options: Command line options.\n",
"\n",
" Return:\n",
" Flow object.\n",
" """\n",
" workdir = options.workdir if (options and options.workdir) else "flow_bse"\n",
" flow = flowtk.Flow(workdir=workdir)\n",
"\n",
" # Build a Work for BSE calculation with the model dielectric function ...\n",
" scf_inp, nscf_inp, bse_inp = make_scf_nscf_bse_inputs()\n",
"\n",
" work = flowtk.BseMdfWork(scf_inp, nscf_inp, bse_inp)\n",
"\n",
" # and add it to the flow\n",
" flow.register_work(work)\n",
"\n",
" return flow\n",
"def build_bse_metallicW_flow(options):\n",
" """\n",
" Build a flow to solve the BSE with metallic screening.\n",
" Note the value of `mdf_epsinf`.\n",
"\n",
" Args:\n",
" options: Command line options.\n",
"\n",
" Return:\n",
" Flow object.\n",
" """\n",
" workdir = options.workdir if (options and options.workdir) else "flow_bse_metallicW"\n",
" flow = flowtk.Flow(workdir=workdir)\n",
"\n",
" # Model dielectric function with metallic screening\n",
" scf_inp, nscf_inp, bse_inp = make_scf_nscf_bse_inputs(ngkpt=(4, 4, 4), ecut=6, ecuteps=3,\n",
" mdf_epsinf=1.0e+12)\n",
"\n",
" flow.register_work(flowtk.BseMdfWork(scf_inp, nscf_inp, bse_inp))\n",
"\n",
" return flow\n",
"def build_bse_kconv_flow(options):\n",
" """\n",
" Build a flow to analyze the convergence of the BSE spectrum wrt k-point sampling.\n",
"\n",
" Args:\n",
" options: Command line options.\n",
"\n",
" Return:\n",
" Flow object.\n",
" """\n",
" workdir = options.workdir if (options and options.workdir) else "flow_bse_kconv"\n",
" flow = flowtk.Flow(workdir=workdir)\n",
"\n",
" # 3 works with differet ngkpt.\n",
" for nk in [4, 6, 8]:\n",
" scf_inp, nscf_inp, bse_inp = make_scf_nscf_bse_inputs(ngkpt=3 * [nk], ecut=6, ecuteps=3)\n",
" work = flowtk.BseMdfWork(scf_inp, nscf_inp, bse_inp)\n",
" flow.register_work(work)\n",
"\n",
" return flow\n",
"