plotly_bandDiagram.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:e3e16318ef005448635d5e2e1f3a481f8911e81fd14ff6d53a749c4fe26d60f1"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Band Diagram and Denisty of States of Silicon"
     ]
    },
    {
     "cell_type": "heading",
     "level": 2,
     "metadata": {},
     "source": [
      "with Plotly's Python API library and PyMatgen"
     ]
    },
    {
     "cell_type": "heading",
     "level": 6,
     "metadata": {},
     "source": [
      "Germain Salvato-Vallverdu (germain.vallverdu@univ-pau.fr)"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This notebook will go over an example for plotting the density of states and the band diagram of Silicon using python with [pymatgen](http://pymatgen.org/) and [plotly](https://plot.ly/) packages.\n",
      "\n",
      "* [pymatgen](http://pymatgen.org/) : (Python Materials Genomics) A robust, open-source Python library for materials analysis.\n",
      "* [plotly](https://plot.ly/) : A platform for publishing beautiful, interactive graphs from Python to the web."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "from pymatgen.io.vaspio.vasp_output import Vasprun    # read vasprun.xml output file of VASP\n",
      "from pymatgen.electronic_structure.core import Spin"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import plotly.plotly as pltly      # plotting functions\n",
      "import plotly.tools as tls         # plotly tools\n",
      "import plotly.graph_objs as go     # plot and configuration tools : Scatter, Line, Layout"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "1) Plot the density of states"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, read projected density of states (DOS) from a VASP calculation on ``\"./DOS/vasprun.xml\"`` file using pymatgen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": true,
     "input": [
      "dosrun = Vasprun(\"./DOS/vasprun.xml\")\n",
      "spd_dos = dosrun.complete_dos.get_spd_dos()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Set up the scatter plots for the total DOS and the contribution of 3s and 3p atomic orbital to the total DOS."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# total DOS\n",
      "trace_tdos = go.Scatter(\n",
      "    x=dosrun.tdos.densities[Spin.up],\n",
      "    y=dosrun.tdos.energies - dosrun.efermi,\n",
      "    mode=\"lines\",\n",
      "    name=\"total DOS\",\n",
      "    line=go.Line(color=\"#444444\"),\n",
      "    fill=\"tozeroy\"\n",
      ")\n",
      "# 3s contribution to the total DOS\n",
      "trace_3s = go.Scatter(\n",
      "    x=spd_dos[\"S\"].densities[Spin.up],\n",
      "    y=dosrun.tdos.energies - dosrun.efermi,\n",
      "    mode=\"lines\",\n",
      "    name=\"3s\",\n",
      "    line=go.Line(color=\"red\")\n",
      ")\n",
      "# 3p contribution to the total DOS\n",
      "trace_3p = go.Scatter(\n",
      "    x=spd_dos[\"P\"].densities[Spin.up],\n",
      "    y=dosrun.tdos.energies - dosrun.efermi,\n",
      "    mode=\"lines\",\n",
      "    name=\"3p\",\n",
      "    line=go.Line(color=\"green\")\n",
      ")\n",
      "dosdata = go.Data([trace_tdos, trace_3s, trace_3p])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Customize axes and general aspect of the plot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dosxaxis = go.XAxis(\n",
      "    title=\"Density of states\",\n",
      "    showgrid=True,\n",
      "    showline=True,\n",
      "    range=[.01, 3],\n",
      "    mirror=\"ticks\",\n",
      "    ticks=\"inside\",\n",
      "    linewidth=2,\n",
      "    tickwidth=2\n",
      ")\n",
      "dosyaxis = go.YAxis(\n",
      "    title=\"$E - E_f \\quad / \\quad \\\\text{eV}$\",\n",
      "    showgrid=True,\n",
      "    showline=True,\n",
      "    ticks=\"inside\",\n",
      "    mirror='ticks',\n",
      "    linewidth=2,\n",
      "    tickwidth=2,\n",
      "    zerolinewidth=2\n",
      ")\n",
      "doslayout = go.Layout(\n",
      "    title=\"Density of states of Silicon\",\n",
      "    xaxis=dosxaxis,\n",
      "    yaxis=dosyaxis\n",
      ")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dosfig = go.Figure(data=dosdata, layout=doslayout)\n",
      "plot_url = pltly.plot(dosfig, filename=\"DOS_Si\", auto_open=False)\n",
      "tls.embed(plot_url)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\"seamless=\"seamless\" src=\"https://plot.ly/~gvallverdu/209.embed\" height=\"525\" width=\"100%\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 6,
       "text": [
        "<plotly.tools.PlotlyDisplay at 0x7f6fc7d6ee48>"
       ]
      }
     ],
     "prompt_number": 6
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "2) Plot the band diagram"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "First, read bands from a VASP calculation on ``\"./Bandes/vasprun.xml\"`` file using pymatgen."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "run = Vasprun(\"./Bandes/vasprun.xml\", parse_projected_eigen = True)\n",
      "bands = run.get_band_structure(\"./Bandes/KPOINTS\", line_mode=True, efermi=dosrun.efermi)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 7
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Look for the boundaries of the band diagram in order to set up y axes range."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "emin = 1e100\n",
      "emax = -1e100\n",
      "for spin in bands.bands.keys():\n",
      "    for band in range(bands.nb_bands):\n",
      "        emin = min(emin, min(bands.bands[spin][band]))\n",
      "        emax = max(emax, max(bands.bands[spin][band]))\n",
      "emin = emin - bands.efermi - 1\n",
      "emax = emax - bands.efermi + 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Each band is plotted using a scatter plot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "kptslist = [k for k in range(len(bands.kpoints))]\n",
      "bandTraces = list()\n",
      "for band in range(bands.nb_bands):\n",
      "    bandTraces.append(\n",
      "        go.Scatter(\n",
      "            x=kptslist,\n",
      "            y=[e - bands.efermi for e in bands.bands[Spin.up][band]],\n",
      "            mode=\"lines\",\n",
      "            line=go.Line(color=\"#666666\"),\n",
      "            showlegend=False\n",
      "        )\n",
      "    )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Add vertical lines at each high symmetry k-points and a label at the bottom."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "labels = [r\"$L$\", r\"$\\Gamma$\", r\"$X$\", r\"$U,K$\", r\"$\\Gamma$\"]\n",
      "step = len(bands.kpoints) / (len(labels) - 1)\n",
      "# vertical lines\n",
      "vlines = list()\n",
      "for i, label in enumerate(labels):\n",
      "    vlines.append(\n",
      "        go.Scatter(\n",
      "            x=[i * step, i * step],\n",
      "            y=[emin, emax],\n",
      "            mode=\"lines\",\n",
      "            line=go.Line(color=\"#111111\", width=1),\n",
      "            showlegend=False\n",
      "        )\n",
      "    )\n",
      "# Labels of highsymetry k-points are added as Annotation object\n",
      "annotations = list()\n",
      "for i, label in enumerate(labels):\n",
      "    annotations.append(\n",
      "        go.Annotation(\n",
      "            x=i * step, y=emin,\n",
      "            xref=\"x1\", yref=\"y1\",\n",
      "            text=label,\n",
      "            xanchor=\"center\", yanchor=\"top\",\n",
      "            showarrow=False\n",
      "        )\n",
      "    )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Customize axes and general aspect of the plot."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bandxaxis = go.XAxis(\n",
      "    title=\"k-points\",\n",
      "    range=[0, len(bands.kpoints)],\n",
      "    showgrid=True,\n",
      "    showline=True,\n",
      "    ticks=\"\", \n",
      "    showticklabels=False,\n",
      "    mirror=True,\n",
      "    linewidth=2\n",
      ")\n",
      "bandyaxis = go.YAxis(\n",
      "    title=\"$E - E_f \\quad / \\quad \\\\text{eV}$\",\n",
      "    range=[emin, emax],\n",
      "    showgrid=True,\n",
      "    showline=True,\n",
      "    zeroline=True,\n",
      "    mirror=\"ticks\",\n",
      "    ticks=\"inside\",\n",
      "    linewidth=2,\n",
      "    tickwidth=2,\n",
      "    zerolinewidth=2\n",
      ")\n",
      "bandlayout = go.Layout(\n",
      "    title=\"Bands diagram of Silicon\",\n",
      "    xaxis=bandxaxis,\n",
      "    yaxis=bandyaxis,\n",
      "    annotations=go.Annotations(annotations)\n",
      ")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "bandfig = go.Figure(data=bandTraces + vlines, layout=bandlayout)\n",
      "plot_url = pltly.plot(bandfig, filename=\"Bands_Si\", auto_open=False)\n",
      "tls.embed(plot_url)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\"seamless=\"seamless\" src=\"https://plot.ly/~gvallverdu/249.embed\" height=\"525\" width=\"100%\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 24,
       "text": [
        "<plotly.tools.PlotlyDisplay at 0x7f6fc5a50f98>"
       ]
      }
     ],
     "prompt_number": 24
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "3) Use subplots to plot DOS and bands together"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We will now make a figure with both the band diagram and the density of states using the ``make_subplots`` facility. First, we set up a figure with two columns, one row. At the end, the two plots will share y axis."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dosbandfig = tls.make_subplots(rows=1, cols=2, shared_yaxes=True)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "This is the format of your plot grid:\n",
        "[ (1,1) x1,y1 ]  [ (1,2) x2,y1 ]\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 25
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "We use the previously defined traces for the band and the densities of state and add them to the figure object.\n",
      "\n",
      "* The bands are plotted on the left subplot (1, 1)\n",
      "* The densities of states are plotted on the right subplot (1, 2)"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# add the bands\n",
      "for btrace in bandTraces:\n",
      "    dosbandfig.append_trace(btrace, 1, 1)\n",
      "# add vlines for specific k-points\n",
      "for vline in vlines:\n",
      "    dosbandfig.append_trace(vline, 1, 1)\n",
      "# add the densities\n",
      "dosbandfig.append_trace(trace_tdos, 1, 2)\n",
      "dosbandfig.append_trace(trace_3s, 1, 2)\n",
      "dosbandfig.append_trace(trace_3p, 1, 2)    "
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Customize axes and general aspect of the plot using previously defined axis and layout options. The ``\"domain\"`` of each subplot is also define."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "dosbandfig[\"layout\"].update(\n",
      "    go.Layout(\n",
      "        title=\"Bands diagram and density of states of Silicon\",\n",
      "        xaxis1=bandxaxis,\n",
      "        yaxis1=bandyaxis,\n",
      "        xaxis2=dosxaxis,\n",
      "        annotations=go.Annotations(annotations)\n",
      "    )\n",
      ")\n",
      "# adjust size of subplots\n",
      "dosbandfig[\"layout\"][\"xaxis1\"][\"domain\"] = [0., 0.7]\n",
      "dosbandfig[\"layout\"][\"xaxis2\"][\"domain\"] = [0.702, 1.]\n",
      "# add some specific options\n",
      "dosbandfig[\"layout\"][\"yaxis1\"][\"mirror\"] = \"allticks\"\n",
      "dosbandfig[\"layout\"][\"xaxis2\"][\"mirror\"] = \"allticks\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 27
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "plot_url = pltly.plot(dosbandfig, filename=\"DOS_bands_Si\", auto_open=False)\n",
      "tls.embed(plot_url)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "html": [
        "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\"seamless=\"seamless\" src=\"https://plot.ly/~gvallverdu/269.embed\" height=\"525\" width=\"100%\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 28,
       "text": [
        "<plotly.tools.PlotlyDisplay at 0x7f6fc5319668>"
       ]
      }
     ],
     "prompt_number": 28
    },
    {
     "cell_type": "heading",
     "level": 3,
     "metadata": {},
     "source": [
      "4) Atomic orbital contributions in bands using a color scale"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As previously done for the density of states, the contribution of 3s and 3p atomic orbital may be highlighted using a color scale. The main idea, from [this example](https://github.com/vossjo/ase-espresso/wiki/Band-structure-calculation-example), is to normalize atomic orbital contributions and build the RGB code of the color from these contributions.\n",
      "\n",
      "Thus, we first compute atomic orbital normalized contributions from projected bands :"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# extract projected bands\n",
      "name = \"Si\"\n",
      "pbands = bands.get_projections_on_elts_and_orbitals({name: [\"s\", \"p\", \"d\"]})\n",
      "\n",
      "# compute contributions\n",
      "contrib = np.zeros((bands.nb_bands, len(bands.kpoints), 3))                        \n",
      "for band in range(bands.nb_bands):\n",
      "    for k in range(len(bands.kpoints)):\n",
      "        sc = pbands[Spin.up][band][k][name][\"s\"]**2\n",
      "        pc = pbands[Spin.up][band][k][name][\"p\"]**2\n",
      "        dc = pbands[Spin.up][band][k][name][\"d\"]**2\n",
      "        tot = sc + pc + dc\n",
      "        if tot != 0.0:\n",
      "            contrib[band, k, 0] = sc / tot\n",
      "            contrib[band, k, 1] = pc / tot\n",
      "            contrib[band, k, 2] = dc / tot"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now for each band, for each couple of consecutive points, we plot a line plot with a specific color, defined from the atomic orbital contributions to the current band."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "colorBands = list() # will contain the list of all lines\n",
      "nkpts = len(bands.kpoints)\n",
      "for band in range(bands.nb_bands):\n",
      "    eband = [e - bands.efermi for e in bands.bands[Spin.up][band]]\n",
      "    for k in range(nkpts - 1):\n",
      "        red, green, blue = [int(255 * (contrib[band, k, i] + contrib[band, k+1, i])/2) for i in range(3)]\n",
      "        colorBands.append(\n",
      "            go.Scatter(\n",
      "                x=[k, k+1],\n",
      "                y=[eband[k], eband[k+1]],\n",
      "                mode=\"lines\",\n",
      "                line=go.Line(color=\"rgb({}, {}, {})\".format(red, green, blue)),\n",
      "                showlegend=False\n",
      "            )\n",
      "        )"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "As previously, two subplots are used to plot the band diagram on the left and the density of states on the right."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# set up a new figure with two subplots\n",
      "colorbandfig = tls.make_subplots(rows=1, cols=2, shared_yaxes=True)\n",
      "# add the bands in the first subplot\n",
      "for btrace in colorBands:\n",
      "    colorbandfig.append_trace(btrace, 1, 1)\n",
      "# add vlines for specific k-points in the first subplot\n",
      "for vline in vlines:\n",
      "    colorbandfig.append_trace(vline, 1, 1)\n",
      "# add the densities in the second subplot\n",
      "colorbandfig.append_trace(trace_tdos, 1, 2)\n",
      "colorbandfig.append_trace(trace_3s, 1, 2)\n",
      "colorbandfig.append_trace(trace_3p, 1, 2)\n",
      "# Layout configuration\n",
      "colorbandfig[\"layout\"].update(\n",
      "    go.Layout(\n",
      "        title=\"Bands diagram and density of states of Silicon\",\n",
      "        xaxis1=bandxaxis,\n",
      "        yaxis1=bandyaxis,\n",
      "        xaxis2=dosxaxis,\n",
      "        annotations=go.Annotations(annotations)\n",
      "    )\n",
      ")\n",
      "# adjust size of subplots\n",
      "colorbandfig[\"layout\"][\"xaxis1\"][\"domain\"] = [0., 0.7]\n",
      "colorbandfig[\"layout\"][\"xaxis2\"][\"domain\"] = [0.702, 1.]\n",
      "# add some specific options\n",
      "colorbandfig[\"layout\"][\"yaxis1\"][\"mirror\"] = \"allticks\"\n",
      "colorbandfig[\"layout\"][\"xaxis2\"][\"mirror\"] = \"allticks\"\n",
      "# add a custom legend\n",
      "legend = go.Legend(\n",
      "    x=.98, y=.98,\n",
      "    xanchor=\"right\", yanchor=\"top\",\n",
      "    bordercolor='#333', borderwidth=1\n",
      ")\n",
      "colorbandfig[\"layout\"][\"legend\"] = legend\n",
      "# do the plot\n",
      "plot_url = pltly.plot(colorbandfig, filename=\"DOS_bands_Si_color\", auto_open=False)\n",
      "tls.embed(plot_url)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "This is the format of your plot grid:\n",
        "[ (1,1) x1,y1 ]  [ (1,2) x2,y1 ]\n",
        "\n"
       ]
      },
      {
       "html": [
        "<iframe id=\"igraph\" scrolling=\"no\" style=\"border:none;\"seamless=\"seamless\" src=\"https://plot.ly/~gvallverdu/306.embed\" height=\"525\" width=\"100%\"></iframe>"
       ],
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 31,
       "text": [
        "<plotly.tools.PlotlyDisplay at 0x7f6fc591a0b8>"
       ]
      }
     ],
     "prompt_number": 31
    }
   ],
   "metadata": {}
  }
 ]
}