Conductancia_Diatomica_IvsV_TD.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Importamos todos los paquetes que necesitamos"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "hbar = 1.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsites = 1002\n",
    "temp = 0.01\n",
    "mu = 0.0\n",
    "beta = -1.0\n",
    "delta = -0.25\n",
    "gamma = -0.10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nsteps = 2100\n",
    "dt = 0.1\n",
    "DeltaV = 0.5"
   ]
  },
  {
   "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": "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": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "DeltaVs = np.linspace(0.0,1.0,20)\n",
    "\n",
    "currents = []\n",
    "\n",
    "for ii, DeltaV in enumerate(DeltaVs):\n",
    "    \n",
    "    # Print something to let us know it's alive\n",
    "    \n",
    "    print('Starting with ',ii,' DeltaV = ',DeltaV)\n",
    "    \n",
    "    # Build Biased H\n",
    "    \n",
    "    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\n",
    "    \n",
    "    # Build biased GS DM\n",
    "    \n",
    "    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)))\n",
    "    \n",
    "    # placeholders for output\n",
    "    \n",
    "    time = []\n",
    "    current = []\n",
    "    \n",
    "    # Boot integrator\n",
    "    \n",
    "    rho = rho_0\n",
    "    rhodot = ((0-1.0j)/hbar)*(np.dot(h,rho)-np.dot(rho,h))\n",
    "    rhoold = rho - dt*rhodot\n",
    "    \n",
    "    # Build propagation H\n",
    "    \n",
    "    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\n",
    "    \n",
    "    # Do time propagation for nsteps\n",
    "    \n",
    "    for step in range(0,nsteps):\n",
    "        if (step % 1000) == 0:\n",
    "            print('Integration step ',step,' of ',nsteps,' in DeltaV ',DeltaV)\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",
    "        current.append(4.0 * h[int(nsites/2)-1,int(nsites/2)].real * rho[int(nsites/2),int(nsites/2)-1].imag)\n",
    "        \n",
    "    # Save current trace for this Delta V\n",
    "    \n",
    "    currents.append(np.array(current))\n",
    "    \n",
    "    # Print something to let us know it's alive\n",
    "    \n",
    "    print('Done with ',i,' DeltaV = ',DeltaV)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "for i,DeltaV in enumerate(DeltaVs):\n",
    "    plt.plot(time,currents[i])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ivsv = []\n",
    "for i,DeltaV in enumerate(DeltaVs):\n",
    "    ivsv.append(np.average(currents[i][-1000:-1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dv_gf = np.loadtxt('deltavs.dat')\n",
    "ivsv_gf = np.loadtxt('ivsv.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(DeltaVs,ivsv)\n",
    "plt.plot(dv_gf,ivsv_gf)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dv_gf_ext = np.loadtxt('deltavs-ext.dat')\n",
    "ivsv_gf_ext = np.loadtxt('ivsv-ext.dat')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(DeltaVs,ivsv)\n",
    "plt.plot(dv_gf_ext,ivsv_gf_ext)\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
}