Small changes to DOS notebook.

osscar-org · Apr 27, 2023 · 19013f3 · 19013f3
1 parent db269f6
commit 19013f3
Showing 1 changed file with 38 additions and 12 deletions.
diff --git a/notebook/band-theory/density_of_states.ipynb b/notebook/band-theory/density_of_states.ipynb
@@ -113,7 +113,9 @@
     "import matplotlib.pyplot as plt\n",
     "from ipywidgets import Button, RadioButtons, Layout, IntSlider, HBox, Text, VBox, Checkbox, Label, FloatSlider, Output, HTML\n",
     "from datetime import datetime\n",
-    "\n",
+    "import logging\n",
+    "logger = logging.getLogger()\n",
+    "logger.setLevel(logging.DEBUG)\n",
     "%matplotlib widget"
    ]
   },
@@ -254,7 +256,7 @@
     "\n",
     "# nkpt = IntSlider(value=4, min=4, max=15,description=\"Number of kpoints per dimension\\n, $N_k$\" , display='flex',\n",
     "#     flex_flow='row', continuous_update=False)\n",
-    "grange = IntSlider(value=0, min=0, max=3, description=\"Gvector range:\", style=style)\n",
+    "grange = IntSlider(value=0, min=0, max=3, description=\"G-vector range:\", style=style)\n",
     "gcov = FloatSlider(value=0.5, min=0.1, max=1.0, description=\"Guassian covariance:\", style=style)\n",
     "\n",
     "def prettify(label):\n",
@@ -463,7 +465,7 @@
     ")\n",
     "nkpt = IntSlider(value=4, min=4, max=15, description=\"\", style={'description_width': 'initial'}, continuous_update=False)\n",
     "nkpt_box = HBox([VBox([line1, line2]), nkpt])\n",
-    "nbin = IntSlider(value=30, min=30, max=500, description=\"Number of bins:\", layout=Layout(width=\"300px\"), style={'description_width': 'initial'})\n",
+    "nbin = IntSlider(value=30, min=30, max=500, description=\"Number of bins/eV (rounded):\", layout=Layout(width=\"300px\"), style={'description_width': 'initial'})\n",
     "gstd = FloatSlider(value=0.01, min=0.01, max=0.1, step=0.01, description=\"Gaussian $\\sigma$ (eV):\", layout=Layout(width=\"300px\"), style={'description_width': 'initial'})\n",
     "\n",
     "#All buttons\n",
@@ -498,7 +500,7 @@
     "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
     "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
     "\n",
-    "    eigs = _compute_dos(kpts, G, 1)\n",
+    "    eigs = _compute_dos(kpts, G, grange.value)\n",
     "\n",
     "    dosx = np.linspace(0, 10, 500)\n",
     "    dosy = lti(real_lattice_bohr, eigs, dosx)\n",
@@ -536,9 +538,13 @@
     "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
     "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
     "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
-    "    eigs = _compute_dos(kpts, G, 1)\n",
-    "    \n",
-    "    hy, hx = np.histogram(eigs.ravel(), bins=nbin.value)\n",
+    "    eigs = _compute_dos(kpts, G, grange.value)\n",
+    "    logging.debug(\"eigs max = {}\".format(eigs.max()))\n",
+    "    num_bins=int((eigs.max()-eigs.min())*(nbin.value))\n",
+    "    logging.debug(\"num bins = {}\".format(num_bins))\n",
+    "\n",
+    "    hy, hx = np.histogram(eigs.ravel(),bins=num_bins )\n",
+    "\n",
     "    hy = hy/np.sum(hy*np.diff(hx))*np.shape(eigs)[-1]\n",
     "    \n",
     "    lhist = ax_dos.barh(hx[:-1]+np.diff(hx)[0], hy, color='yellow', edgecolor='black', \n",
@@ -575,13 +581,33 @@
     "    shape = (nkpt.value, nkpt.value, nkpt.value)\n",
     "    G = Cell(real_lattice_bohr).reciprocal()*2*np.pi\n",
     "    kpts = np.dot(monkhorst_pack(shape), G).reshape(shape + (3,))\n",
-    "    eigs = _compute_dos(kpts, G, 1)\n",
+    "    eigs = _compute_dos(kpts, G, grange.value)\n",
     "    \n",
-    "    gx = np.linspace(-0.03, 1.5, 70)\n",
+    "    # gx = np.linspace(-0.03, 1.5, 70)\n",
+    "    gx = np.linspace(-0.03, 0.5, 70)\n",
+    "\n",
     "    gy = 0.0*gx\n",
     "    \n",
+    "    knum=nkpt.value\n",
+    "    Gnum=grange.value\n",
+    "    \n",
+    "    start = time.time()\n",
+    "\n",
+    "    # pdf_vals = norm(eigs.ravel().reshape(-1,1), scale=gstd).pdf(gx)\n",
+    "    # gy = np.sum(pdf_vals, axis=0)\n",
+    "    # pdf_vals = norm(eigs.ravel(), scale=gstd).pdf(gx)\n",
+    "    # gy = np.einsum('ijkg,i->j', pdf_vals.reshape(nkpt.value, nkpt.value, nkpt.value, grange.value, -1), np.ones(len(gx)))\n",
+    "\n",
     "    for eig in eigs.ravel():\n",
     "        gy += norm(eig, gstd.value).pdf(gx)\n",
+    "    # pdf_vals = np.array([norm(eigs.ravel(), scale=gstd).pdf(x) for x in gx])\n",
+    "    # gy = np.sum(pdf_vals, axis=-1)\n",
+    "    end = time.time()\n",
+    "    with output:\n",
+    "        print(\"(eigs.shape)={}. gy.shape={}\".format(eigs.shape,gy.shape))\n",
+    "\n",
+    "        print(\"time taken =\")\n",
+    "        print(end - start)\n",
     " \n",
     "    gy = gy/np.size(eigs)*np.shape(eigs)[-1]\n",
     "    lgas, = ax_dos.plot(gy, gx, 'k--', label=\"Gaussian smearing\")\n",
@@ -756,7 +782,7 @@
     {
      "data": {
       "application/vnd.jupyter.widget-view+json": {
-       "model_id": "744f6dfc398849a584678ca1fe3d5602",
+       "model_id": "7047ce7ca26549aa803d7999620554d6",
        "version_major": 2,
        "version_minor": 0
       },
@@ -776,7 +802,7 @@
     "method4 = HBox([btclear, Label(value=\"(Clear the calculated results)\")])\n",
     "\n",
     "label1 = HTML(value = f\"<b><font color='red'>Choose a method to calculate the DOS:</b>\")\n",
-    "display(VBox([output,VBox([ nkpt_box, label1, method1, method2, method3])]))\n",
+    "display(VBox([output,VBox([ nkpt_box,grange, label1, method1, method2, method3])]))\n",
     "# display(HBox([VBox([output_bs, nkpt_box, label1, method1, method2, method3]), VBox([output, method4])]))"
    ]
   },
@@ -835,7 +861,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.10"
+   "version": "3.10.6"
   },
   "voila": {
    "authors": "Dou Du and Giovanni Pizzi",