solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Live demo of schrpoisson-2dmaterials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This demo showcases the code developed as part of the following research work: \n",
    "\n",
    "A. Bussy, G. Pizzi, M. Gibertini, *Strain-induced polar discontinuities in 2D materials from combined first-principles and Schrödinger-Poisson simulations*, **Phys. Rev. B 96**, 165438 (2017), [doi:10.1103/PhysRevB.96.165438](http://doi.org/10.1103/PhysRevB.96.165438).\n",
    "\n",
    "Open-access: [arXiv:1705.01313](http://arxiv.org/abs/1705.01303)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "sys.path.append(os.path.join(os.path.split(os.path.abspath(__name__))[0],'code'))\n",
    "import schrpoisson2D as sp2d\n",
    "import ipywidgets as widgets\n",
    "%matplotlib notebook\n",
    "import matplotlib\n",
    "import pylab as pl\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "display(widgets.HTML(\"Using the schroedinger-Poisson 2D solver at version: {}<br>Material: SnSe\".format(sp2d.__version__)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "style = {'description_width': 'initial'}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def on_in_KbT(event):\n",
    "    # Same conversion factor used in the code, approximate\n",
    "    smearing_label_right.value = \"(Smearing in meV, ~ {:6.3f} Ry)\".format(in_KbT.value / 13.6 / 1000.)\n",
    "\n",
    "in_KbT = widgets.FloatSlider(\n",
    "    value=100,\n",
    "    min=0,\n",
    "    max=700,\n",
    "    step=10,\n",
    "    description='KbT:',\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='5.0f',\n",
    "    style=style\n",
    ")\n",
    "smearing_label_right = widgets.Label()\n",
    "display(widgets.HBox([in_KbT, smearing_label_right]))\n",
    "\n",
    "in_KbT.observe(on_in_KbT, names='value')\n",
    "on_in_KbT(None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_central = widgets.FloatSlider(\n",
    "    value=15,\n",
    "    min=1,\n",
    "    max=40,\n",
    "    step=0.1,\n",
    "    description=r'Central region length (angstrom):',\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='5.1f',\n",
    "    style=style,\n",
    "    layout=widgets.Layout(width='50%')    \n",
    ")\n",
    "in_central_strain = widgets.Dropdown(\n",
    "    options=[('0%', 0.), ('0.5%', 0.005), ('1%', 0.01), ('10%', 0.1)],\n",
    "    value=0.,\n",
    "    description='Strain:',\n",
    "    layout=widgets.Layout(width='20%')\n",
    ")\n",
    "display(widgets.HBox([in_central, in_central_strain]))\n",
    "\n",
    "in_outer = widgets.FloatSlider(\n",
    "    value=18,\n",
    "    min=1,\n",
    "    max=40,\n",
    "    step=0.1,\n",
    "    description=r'Outer region length (angstrom):',\n",
    "    disabled=False,\n",
    "    continuous_update=False,\n",
    "    orientation='horizontal',\n",
    "    readout=True,\n",
    "    readout_format='5.1f',\n",
    "    style=style,\n",
    "    layout=widgets.Layout(width='50%')\n",
    ")\n",
    "in_outer_strain = widgets.Dropdown(\n",
    "    options=[('0%', 0.), ('0.5%', 0.005), ('1%', 0.01), ('10%', 0.1)],\n",
    "    value=0.1,\n",
    "    description='Strain:',\n",
    "    layout=widgets.Layout(width='20%')\n",
    ")\n",
    "display(widgets.HBox([in_outer, in_outer_strain]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_input(KbT, central_length,outer_length, central_strain, outer_strain):\n",
    "    return {\n",
    "        \"calculation\": \"single_point\",\n",
    "        \"smearing\" : True,\n",
    "        \"KbT\" : KbT, \n",
    "        \"delta_x\" : 0.2,\n",
    "        \"max_iteration\" : 1000,\n",
    "        \"nb_of_states_per_band\" : 2,\n",
    "        \"plot_fit\" : False,\n",
    "        \"out_dir\" : \"single_point_output\",\n",
    "        \"setup\" : { \"slab1\" : { \"strain\" : central_strain,\n",
    "                                \"width\" : outer_length/2.,\n",
    "                                \"polarization\" : \"positive\",\n",
    "                            },\n",
    "                    \"slab2\" : { \"strain\" : outer_strain,\n",
    "                                \"width\" : central_length,\n",
    "                                \"polarization\" : \"positive\",\n",
    "                            },\n",
    "                    \"slab3\" : { \"strain\" : central_strain,\n",
    "                                \"width\" : outer_length/2.,\n",
    "                                \"polarization\" : \"positive\",\n",
    "                            },\n",
    "                },\n",
    "    }\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "material_SnSne = {\n",
    "    \"0.00\": { \"alpha_xx\": 10.22,\n",
    "       \"x_lat\": 4.408,\n",
    "       \"y_lat\": 4.288,\n",
    "       \"polarization_charge\": 0.0,\n",
    "       \"vacuum_level\": 3.471695,\n",
    "       \"valence_gamma\": {\"energy\": -1.508255819,\n",
    "                         \"conf_mass\": 1.755925079,\n",
    "                         \"DOS_mass\": 2.7330924,\n",
    "                         \"degeneracy\": 1\n",
    "                         },\n",
    "       \"valence_gamma-X\": {\"energy\": -0.8832657543,\n",
    "                           \"conf_mass\": 0.125168454,\n",
    "                           \"DOS_mass\": 0.159070572,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"valence_gamma-Y\": {\"energy\": -1.064317364,\n",
    "                           \"conf_mass\": 0.109554071,\n",
    "                           \"DOS_mass\": 0.159511425,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"conduction_gamma\": {\"energy\": 0.6090962485,\n",
    "                            \"conf_mass\": 2.741100107,\n",
    "                            \"DOS_mass\": 2.994158008,\n",
    "                            \"degeneracy\": 1\n",
    "                           },\n",
    "       \"conduction_gamma-X\": {\"energy\": 0.0902128713,\n",
    "                              \"conf_mass\": 0.110807373,\n",
    "                              \"DOS_mass\": 0.190387132,\n",
    "                              \"degeneracy\": 2\n",
    "                             },\n",
    "       \"conduction_gamma-Y\": {\"energy\": 0.05720314215,\n",
    "                              \"conf_mass\": 0.131542031,\n",
    "                              \"DOS_mass\": 0.130378261,\n",
    "                              \"degeneracy\": 2\n",
    "                             }\n",
    "     },\n",
    "    \"0.005\": { \"alpha_xx\":  9.85,\n",
    "       \"polarization_charge\": 0.0296,\n",
    "       \"vacuum_level\": 3.456272,\n",
    "       \"valence_gamma\": {\"energy\": -1.53133791,\n",
    "                         \"conf_mass\": 1.663041043,\n",
    "                         \"DOS_mass\": 2.721339909,\n",
    "                         \"degeneracy\": 1\n",
    "                         },\n",
    "       \"valence_gamma-X\": {\"energy\": -0.9222633399,\n",
    "                           \"conf_mass\": 0.128689309,\n",
    "                           \"DOS_mass\": 0.173073867,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"valence_gamma-Y\": {\"energy\": -1.129934252,\n",
    "                           \"conf_mass\": 0.11396212,\n",
    "                           \"DOS_mass\": 0.169760624,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"conduction_gamma\": {\"energy\": 0.5952654761,\n",
    "                            \"conf_mass\": 2.57832201,\n",
    "                            \"DOS_mass\": 2.822577572,\n",
    "                            \"degeneracy\": 1\n",
    "                           },\n",
    "       \"conduction_gamma-X\": {\"energy\": 0.07651969791,\n",
    "                              \"conf_mass\": 0.113872233,\n",
    "                              \"DOS_mass\": 0.195572846,\n",
    "                              \"degeneracy\": 2\n",
    "                             },\n",
    "       \"conduction_gamma-Y\": {\"energy\": 0.03308521035,\n",
    "                              \"conf_mass\": 0.134383779,\n",
    "                              \"DOS_mass\": 0.136358334,\n",
    "                              \"degeneracy\": 2\n",
    "                             }\n",
    "    },\n",
    "    \"0.01\": {  \"alpha_xx\":  9.43,\n",
    "       \"polarization_charge\": 0.06357,\n",
    "       \"vacuum_level\": 3.4408595,\n",
    "       \"valence_gamma\": {\"energy\": -1.54984677,\n",
    "                         \"conf_mass\": 1.56510041,\n",
    "                         \"DOS_mass\": 2.70985957,\n",
    "                         \"degeneracy\": 1\n",
    "                         },\n",
    "       \"valence_gamma-X\": {\"energy\": -0.964184483,\n",
    "                           \"conf_mass\": 0.13291259,\n",
    "                           \"DOS_mass\": 0.191355,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"valence_gamma-Y\": {\"energy\": -1.198008621,\n",
    "                           \"conf_mass\": 0.11936355,\n",
    "                           \"DOS_mass\": 0.18232315,\n",
    "                           \"degeneracy\": 2\n",
    "                          },\n",
    "       \"conduction_gamma\": {\"energy\": 0.5788854169,\n",
    "                            \"conf_mass\": 2.38658757,\n",
    "                            \"DOS_mass\": 2.66126322,\n",
    "                           \"degeneracy\": 1\n",
    "                           },\n",
    "       \"conduction_gamma-X\": {\"energy\": 0.06282279948,\n",
    "                              \"conf_mass\": 0.117584668,\n",
    "                              \"DOS_mass\": 0.20252247,\n",
    "                              \"degeneracy\": 2\n",
    "                             },\n",
    "       \"conduction_gamma-Y\": {\"energy\": 0.01621894616,\n",
    "                              \"conf_mass\": 0.137705331,\n",
    "                              \"DOS_mass\": 0.143839918,\n",
    "                              \"degeneracy\": 2\n",
    "                             }\n",
    "    },\n",
    "     \"0.10\": {\n",
    "         \"epsilon_xx\" : 4.390845042,\n",
    "   \"epsilon_yy\": 4.609306158,\n",
    "   \"alpha_xx\":  5.40,\n",
    "   \"alpha_yy\":  5.74 ,\n",
    "   \"polarization\": -0.1774925,\n",
    "   \"polarization_mod\": 0.3733465,\n",
    "   \"polarization_units\": \"C/m^2\",\n",
    "   \"polarization_charge\": 0.38643,\n",
    "   \"vacuum_level\":  3.183275,\n",
    "   \"band_gap\": 1.425,\n",
    "   \"valence_gamma\": {\"energy\": -1.846872487,\n",
    "                     \"position\": [0.0,0.0],\n",
    "                     \"conf_mass\": 0.804425486,\n",
    "                     \"DOS_mass\": 2.747762669,\n",
    "                     \"degeneracy\": 1\n",
    "                     },\n",
    "   \"valence_gamma-X\": {\"energy\": -1.639183729,\n",
    "                       \"position\": [0.5126292283,0.0],\n",
    "                       \"conf_mass\": 0.205658144,\n",
    "                       \"DOS_mass\": 1.059071873,\n",
    "                       \"degeneracy\": 2\n",
    "                      },\n",
    "   \"valence_gamma-Y\": {\"energy\": -2.150722701,\n",
    "                       \"position\": [0.0,0.5741904823],\n",
    "                       \"conf_mass\": 0.271355674,\n",
    "                       \"DOS_mass\": 0.688994731,\n",
    "                       \"degeneracy\": 2\n",
    "                      },\n",
    "   \"conduction_gamma\": {\"energy\": 0.3081798586,\n",
    "                        \"position\": [0.0,0.0],\n",
    "                        \"conf_mass\": 0.98472771,\n",
    "                        \"DOS_mass\": 1.593144026,\n",
    "                        \"degeneracy\": 1\n",
    "                       },\n",
    "   \"conduction_gamma-X\": {\"energy\": -0.2144276807,\n",
    "                          \"position\": [0.5231571587,0.0],\n",
    "                          \"conf_mass\": 0.182798895,\n",
    "                          \"DOS_mass\": 0.213362934,\n",
    "                          \"degeneracy\": 2\n",
    "                         },\n",
    "   \"conduction_gamma-Y\": {\"energy\": -0.205542465 ,\n",
    "                          \"position\": [0.0,0.6099056797],\n",
    "                          \"conf_mass\": 0.205831494,\n",
    "                          \"DOS_mass\": 0.3196511,\n",
    "                          \"degeneracy\": 2\n",
    "                         }\n",
    "    }\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def onChange(event):\n",
    "    #pl.get_current_fig_manager().toolbar.mode == '': # Not in Zoom-Mode\n",
    "    xl = total_free_charge_plot.get_xlim()\n",
    "    yl = total_free_charge_plot.get_ylim()\n",
    "    current_aspect = (yl[1] - yl[0]) / (xl[1] - xl[0])\n",
    "    new_xlim = free_charge_plot.get_xlim()\n",
    "    total_free_charge_plot.set_xlim(new_xlim)\n",
    "    total_free_charge_plot.set_ylim([0,(new_xlim[1] - new_xlim[0]) * current_aspect])    \n",
    "    #with calc_output:\n",
    "    #    print(event, dir(event))\n",
    "    #    print(event.canvas, event.guiEvent, event.name, event.renderer)\n",
    "\n",
    "\n",
    "def show_output(retval):\n",
    "    global total_free_charge_plot, free_charge_plot\n",
    "    cmap = matplotlib.cm.RdYlBu\n",
    "    graphs_output.clear_output()\n",
    "    graphs2_output.clear_output()    \n",
    "    with graphs_output:\n",
    "        density_profile = retval['out_files']['density_profile']['data']\n",
    "        pl.figure()\n",
    "        free_charge_plot = pl.subplot(2,1,1)\n",
    "        x = density_profile[:,0]\n",
    "        pl.title(\"Free charge density profile (you can zoom and pan)\")\n",
    "        pl.ylabel(\"Charge density (cm$^{-2}$)\")\n",
    "        pl.plot(x, density_profile[:,1], label=\"Free electrons\", color=cmap(0))\n",
    "        pl.plot(x, density_profile[:,2], label=\"Free holes\", color=cmap(cmap.N))        \n",
    "        pl.plot(x, np.zeros(len(x)), color=(0.3,0.3,0.3,1.))\n",
    "        free_charge_plot.set_xlim([x[0], x[-1]])\n",
    "        pl.axvline(0., color=(0.9,0.9,0.9,1.))\n",
    "        pl.axvline(in_outer.value/2., color=(0.7,0.7,0.7,1.))\n",
    "        pl.axvline(in_outer.value/2. + in_central.value, color=(0.7,0.7,0.7,1.))\n",
    "        pl.axvline(in_outer.value + in_central.value, color=(0.9,0.9,0.9,1.))\n",
    "        pl.legend()        \n",
    "                \n",
    "        total_free_charge_plot = pl.subplot(2,1,2)\n",
    "        pl.title(\"Net free charge\")\n",
    "        net_charge = density_profile[:,2] - density_profile[:,1]\n",
    "        pl.imshow(np.vstack([net_charge] * 2), \n",
    "                  extent = (x[0], x[-1], 0, (x[-1]-x[0])/10.), \n",
    "                  #aspect=10.,\n",
    "                  interpolation='bilinear',\n",
    "                  cmap=cmap,\n",
    "                 )\n",
    "        total_free_charge_plot.yaxis.set_visible(False)\n",
    "        total_free_charge_plot.axes.can_zoom = False\n",
    "        total_free_charge_plot.axes.can_pan = False        \n",
    "        \n",
    "        pl.axvline(0., color=(0.9,0.9,0.9,1.))\n",
    "        pl.axvline(in_outer.value/2., color=(0.7,0.7,0.7,1.))\n",
    "        pl.axvline(in_outer.value/2. + in_central.value, color=(0.7,0.7,0.7,1.))\n",
    "        pl.axvline(in_outer.value + in_central.value, color=(0.9,0.9,0.9,1.))\n",
    "\n",
    "        pl.xlabel(\"Position ($\\AA$)\")\n",
    "        \n",
    "        pl.show()\n",
    "        pl.connect('draw_event', onChange)  \n",
    "\n",
    "    with graphs2_output:\n",
    "        #line_colors = {\n",
    "        #    'conduction': (252/255.,141/255.,89/255.,1.),\n",
    "        #    'valence': (145/255.,191/255.,219/255.,1.),\n",
    "        #}\n",
    "        line_colors = {\n",
    "            'conduction': ['#fef0d9','#fdd49e','#fdbb84','#fc8d59','#ef6548','#d7301f','#990000'][::-1],\n",
    "            'valence': ['#f1eef6','#d0d1e6','#a6bddb','#74a9cf','#3690c0','#0570b0','#034e7b'][::-1]\n",
    "        }\n",
    "        band_data = retval['out_files']['band_data']['metadata']\n",
    "        x = band_data['x']\n",
    "        pl.figure()\n",
    "        s = pl.subplot(1,1,1)\n",
    "        pl.plot(x, [0.] * len(x), color=(171/255.,221/255.,164/255.,1.), linestyle='-', linewidth=0.5) # Fermi energy\n",
    "\n",
    "        for band_type in ['valence', 'conduction']:\n",
    "            for i, the_band in enumerate(band_data[band_type]):\n",
    "                color = line_colors[band_type][(i*2)%len(line_colors[band_type])]\n",
    "                pl.plot(x, np.array(the_band['profile']) - band_data['e_fermi'], color=color, linewidth=2)\n",
    "                # print(band_type, len(the_band['states']))\n",
    "                for state in the_band['states']:\n",
    "                    #print(state['energy'] - band_data['e_fermi'])\n",
    "                    pl.plot(x, np.array(state['profile']) - band_data['e_fermi'], color=color, linewidth=0.5)\n",
    "                    pl.plot(x, [state['energy'] - band_data['e_fermi']]*len(x), color=color, linewidth=0.5, linestyle='--')                \n",
    "\n",
    "\n",
    "        pl.xlabel(\"Position ($\\AA$)\")\n",
    "        pl.ylabel(\"Band energy (eV)\")\n",
    "        pl.title(\"Band profiles and eigenstates\")\n",
    "        s.set_xlim([x[0], x[-1]])\n",
    "        pl.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def my_callback(**kwargs):\n",
    "    status_output.value = \"Simulation step {} completed.\".format(kwargs['step'])\n",
    "\n",
    "def on_run(event):    \n",
    "    calc_output.clear_output()\n",
    "    graphs_output.clear_output()\n",
    "    graphs2_output.clear_output()    \n",
    "    run_btn.disabled = True\n",
    "    with calc_output:\n",
    "        retval = sp2d.main_run(\n",
    "            matprop=material_SnSne, \n",
    "            input_dict=create_input(\n",
    "                KbT=in_KbT.value / 1000. / 13.6, # The code wants Ry, internally uses this conversion factor \n",
    "                central_length=in_central.value,\n",
    "                outer_length=in_outer.value,\n",
    "                central_strain=in_central_strain.value,\n",
    "                outer_strain=in_outer_strain.value,\n",
    "            ),\n",
    "            callback=my_callback)\n",
    "        run_btn.disabled = False\n",
    "        status_output.value += \"<br><strong>Simulation completed!\"\n",
    "        show_output(retval)\n",
    "        #print(retval)\n",
    "\n",
    "run_btn = widgets.Button(\n",
    "    description='Run simulation',\n",
    "    disabled=False,\n",
    "    button_style='success', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Run with the 2D Schroedinger-Poisson solver',\n",
    "    icon='check'\n",
    ")\n",
    "run_btn.on_click(on_run)\n",
    "display(run_btn)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def on_show_raw(event):\n",
    "    if len(raw_group.children) > 0:\n",
    "        raw_group.children = []\n",
    "        show_raw_btn.description = \"Show raw output\"\n",
    "        raw_group.box_style = ''\n",
    "    else:\n",
    "        raw_group.children = [calc_output]\n",
    "        show_raw_btn.description = \"Hide raw output\"\n",
    "        raw_group.box_style = 'warning'\n",
    "    \n",
    "status_output = widgets.HTML()\n",
    "display(status_output)\n",
    "\n",
    "show_raw_btn =  widgets.Button(\n",
    "    description='Show raw output',\n",
    "    disabled=False,\n",
    "    button_style='', # 'success', 'info', 'warning', 'danger' or ''\n",
    "    tooltip='Show raw output',\n",
    "    icon=''\n",
    ")\n",
    "show_raw_btn.on_click(on_show_raw)\n",
    "display(show_raw_btn)\n",
    "calc_output = widgets.Output()\n",
    "raw_group = widgets.VBox([], box_style='')\n",
    "display(raw_group)\n",
    "graphs_output = widgets.Output()\n",
    "display(graphs_output)\n",
    "graphs2_output = widgets.Output()\n",
    "display(graphs2_output)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.7.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}