{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "# Introduction to band gaps and band structures\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "In this exercise we study some of the key properties of materials for photovoltaic applications. In the first part of the exercise you will be given a semiconductor material and you are requested to investigate:\n",
        "\n",
        "* atomic structure\n",
        "* band gap\n",
        "* band gap position\n",
        "* band structure\n",
        "* compare how different exchange correlation functionals perform\n",
        "\n",
        "We will use ASE and GPAW packages and at the end of this notebook, you will be requested to write your own scripts and submit them to the supercomputer. You will be asked to compare your results to each other and to discuss your results with other groups studying different materials.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Atomic structure\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "As you have already learnd in the previous sesion, when investigating the electronic structure of a material, the first thing to be done is to find the atomic positions by relaxing the forces.\n",
        "\n",
        "Here is some information to help you to build the ase.Atoms object:\n",
        "* Silicon crystalizes in the diamond structure with lattice constant a=5.43 \u00c5\n",
        "* Germanium crystalizes in the diamond structure with lattice constant a=5.66 \u00c5\n",
        "* Diamond has diamond structure (!) with lattice constant a=3.56 \u00c5\n",
        "* CdTe crystalizes in the zincblende structure with lattice constant a=6.48 \u00c5\n",
        "* GaAs crystalizes in the zincblende structure with lattice constant a=5.65 \u00c5\n",
        "* Monolayer BN centered in a hexagonal unit cell with a=2.5 \u00c5 (and 7 \u00c5 of vacuum at each side to prevent it from interacting with its periodic copies) and a basis of (0,0) and (0,$a / \\sqrt{3}$)\n",
        "\n",
        "The first thing you should do is to create an ase.Atoms object. In order to do so, you might find useful to use one of the crystal structures included in ase.build.bulk (hint, if you have an element of the IV group you might be interested on this link\n",
        "https://wiki.fysik.dtu.dk/ase/ase/build/build.html#module-ase.build).\n",
        "\n",
        "or you might have to create a list/array for the atomic positions and another one for the unit cell and then create an atoms object (hint: see above). You can find an example of how to build an ase.Atoms object here:\n",
        "https://wiki.fysik.dtu.dk/ase/ase/atoms.html\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# basic imports\n",
        "import numpy as np\n",
        "from ase import Atoms\n",
        "from ase.build import bulk\n",
        "from ase.visualize import view"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "atoms = ???\n",
        "label = '???'\n",
        "\n",
        "view(atoms)\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We are now going to relax the structure. To do so, we need to add a calculator, GPAW, to get DFT energies, and forces. We are going to use PBE exchange correlation functional.\n",
        "\n",
        "Since we are going to relax the unit cell, we need to use the plane wave mode, since it is the only that includes the stress-tensor. In order to do so, remember this mode requires you to specify the plane wave cut-off\n",
        "(hint: We recommend plane wave cut-off of 600 eV and the k-point mesh size could be (6,6,6) if you want it to run reasonably fast and to get a reasonable result). We will discuss convergence further in the next section.\n",
        "\n",
        "The materials we are looking at are semiconductors. Thus, the default value for the Fermi-Dirac smearing (i.e. occupations function) is too high (it is set up to 0.1 eV to work with metals). We recommend setting it to 0.01eV\n",
        "\n",
        "These links might be helpful for you:\n",
        "* https://wiki.fysik.dtu.dk/gpaw/documentation/basic.html#manual-mode\n",
        "* https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/lattice_constants/lattice_constants.html\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from gpaw import GPAW, PW, FermiDirac\n",
        "\n",
        "# We create the calculator object\n",
        "calc = GPAW(xc = 'PBE',\n",
        "            txt = label + '_relax.txt',\n",
        "            occupations = ???\n",
        "            mode = ???\n",
        "            kpts = ???\n",
        ")\n",
        "\n",
        "atoms.calc = calc"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We are going to relax the atomic positions and the unit cell at the same time. To do so, we are going to use the UnitCellFilter (see https://wiki.fysik.dtu.dk/ase/ase/constraints.html#ase.constraints.UnitCellFilter) and the BFGS (or QuasiNewton) optimizer.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.filters import UnitCellFilter\n",
        "from ase.optimize import BFGS\n",
        "\n",
        "filt = ???\n",
        "op = ???\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Make sure that you have understand the difference of optimizing a bare atoms object and using a filter!\n",
        "**Bonus**: Would you like to visualize the trajectory using ase gui, you can attach a trajectory file now by creating a new cell. After you execute the op.run cell, you can create another cell saying\n",
        "! ase gui filename.traj\n",
        "and execute it\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Run the optimization. This will take some time, do not get nervous.\n",
        "# Only if it takes longer than 4-5 minutes or if it does not print anything\n",
        "# contact us :)\n",
        "op.run(fmax=0.05)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# save the results in a file\n",
        "calc.write(label + '_gs.gpw')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Band gap and band structure\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We are now going to illustrate how to compute the band gap and obtain the band structure for a toy example with GPAW. You will be writing and submitting scripts doing your own meaningful calculations in the next section of this exercise, so do not worry now about the parameters, we know they are not a good choice and band structures look ugly :).\n",
        "\n",
        "The starting point for this section (and you might also want to use it in your own scripts) will be the PBE relaxed structure from the previous section:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "from ase.io import read\n",
        "# Add GPAW's relevant submodules again if you have restarted\n",
        "# the kernel or if you are copy and pasting to a script\n",
        "\n",
        "# read only the structure\n",
        "atoms = read(label + '_gs.gpw')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "We are now going to restart the calculator and recompute the ground state, saving it to a new gpw file. As we are dealing with small bulk system, plane wave mode is the most appropriate here.\n",
        "It is generally a good idea to choose a finer kpoint mesh for the band structure, but we are going to make the opposite choice here.\n",
        "We are also going to use LDA, which is faster but not very good at predicting bandgaps (yes, we know, you are going to get a silly value here).\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# self consistency in LDA\n",
        "calc = GPAW(mode=PW(200),\n",
        "            xc='LDA',\n",
        "            kpts=(2, 2, 2),\n",
        "            occupations=FermiDirac(0.01))\n",
        "atoms.calc = calc"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Lets use this calculator to get the energy, the *valence band maximum*, the *conduction band minimum*, and the *band gap*, as the difference of the two of the VBM and the CBM.\n",
        "\n",
        "For the VBM and CBM, we are going to use the get_homo_lumo method of the calculator. This method returns the energy of the highest Kohn-Sham occupied orbital (called HOMO here) and the energy lowest Kohn-Sham unoccupied orbital (the LUMO). We are going to compute the band gap at this level of theory from the difference between both.\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Run this cell to see the documentation of the get_homo_lumo method\n",
        "calc.get_homo_lumo?\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Potential energy\n",
        "E = atoms.get_potential_energy()\n",
        "vbm, cbm = calc.get_homo_lumo()\n",
        "\n",
        "print('E=', E)\n",
        "print('VBM=', vbm, 'CBM=', cbm)\n",
        "print('band gap=', cbm - vbm)"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Save the ground state to file\n",
        "calc.write(label + '_gs_LDA.gpw')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "### Band structure:\n",
        "Next, we calculate eigenvalues along a high symmetry path in the Brillouin zone. You can find the definition of the high symmetry k-points for the fcc lattice here:\n",
        "\n",
        "https://wiki.fysik.dtu.dk/ase/ase/dft/kpoints.html#ase.dft.kpoints.special_points\n",
        "\n",
        "If your system is in the fcc or the diamond structures, then, your path may look something like 'GXWKL'. For BN, 'GMKG'.\n",
        "\n",
        "For the band structure calculation, the density is fixed to the previously calculated ground state density, and as we want to calculate all k-points, symmetry is not used (symmetry='off').\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Restart from ground state and fix potential:\n",
        "calc = GPAW(label + '_gs_LDA.gpw').fixed_density(\n",
        "    nbands = ?\n",
        "    symmetry='off',\n",
        "    kpts={'path': ???,  # write your path here e.g. GXWKL/GMKG\n",
        "          'npoints': 60},\n",
        "    convergence=???\n",
        "    )"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Finally, we compute the band structure using ASE's band structure method, whose documentation you can find here:\n",
        "https://wiki.fysik.dtu.dk/ase/ase/dft/kpoints.html#ase.dft.band_structure.BandStructure\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Have a look at the documentation of the band structure method\n",
        "calc.band_structure?\n"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "bs = calc.band_structure()\n",
        "bs.plot(filename=label + '_bandstructure_LDA.png', show=True) # emax=10.0"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "code",
      "metadata": {},
      "source": [
        "# Save the band structure data, to discuss it the last day.\n",
        "bs.write(label + '_bandstructure_LDA.json')"
      ],
      "outputs": [],
      "execution_count": null
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## Convergence (optional but recommended)\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "What happened to the results in the previous section? Did they look reasonable, or can you tell something went wrong?\n",
        "In this section, we study the convergence of the results with the parameters that improve the completeness of the basis.\n",
        "\n",
        "**Note**: If your are running out of time to complete the exercise (i.e., you are left 20 minutes), contact us, we will help you to jump to the next section. You might be able to come back to discuss convergence in DFT in the last day of the summer school.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "Numerical convergence of DFT calculations should always be checked to avoid obtaining spurious results that are caused by a very coarse discretization. In this tutorial you can find an example on how to find a converged lattice constant for aluminum:\n",
        "https://wiki.fysik.dtu.dk/gpaw/tutorialsexercises/structureoptimization/lattice_constants/lattice_constants.html\n",
        "\n",
        "The k-point mesh and the plane wave energy cut-off in the previous section were too low.\n",
        "\n",
        "We suggest that you play around with the number of k-points and the plane wave cutoff rerunning the previous cells. Increasing their value produces better results, but also increases the computation time.\n",
        "\n",
        "Finally, we suggest you to explore the convergence of the band gap as in the tutorial for the lattice constant. To do so, You will have to write a script. You may want to use the tutorial and the previous cells as a guide.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n",
        "## The band gap with different exchange correlation functionals\n",
        "\n",
        "You are now about to complete the last part of the exercise. Now that you know how to do ground state plane wave calculations and to find the band structure of a semiconductor, we ask you to discuss the effect of choosing a functional at a given level of theory.\n",
        "We propose you to study the results with the following functionals:\n",
        "* LDA (the one you have just used)\n",
        "* PBE\n",
        "* RPBE\n",
        "* mBEEF\n",
        "\n",
        "mBEEF is an meta GGA exchange correlation functional inspired from Bayesian statistics. An essential feature of these functionals is an ensemble of functionals around the optimum one, which allows an estimate of the computational error to be easily calculated in a non-self-consistent fashion. Further description can be found in:\n",
        "* J. J. Mortensen, K. Kaasbjerg, S. L. Frederiksen, J. K. N\u00f8rskov, J. P. Sethna, and K. W. Jacobsen (2005). Phys. Rev. Lett. 95, 216401\n",
        "* Wellendorff, J., Lundgaard, K. T., Jacobsen, K. W., & Bligaard, T. (2014). The Journal of Chemical Physics, 140(14), 144107.\n",
        "\n",
        "To complete this part of the exercise, we suggest that you write scripts and submit them (i.e. write one script for each functional) so that they can run in parallel. As a guide, you can use the LDA calculations you have already done.\n",
        "\n",
        "Remember to choose a k-point mesh and an plane energy cut-off that make sense. In case of doubt, just ask.\n",
        "\n",
        "You do not need to relax the structure again, the PBE relaxed structure you have saved to a gpw file in the beginning is good enough. Your script can just read it. :)\n",
        "\n",
        "Your script should contain a ground state calculation with the functional you are studying. We suggest you print E, VBM, CBM and the gap to a file, together with the functional, so that you can use them for the discussion the last day of the summer school. You should save the ground state to a gpw file (it is a good idea to give different names to the files of different xc functionals) and restart the calculator from that file to compute the band structure. Save the plot to a png file and the data to a .json.\n"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "CAMD2022",
      "language": "python",
      "name": "camd2022"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.6.1"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}