Wire 1D + diatomica TD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cable Unidimensional con diatómica (caso TD)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importamos los módulos que necesitamos, para dibujar y hacer álgebra matricial"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import math\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trabajamos en unidades atómicas pero nos gusta que aparezcan las constantes en los lugares que esperamos, hacemos $\\hbar =1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hbar = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parámetros generales de la simulación que definen el sistema:\n",
    "* `nsites` es la longitud del cable en sitios, lo hacemos par por conveniencia, esto nos permite poner la diatómica en el medio.\n",
    "* `temp` es la temperatura en unidades de $k_b T$, la usamos para construír la matriz densidad inicial.\n",
    "* `mu` es el potencial electroquímico de los electrones para el sistema en equilibrio.\n",
    "* `beta` es el hopping entre sitios en el cable $\\beta$\n",
    "* `gamma` es el hopping de contacto entre la diatómica y el lead $\\gamma$\n",
    "* `delta` es el hoppoing en la diatómica, la separación HOMO-LUMO es igual a $|2\\delta|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsites = 502\n",
    "temp = 0.01\n",
    "mu = 0.0\n",
    "beta = -1.0\n",
    "delta = -0.25\n",
    "gamma = -0.10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parámetros generales que definen **esta** simulación:\n",
    "* `nsteps` es el número de pasos que queremos correr.\n",
    "* `dt` es el paso de tiempo $\\Delta t$.\n",
    "* `DeltaV` es la diferencia de potencial inicial (bias) $\\Delta V$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsteps = 2100\n",
    "dt = 0.1\n",
    "DeltaV = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definimos la $\\delta$ de Kronecker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def kdelta(k, l):\n",
    "    if k == l: return True\n",
    "    else: return False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construímos el Hamiltoniano con el que armaremos el estado inicial del sistema, es tridiagonal, y tiene incorporado el bias ($\\Delta V$), a la izquierda subimos el potencial local de los sitios en $\\frac{\\Delta V}{2}$ y a la derecha lo bajamos en $\\frac{\\Delta V}{2}$. Dejamos nuestra \"impureza\" sin potencial. **IMPORTANTE** Esta forma de armar el Hamiltoniano no está escrita para que sea eficiente (y mucho menos en Python) sino para que sea fácil de leer. Nótese por ejemplo que recorremos **toda** la matriz para hacer sólo unos muy poquitos elementos distintos de cero. En la \"vida real\" uno haría las cosas de otra forma mejor, hay, como siempre, muchas posible y no entraremos en esto."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "h = np.zeros((nsites,nsites))\n",
    "\n",
    "for i in range(0,nsites):\n",
    "    for j in range (0,nsites):\n",
    "        if kdelta(i,j-1):\n",
    "            h[i,j] = beta\n",
    "        elif kdelta(i,j+1):\n",
    "            h[i,j] = beta\n",
    "        elif kdelta(i,j) and i < int(nsites/2) - 1:\n",
    "            h[i,j] = DeltaV / 2\n",
    "        elif kdelta(i,j) and i > int(nsites/2) + 0:\n",
    "            h[i,j] = - DeltaV / 2\n",
    "            \n",
    "h[int(nsites/2)-2,int(nsites/2)-1] = gamma\n",
    "h[int(nsites/2)-1,int(nsites/2)  ] = delta\n",
    "h[int(nsites/2)  ,int(nsites/2)+1] = gamma\n",
    "\n",
    "h[int(nsites/2)-1,int(nsites/2)-2] = gamma\n",
    "h[int(nsites/2)  ,int(nsites/2)-1] = delta\n",
    "h[int(nsites/2)+1,int(nsites/2)  ] = gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "h[int(nsites/2)-5:int(nsites/2)+5,int(nsites/2)-5:int(nsites/2)+5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Definimos la función de Fermi \n",
    "$$f(\\epsilon) = \\frac{1}{1-e^{\\frac{\\epsilon-\\mu}{k_b T}}}$$\n",
    "nótese que revisamos el valor del exponente para no caer en una situación de under/overflow antes de evaluar la exponencial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def fermif(energy, mu, temp):\n",
    "    expo = (energy-mu)/temp\n",
    "    if expo < -100.0:\n",
    "        return 1.0\n",
    "    elif expo > 100.0:\n",
    "        return 0.0\n",
    "    else:\n",
    "        return 1.0/(math.exp(expo)+1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Construímos la matriz densidad del estado inicial $\\hat\\rho_0$. Para eso primero la armamos en la base de autofunciones de $\\hat{H}$, en la que es diagonal y vale $\\hat\\rho_0^d = f(\\hat H)$, con $f$ la función de Fermi. Para eso diagonalizamos $\\hat H$, `w` contiene los autovalores y `V` la matriz de autovectores. Usamos le aplicamos $f$ a `w` para armar la $\\hat\\rho_0$ en la base de autofunciones de $\\hat H$ y después la transformamos a la base de sitios haciendo $\\hat\\rho_0 = V \\hat\\rho_0^d V^\\dagger$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "w, v = np.linalg.eig(h.real)\n",
    "rho_diag = np.zeros((nsites, nsites))\n",
    "for i in range(0, nsites):\n",
    "    rho_diag[i, i] = fermif(w[i], mu, temp)\n",
    "rho_0 = np.dot(v, np.dot(rho_diag, np.transpose(v)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Miremos la pinta de las ocupaciones"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.scatter(w,np.diag(rho_diag))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Declaro una serie de arreglos vacíos que contendrán cosas que pretendo ir guardando durante la simulación, el tiempo, la traza, energía, corriente, etc..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "time = []\n",
    "trace = []\n",
    "energy = []\n",
    "current = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Integración\n",
    "\n",
    "El algoritmo de integración que usamos es el siguiente:\n",
    "$$ \\rho(t+\\Delta t)=\\rho(t-\\Delta t) + 2 \\Delta t\\ \\dot\\rho(t)$$\n",
    "con\n",
    "$$ \\dot\\rho(t) = -\\frac{i}{\\hbar} [H,\\rho] $$\n",
    "el algoritmo es un *leapfrog* que usa info del paso anterior y la derivada en el presente para encontrar la $\\rho$ futura. Para inicializarlo hacemos un paso de Euler hacia atrás\n",
    "$$\\rho_{\\mathrm old} = \\rho_0 - \\Delta t\\ \\dot\\rho(t=0)$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "rho = rho_0\n",
    "rhodot = ((0-1.0j)/hbar)*(np.dot(h,rho)-np.dot(rho,h))\n",
    "rhoold = rho - dt*rhodot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hasta ahora nuestra $\\rho$ conmuta (por construcción) con el $H$ que armamos con bias, para lograr que evolucione armamos un nuevo Hamiltoniano eliminando el bias. Con esto tenemos un cable con un exceso de electrones del lado izquierdo y un defecto en el lado derecho.Con este nuevo Hamiltoniano deberíamos generar una corriente hacia la derecha. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.fill_diagonal(h,0.0)\n",
    "\n",
    "h[int(nsites/2)-2,int(nsites/2)-1] = gamma\n",
    "h[int(nsites/2)-1,int(nsites/2)  ] = delta\n",
    "h[int(nsites/2)  ,int(nsites/2)+1] = gamma\n",
    "\n",
    "h[int(nsites/2)-1,int(nsites/2)-2] = gamma\n",
    "h[int(nsites/2)  ,int(nsites/2)-1] = delta\n",
    "h[int(nsites/2)+1,int(nsites/2)  ] = gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hacemos evolucionar $\\rho$ por `nsteps` guardando a cada paso la información que nos interese. De nuevo este código no es la forma más elegante ni eficiente, está escrito para que sea fácil de entender."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for step in range(0,nsteps):\n",
    "    rhodot = ((0-1.0j)/hbar)*(np.dot(h,rho)-np.dot(rho,h))\n",
    "    rhonew = rhoold + 2.0 * dt * rhodot\n",
    "    rhoold = rho\n",
    "    rho = rhonew\n",
    "    # calculate observables\n",
    "    time.append(step*dt)\n",
    "    trace.append(2.0*np.trace(rho).real)\n",
    "    energy.append(np.trace(np.dot(h,rho)).real)\n",
    "    current.append(4.0 * h[int(nsites/2)-1,int(nsites/2)].real * rho[int(nsites/2),int(nsites/2)-1].imag )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Nuestro algoritmo conserva la traza de forma numéricamente exacta."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time,trace)\n",
    "ax = plt.gca()\n",
    "ax.ticklabel_format(useOffset=False)\n",
    "plt.ylim((2.0*np.trace(rho_0)-0.000001,2.0*np.trace(rho_0)+0.000001))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y también debería conservar la energía."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time,energy)\n",
    "ax = plt.gca()\n",
    "ax.ticklabel_format(useOffset=False)\n",
    "plt.ylim(np.trace(np.dot(h,rho)).real-0.000001,np.trace(np.dot(h,rho)).real+0.000001)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Y acá está ahora lo más interesante de todo! Qué pasa con la corriente?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(time,current)\n",
    "ax = plt.gca()\n",
    "ax.ticklabel_format(useOffset=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}