fermi_integrate_np.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "https://notebook.rcc.uchicago.edu/files/acs.chemmater.9b05047/Data/bulk/dft/mu.py\n",
    "\n",
    "numpy-discussion, Oct 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.06 ms ± 64.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "def fermi(E,mu,kT):\n",
    "    return 1./(np.exp((E-mu)/kT) + 1)\n",
    "\n",
    "#data = genfromtxt(\"./dos/BiVO4.dos-202018-sm003\")\n",
    "#energy=data[:,0] #in eV\n",
    "#DOS=data[:,1]    #states/unit cell/ eV\n",
    "\n",
    "## cutoff some portion of DOS\n",
    "#index = argmax(energy > -11.0)\n",
    "#energy=energy[index:]\n",
    "#DOS=DOS[index:]\n",
    "#print(energy,DOS)\n",
    "\n",
    "\n",
    "energy = np.array([1, 2, 3, 4, 5], dtype=float)\n",
    "DOS = np.array([0.1, 0.2, 0.3, 0.4, 0.5])\n",
    "index = 0\n",
    "\n",
    "\n",
    "##k=1.3806488E-23 #J/K\n",
    "#k=8.6173324e-5 #eV/K\n",
    "#T= 300 # K\n",
    "#volume = 314.65333308 # ang^3, BiVO4 tet, pbe \n",
    "#volume_cm = volume*1e-24\n",
    "#numelect = 136 # BiVO4, tet, pbe\n",
    "#kT=k*T #eV\n",
    "\n",
    "kT = 1\n",
    "numelect = 2\n",
    "volume_cm = 1\n",
    "\n",
    "\n",
    "# range of mus\n",
    "#mu_init=10.3 #eV\n",
    "#mu_step=0.0001 #eV\n",
    "#res=np.zeros(npts)\n",
    "#mu_all=np.zeros(npts)\n",
    "\n",
    "\n",
    "npts = 17\n",
    "mu_init = 0\n",
    "mu_step = 0.01\n",
    "mu_all = np.zeros(npts)  #np.arange(0, npts)*mu_step + mu_init\n",
    "res = np.zeros(npts)\n",
    "\n",
    "for i in np.arange(0, npts):\n",
    "    mu=mu_init+i*mu_step\n",
    "    mu_all[i]=mu\n",
    "    integrand=DOS*fermi(energy,mu,kT)\n",
    "    res[i]=(np.trapz(integrand,energy)-numelect)/volume_cm #in cm-3\n",
    "\n",
    "allres_baseline = np.column_stack((mu_all,res))\n",
    "\n",
    "allres_baseline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# tweak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "97.8 µs ± 3.27 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n"
     ]
    }
   ],
   "source": [
    "%%timeit\n",
    "\n",
    "def fermi(E,mu,kT):\n",
    "    return 1./(np.exp((E-mu)/kT) + 1)\n",
    "\n",
    "#data = genfromtxt(\"./dos/BiVO4.dos-202018-sm003\")\n",
    "#energy=data[:,0] #in eV\n",
    "#DOS=data[:,1]    #states/unit cell/ eV\n",
    "\n",
    "## cutoff some portion of DOS\n",
    "#index = argmax(energy > -11.0)\n",
    "#energy=energy[index:]\n",
    "#DOS=DOS[index:]\n",
    "#print(energy,DOS)\n",
    "\n",
    "\n",
    "energy = np.array([1, 2, 3, 4, 5], dtype=float)\n",
    "DOS = np.array([0.1, 0.2, 0.3, 0.4, 0.5])\n",
    "index = 0\n",
    "\n",
    "\n",
    "##k=1.3806488E-23 #J/K\n",
    "#k=8.6173324e-5 #eV/K\n",
    "#T= 300 # K\n",
    "#volume = 314.65333308 # ang^3, BiVO4 tet, pbe \n",
    "#volume_cm = volume*1e-24\n",
    "#numelect = 136 # BiVO4, tet, pbe\n",
    "#kT=k*T #eV\n",
    "\n",
    "kT = 1\n",
    "numelect = 2\n",
    "volume_cm = 1\n",
    "\n",
    "\n",
    "# range of mus\n",
    "#mu_init=10.3 #eV\n",
    "#mu_step=0.0001 #eV\n",
    "#res=np.zeros(npts)\n",
    "#mu_all=np.zeros(npts)\n",
    "\n",
    "\n",
    "npts = 17\n",
    "mu_init = 0\n",
    "mu_step = 0.01\n",
    "mu_all = np.arange(0, npts)*mu_step + mu_init\n",
    "res = np.zeros(npts)\n",
    "\n",
    "#for i in np.arange(0, npts):\n",
    "#    integrand=DOS*fermi(energy,mu_all[i],kT)\n",
    "#    res[i]=(np.trapz(integrand,energy)-numelect)/volume_cm #in cm-3\n",
    "#\n",
    "#allres = np.column_stack((mu_all,res))\n",
    "\n",
    "integrand = DOS * fermi(energy, mu_all[:, None], kT)\n",
    "res = (np.trapz(integrand, energy, axis=1)  - numelect)/ volume_cm\n",
    "allres = np.column_stack((mu_all,res))\n",
    "\n",
    "#print(allres - allres_baseline)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.1, 0.2, 0.3, 0.4, 0.5])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.02689414, 0.02384058, 0.01422776, 0.00719448, 0.00334643])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fermi(energy, mu_all[0], kT) * DOS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.26894142, 0.11920292, 0.04742587, 0.01798621, 0.00669285],\n",
       "       [0.27091208, 0.12025686, 0.04787969, 0.01816369, 0.00675966],\n",
       "       [0.27289178, 0.12131884, 0.04833763, 0.01834289, 0.00682713],\n",
       "       [0.2748805 , 0.12238889, 0.04879972, 0.01852382, 0.00689527],\n",
       "       [0.27687819, 0.12346705, 0.04926601, 0.01870651, 0.00696409],\n",
       "       [0.27888482, 0.12455336, 0.04973651, 0.01889096, 0.00703359],\n",
       "       [0.28090034, 0.12564786, 0.05021127, 0.0190772 , 0.00710377]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fermi(energy, mu_all[:, None], kT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.02689414, 0.02384058, 0.01422776, 0.00719448, 0.00334643],\n",
       "       [0.02709121, 0.02405137, 0.01436391, 0.00726548, 0.00337983],\n",
       "       [0.02728918, 0.02426377, 0.01450129, 0.00733716, 0.00341357],\n",
       "       [0.02748805, 0.02447778, 0.01463992, 0.00740953, 0.00344764],\n",
       "       [0.02768782, 0.02469341, 0.0147798 , 0.0074826 , 0.00348204],\n",
       "       [0.02788848, 0.02491067, 0.01492095, 0.00755638, 0.00351679],\n",
       "       [0.02809003, 0.02512957, 0.01506338, 0.00763088, 0.00355189]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DOS * fermi(energy, mu_all[:, None], kT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.06038311414210114,\n",
       " 0.06091627487011181,\n",
       " 0.061453584761894324,\n",
       " 0.061995067643239914,\n",
       " 0.06254074735821567,\n",
       " 0.06309064776747156,\n",
       " 0.06364479274653036]"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[np.trapz(fermi(energy, mu_all[_], kT) * DOS) for _ in range(7)]"
   ]
  },
  {
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}