{ "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", "
\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",
"