2015-05 Effective Mass Polaron Calculator.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0e-10"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Constants\n",
    "ε_0 = 8.854E-12 #Units: C2N−1m−2, permittivity of free space\n",
    "q = eV = 1.60217657E-19\n",
    "c = 3E8 # Speed of light in this universe\n",
    "m = 9.10938291E-31 # Mass of electron (kg)\n",
    "\n",
    "hbar= 1.054E-34\n",
    "h=    6.62606957E-34\n",
    "\n",
    "# Units\n",
    "Å=1E-10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "6.2e-10"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Material specifics\n",
    "d=6.2Å # lattice spacing, MAPbI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      " Cross check vs. literature values.\n",
      " NaCl Frohlich paper α="
     ]
    },
    {
     "data": {
      "text/plain": [
       "2.395939683378253"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "5.26671687867157\n",
      " CdTe  α=0.350834424140526\n",
      " GaAs  α=0.07087834402425963\n",
      "\n",
      " Guess at PCBM: 4.0, 6.0 ; α=33.826242488401434\n",
      " MAPI:\n",
      " MAPI  4.5, 24.1, 9THz ; α=1.1979698416891265\n",
      " MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=2.395939683378253\n",
      " MAPI  6.0, 25.7, 9THz ; α=0.8468397250537202\n",
      " MAPI  6.0, 36.0, 9THz ; α=0.9206337112470648\n",
      " MAPI  6.0, 36.0, 1THz ; α=2.7619011337411945"
     ]
    }
   ],
   "source": [
    "# via Feynman 1955\n",
    "#   http://dx.doi.org/10.1103/PhysRev.97.660\n",
    "function feynmanalpha(ε_Inf,ε_S,freq,m_eff)\n",
    "    ω=freq*2*pi #frequency to angular velocity\n",
    "    # Note to self - you've introduced the 4*pi factor into the dielectric constant; \n",
    "    # This gives numeric agreement with literature values, but I'm uncertain of the justification.\n",
    "    # Such a factor also seems to be necessary for calculation of Exciton binding energies (see below).\n",
    "    # Maybe scientists in the 50s assumed that the Epsilon subsumed the 4*pi factor?\n",
    "    α=0.5* (1/ε_Inf - 1/ε_S)/(4*pi*ε_0) * (q^2/(hbar*ω)) * (2*m*m_eff*ω/hbar)^0.5\n",
    "    return (α)\n",
    "end\n",
    "\n",
    "print(\"\\n Cross check vs. literature values.\")\n",
    "print(\"\\n NaCl Frohlich paper α=\",feynmanalpha(2.3, 5.6, (4.9E13/(2*pi)), 1.0)) # should be ~about 5 (Feynman1955)\n",
    "print(\"\\n CdTe  α=\",feynmanalpha(7.1,   10.4,  5.08E12, 0.095)) #Stone 0.39 / Devreese 0.29\n",
    "print(\"\\n GaAs  α=\",feynmanalpha(10.89, 12.9,  8.46E12, 0.063)) # Devreese 0.068\n",
    "\n",
    "print(\"\\n\")\n",
    "print(\"\\n Guess at PCBM: 4.0, 6.0 ; α=\",feynmanalpha(4.0,6.0, 1E12, 50))\n",
    "print(\"\\n MAPI:\")\n",
    "print(\"\\n MAPI  4.5, 24.1, 9THz ; α=\",feynmanalpha(4.5,   24.1,  9.0E12,    0.12))\n",
    "print(\"\\n MAPI  4.5, 24.1, 2.25THz - 75 cm^-1 ; α=\",feynmanalpha(4.5,   24.1,  2.25E12,    0.12))\n",
    "print(\"\\n MAPI  6.0, 25.7, 9THz ; α=\",feynmanalpha(6.0,   25.7,  9.0E12,    0.12))\n",
    "print(\"\\n MAPI  6.0, 36.0, 9THz ; α=\",feynmanalpha(6.0,   36,  9.0E12,    0.12))\n",
    "print(\"\\n MAPI  6.0, 36.0, 1THz ; α=\",feynmanalpha(6.0,   36,  1.0E12,    0.12))\n",
    "\n",
    "\n",
    "m_eff = 0.12\n",
    "#This value used for further calculations...\n",
    "α=feynmanalpha(4.5, 24.1 ,  2.25E12,    m_eff)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bare mass (no polaronic...) 0.12\n",
      "Effective mass Large polaron limit, *(1+α/6) , 0.16791879366756507\n",
      "Effective mass Small polaron limit, *0.02*α^6 , 0.45401115747646864\n"
     ]
    }
   ],
   "source": [
    "#α=1.1932891264417276 \n",
    "m_eff=0.12 #from DFT calc\n",
    "\n",
    "println(\"Bare mass (no polaronic...) \",m_eff)\n",
    "# By Feynman approximation for polaron mass. See Feynman 1955 for original literature.\n",
    "#  Also Eqn (9.6.21); p. 955\n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "println(\"Effective mass Large polaron limit, *(1+α/6) , \",m_eff*(1+α/6))\n",
    "\n",
    "println(\"Effective mass Small polaron limit, *0.02*α^6 , \",m_eff*0.02*α^6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ε_Red: 0."
     ]
    },
    {
     "data": {
      "text/plain": [
       "1.1979698416891265"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18072844628861226\n",
      "m_eff: 0.16791879366756507\n",
      "Large polaron radius rp = 2.735992793758147e-9 m\n",
      "Large polaron radius rp = 27.35992793758147 Å = 4.412891602835721 Lattice units\n",
      "Large Polaron binding energy: 4.7877648241119226e-20 J 0.2988287878976986 eV\n"
     ]
    }
   ],
   "source": [
    "# Large polaron radius via p.942 in Theoretical Solid State Physics, Vol 2\n",
    "ε_Inf=4.5\n",
    "ε_S=24.1\n",
    "ε_Red=(1/ε_Inf-1/ε_S)\n",
    "print(\"ε_Red: \",ε_Red,\"\\n\")\n",
    "\n",
    "m_eff=0.12 #from DFT calc, MAPI\n",
    "#m_eff=50 # Guess at PCBM?\n",
    "\n",
    "# Not sure if you use dressed polaron effective mass here; or direct band structure one\n",
    "m_eff=m_eff*(1+α/6)\n",
    "print(\"m_eff: \",m_eff,\"\\n\")\n",
    "\n",
    "# Eqn (9.3.7); p.942; \n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "rp=2*pi^2*hbar^2*ε_0*ε_Red^-1/(q*q*m*m_eff)\n",
    "\n",
    "println(\"Large polaron radius rp = \",rp,\" m\")\n",
    "println(\"Large polaron radius rp = \",rp/Å,\" Å = \",rp/d,\" Lattice units\")\n",
    "\n",
    "# Eqn (9.3.8); p.942; \n",
    "#  Jones and March, Theoretical solid state physics, Vol 2\n",
    "W=(0.25*q^2/rp)*ε_Red/ε_0\n",
    "println(\"Large Polaron binding energy: $W J \",W/q,\" eV\")\n",
    "\n",
    "feynmanalpha(4.5,   24.1,  9.0E12,    0.12)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Plasma frequency of free carriers screened by dielectric: ωp = 1"
     ]
    },
    {
     "data": {
      "text/plain": [
       "4.5"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ".0011502323062496e-8"
     ]
    }
   ],
   "source": [
    "Nc=10^22\n",
    "ε_Inf=4.5\n",
    "m_eff=0.12 #from DFT calc, MAPI\n",
    "\n",
    "ωp=4*pi*Nc*q^2 /(4*pi*m_eff*ε_Inf*ε_0)\n",
    "\n",
    "print(\"Plasma frequency of free carriers screened by dielectric: ωp = \",ωp )\n",
    "\n",
    "γ=1\n",
    "ω=1\n",
    "εp=ε_Inf*(1-ωp^2/(ω^2+γ^2))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Effective masses: me: 1.0931259492e-31 mh: 1.3664074365e-31 mu: 6.07292194e-32\n",
      "\n",
      "Exciton radius a0: 3.5679985202245875e-9 m 35.679985202245874 Å 5.754836322942883 Lattice Spacings\n",
      "\n",
      "n: 1 k: 0 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 0 En(k,n): -11.210728328591374 meV\n",
      "n: 3 k: 0 En(k,n): -4.982545923818388 meV\n",
      "n: 4 k: 0 En(k,n): -2.8026820821478435 meV\n",
      "n: 5 k: 0 En(k,n): -1.7937165325746196 meV\n",
      "n: 6 k: 0 En(k,n): -1.245636480954597 meV\n",
      "n: 7 k: 0 En(k,n): -0.9151614962115404 meV\n",
      "n: 8 k: 0 En(k,n): -0.7006705205369609 meV\n",
      "n: 9 k: 0 En(k,n): -0.5536162137575987 meV\n",
      "n: 10 k: 0 En(k,n): -0.4484291331436549 meV\n",
      "\n",
      "n: 1 k: 1 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 1 En(k,n): -11.210728328591374 meV\n",
      "n: 3 k: 1 En(k,n): -4.982545923818388 meV\n",
      "n: 4 k: 1 En(k,n): -2.8026820821478435 meV\n",
      "n: 5 k: 1 En(k,n): -1.7937165325746196 meV\n",
      "n: 6 k: 1 En(k,n): -1.2456364809545968 meV\n",
      "n: 7 k: 1 En(k,n): -0.9151614962115403 meV\n",
      "n: 8 k: 1 En(k,n): -0.7006705205369608 meV\n",
      "n: 9 k: 1 En(k,n): -0.5536162137575985 meV\n",
      "n: 10 k: 1 En(k,n): -0.4484291331436548 meV\n",
      "\n",
      "n: 1 k: 2 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 2 En(k,n): -11.210728328591374 meV\n",
      "n: 3 k: 2 En(k,n): -4.982545923818387 meV\n",
      "n: 4 k: 2 En(k,n): -2.802682082147843 meV\n",
      "n: 5 k: 2 En(k,n): -1.7937165325746192 meV\n",
      "n: 6 k: 2 En(k,n): -1.2456364809545963 meV\n",
      "n: 7 k: 2 En(k,n): -0.9151614962115399 meV\n",
      "n: 8 k: 2 En(k,n): -0.7006705205369603 meV\n",
      "n: 9 k: 2 En(k,n): -0.5536162137575981 meV\n",
      "n: 10 k: 2 En(k,n): -0.44842913314365435 meV\n",
      "\n",
      "n: 1 k: 3 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 3 En(k,n): -11.210728328591372 meV\n",
      "n: 3 k: 3 En(k,n): -4.982545923818386 meV\n",
      "n: 4 k: 3 En(k,n): -2.802682082147842 meV\n",
      "n: 5 k: 3 En(k,n): -1.7937165325746185 meV\n",
      "n: 6 k: 3 En(k,n): -1.2456364809545957 meV\n",
      "n: 7 k: 3 En(k,n): -0.9151614962115392 meV\n",
      "n: 8 k: 3 En(k,n): -0.7006705205369596 meV\n",
      "n: 9 k: 3 En(k,n): -0.5536162137575974 meV\n",
      "n: 10 k: 3 En(k,n): -0.4484291331436536 meV\n",
      "\n",
      "n: 1 k: 4 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 4 En(k,n): -11.210728328591372 meV\n",
      "n: 3 k: 4 En(k,n): -4.982545923818385 meV\n",
      "n: 4 k: 4 En(k,n): -2.8026820821478413 meV\n",
      "n: 5 k: 4 En(k,n): -1.7937165325746174 meV\n",
      "n: 6 k: 4 En(k,n): -1.2456364809545948 meV\n",
      "n: 7 k: 4 En(k,n): -0.9151614962115383 meV\n",
      "n: 8 k: 4 En(k,n): -0.7006705205369587 meV\n",
      "n: 9 k: 4 En(k,n): -0.5536162137575964 meV\n",
      "n: 10 k: 4 En(k,n): -0.4484291331436527 meV\n",
      "\n",
      "n: 1 k: 5 En(k,n): -44.8429133143655 meV\n",
      "n: 2 k: 5 En(k,n): -11.210728328591369 meV\n",
      "n: 3 k: 5 En(k,n): -4.9825459238183845 meV\n",
      "n: 4 k: 5 En(k,n): -2.80268208214784 meV\n",
      "n: 5 k: 5 En(k,n): -1.793716532574616 meV\n",
      "n: 6 k: 5 En(k,n): -1.2456364809545932 meV\n",
      "n: 7 k: 5 En(k,n): -0.915161496211537 meV\n",
      "n: 8 k: 5 En(k,n): -0.7006705205369573 meV\n",
      "n: 9 k: 5 En(k,n): -0.5536162137575952 meV\n",
      "n: 10 k: 5 En(k,n): -0.4484291331436514 meV\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# EXCITONS\n",
    "# Following p.973 Jones & March, Theoretical Solid State Phys Vol 2.\n",
    "# Hygronic like Wannier exciton; Eqn (9.8.10); p.973\n",
    "\n",
    "# Physical Constants\n",
    "hbar= 1.054E-34\n",
    "q = eV = 1.60217657E-19\n",
    "melectron = 9.10938291E-31 # Mass of electron, kg\n",
    "Å=1E-10\n",
    "ε_0 = 8.854E-12 #Units: C2N−1m−2\n",
    "\n",
    "ε_Inf=4.5\n",
    "ε_S=24.1\n",
    "\n",
    "ε=ε_Inf * ε_0 * 4*pi \n",
    "# term of 4pi appears because in the derivation of Dresselhaus 1956, \n",
    "# the Epsilon is used as a general prefactor in the denominator of a Coulomb term in the Hamiltonian\n",
    "\n",
    "me=0.12 * melectron\n",
    "mh=0.15 * melectron\n",
    "\n",
    "# reduced mass\n",
    "mu=1/(1/me + 1/mh)\n",
    "\n",
    "# From http:/dx.doi.org/10.1103/PhysRev.116.573 - Thomas and Hopfield, PhysRev 1959 CdS\n",
    "# Eb=28 meV; Dielectric = 9.3; Mu = 0.18\n",
    "# mu=0.18 * melectron\n",
    "# ε_Inf=9.3\n",
    "# ε=ε_Inf * ε_0 * 4*pi \n",
    "# ---> n: 1 k: 0 En(k,n): -28.34762730330596 meV\n",
    "#   == Happy Jarvist 8^)\n",
    "\n",
    "println(\"Effective masses: me: $me mh: $mh mu: $mu\\n\")\n",
    "\n",
    "En(k,n) = -mu*q^4/(2*hbar^2*ε^2*n^2) + hbar^2*k^2/(2*(me+mh))\n",
    "a0 = (ε*hbar^2)/(mu*q^2)\n",
    "d=6.2Å # lattice spacing\n",
    "println(\"Exciton radius a0: \",a0,\" m \",a0/Å,\" Å \",a0/d,\" Lattice Spacings\\n\")\n",
    "\n",
    "\n",
    "for k in 0:5\n",
    "    for n in 1:10\n",
    "        println(\"n: \",n,\" k: \",k,\" En(k,n): \",1000*En(k,n)/eV,\" meV\")\n",
    "    end\n",
    "    println()\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: replacing module PyPlot\n",
      "WARNING: Method definition show() in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:450 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:450.\n",
      "WARNING: Method definition connect(Union{Symbol, AbstractString}, Function) in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:469 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:469.\n",
      "WARNING: Method definition step(Any, Any) in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:464 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:464.\n",
      "WARNING: Method definition step(Array, Any, Any) in module PyPlot overwritten in module PyPlot.\n",
      "WARNING: Method definition close() in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:467 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:467.\n",
      "WARNING: New definition \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466\n",
      "is ambiguous with: \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466.\n",
      "To fix, define \n",
      "    close(Union{Integer, Symbol, AbstractString})\n",
      "before the new definition.\n",
      "WARNING: New definition \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466\n",
      "is ambiguous with: \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466.\n",
      "To fix, define \n",
      "    close(Union{Integer, Symbol, AbstractString})\n",
      "before the new definition.\n",
      "WARNING: New definition \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466\n",
      "is ambiguous with: \n",
      "    close(Union{Integer, PyPlot.Figure, Symbol, AbstractString}) at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:466.\n",
      "To fix, define \n",
      "    close(Union{Integer, Symbol, AbstractString})\n",
      "before the new definition.\n",
      "WARNING: Method definition fill(AbstractArray, AbstractArray, Any...) in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:471 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:471.\n",
      "WARNING: Method definition fill(Array, AbstractArray, AbstractArray, Any...) in module PyPlot overwritten in module PyPlot.\n",
      "WARNING: Method definition pygui(Bool) in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:285 overwritten in module PyPlot at /Users/jarvist/.julia/Pyplot/src/Pyplot.jl:285.\n",
      "WARNING: requiring \"Pyplot\" in module \"Main\" did not define a corresponding module.\n"
     ]
    }
   ],
   "source": [
    "using Pyplot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "LoadError: MethodError: `*` has no method matching *(::Array{Int64,1}, ::Array{Int64,1})\nClosest candidates are:\n  *(::Any, ::Any, !Matched::Any, !Matched::Any...)\n  *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(!Matched::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{DenseArray{S,1},SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n  *{TA,TB}(!Matched::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{DenseArray{TB,1},DenseArray{TB,2},SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n  ...\nwhile loading In[30], in expression starting on line 4",
     "output_type": "error",
     "traceback": [
      "LoadError: MethodError: `*` has no method matching *(::Array{Int64,1}, ::Array{Int64,1})\nClosest candidates are:\n  *(::Any, ::Any, !Matched::Any, !Matched::Any...)\n  *{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},S}(!Matched::Union{DenseArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2},SubArray{T<:Union{Complex{Float32},Complex{Float64},Float32,Float64},2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}}, ::Union{DenseArray{S,1},SubArray{S,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n  *{TA,TB}(!Matched::Base.LinAlg.AbstractTriangular{TA,S<:AbstractArray{T,2}}, ::Union{DenseArray{TB,1},DenseArray{TB,2},SubArray{TB,1,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD},SubArray{TB,2,A<:DenseArray{T,N},I<:Tuple{Vararg{Union{Colon,Int64,Range{Int64}}}},LD}})\n  ...\nwhile loading In[30], in expression starting on line 4",
      "",
      " in power_by_squaring at intfuncs.jl:78",
      " in En at In[28]:35"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ",1,2,3,4,5,6,7,8,9,10]\n"
     ]
    }
   ],
   "source": [
    "k=0\n",
    "n=[0:10]\n",
    "println(n)\n",
    "E=En(n,k)\n",
    "plot(n,En(n,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "r: 0.0 H(r) -Inf J -Inf eV\n",
      "r: 5.263157894736843e-10 H(r) -4.383539981795058e-19 J -2.7359905667544857 eV\n",
      "r: 1.0526315789473685e-9 H(r) -2.191769990897529e-19 J -1.3679952833772429 eV\n",
      "r: 1.578947368421053e-9 H(r) -1.461179993931686e-19 J -0.9119968555848286 eV\n",
      "r: 2.105263157894737e-9 H(r) -1.0958849954487645e-19 J -0.6839976416886214 eV\n",
      "r: 2.631578947368421e-9 H(r) -8.767079963590117e-20 J -0.5471981133508972 eV\n",
      "r: 3.157894736842106e-9 H(r) -7.30589996965843e-20 J -0.4559984277924143 eV\n",
      "r: 3.6842105263157898e-9 H(r) -6.262199973992941e-20 J -0.3908557952506409 eV\n",
      "r: 4.210526315789474e-9 H(r) -5.479424977243822e-20 J -0.3419988208443107 eV\n",
      "r: 4.736842105263158e-9 H(r) -4.870599979772287e-20 J -0.3039989518616096 eV\n",
      "r: 5.263157894736842e-9 H(r) -4.3835399817950583e-20 J -0.2735990566754486 eV\n",
      "r: 5.789473684210527e-9 H(r) -3.9850363470864167e-20 J -0.24872641515949873 eV\n",
      "r: 6.315789473684212e-9 H(r) -3.652949984829215e-20 J -0.22799921389620714 eV\n",
      "r: 6.842105263157895e-9 H(r) -3.371953832150045e-20 J -0.21046081282726817 eV\n",
      "r: 7.3684210526315796e-9 H(r) -3.1310999869964703e-20 J -0.19542789762532045 eV\n",
      "r: 7.894736842105263e-9 H(r) -2.9223599878633724e-20 J -0.18239937111696575 eV\n",
      "r: 8.421052631578948e-9 H(r) -2.739712488621911e-20 J -0.17099941042215536 eV\n",
      "r: 8.947368421052632e-9 H(r) -2.5785529304676813e-20 J -0.1609406215737933 eV\n",
      "r: 9.473684210526316e-9 H(r) -2.4352999898861435e-20 J -0.1519994759308048 eV\n",
      "r: 1.0e-8 H(r) -2.3071263062079254e-20 J -0.143999503513394 eV\n"
     ]
    }
   ],
   "source": [
    "V(r)=-q^2/(4*pi*ε_0*r) #- q^2/()\n",
    "\n",
    "for r in linspace(0,100*Å,20)\n",
    "    println(\"r: $r H(r) \",V(r),\" J \",V(r)/eV,\" eV\")\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: [a] concatenation is deprecated; use collect(a) instead\n"
     ]
    },
    {
     "ename": "LoadError",
     "evalue": "LoadError: UndefVarError: plot not defined\nwhile loading In[32], in expression starting on line 21",
     "output_type": "error",
     "traceback": [
      "LoadError: UndefVarError: plot not defined\nwhile loading In[32], in expression starting on line 21",
      "",
      " [inlined code] from In[32]:22",
      " in anonymous at no file:0"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      " in depwarn at deprecated.jl:73\n",
      " in oldstyle_vcat_warning at /Applications/Julia-0.4.2.app/Contents/Resources/julia/lib/julia/sys.dylib\n",
      " in anonymous at no file\n",
      " in include_string at loading.jl:266\n",
      " in execute_request_0x535c5df2 at /Users/jarvist/.julia/IJulia/src/execute_request.jl:177\n",
      " in eventloop at /Users/jarvist/.julia/IJulia/src/IJulia.jl:141\n",
      " in anonymous at task.jl:447\n",
      "while loading In[32], in expression starting on line 21\n"
     ]
    }
   ],
   "source": [
    "# Following: http://www.cond-mat.de/teaching/DFT/qm1d.html\n",
    "# 1D numeric Shrodinger equation solver, by discretisation to an Eigenvalue problem\n",
    "#  $f''(x_i) \\approx (f(x_{i-1})-2*f(x_i)+f(x_{i+1}))/dx^2$\n",
    "\n",
    "N=99\n",
    "\n",
    "#V(r)=-3/abs(r) #- q^2/()   #... ish\n",
    "#V(r)=1E-3*r^2  # Harmonic well\n",
    "V(r)= 0.01*r^2+0.1*r  # an-Harmonic well\n",
    "\n",
    "dx=20/(N-1)\n",
    "\n",
    "# PE terms on the trace\n",
    "diagonal = [(2.0/dx^2 + V(r))::Float64 for r in -N/2:N/2]\n",
    "# KE terms on the tridiagonals\n",
    "updiagonal = [(-1/dx^2)::Float64 for r in 1:N]\n",
    "\n",
    "H =diagm(diagonal,0) + diagm(updiagonal,1) + diagm(updiagonal,-1)\n",
    "evals,evec=eig(H)\n",
    "\n",
    "for i in [1:7]\n",
    "    plot(evec[:,i]+evals[i])\n",
    "    plot(evals[i])\n",
    "end\n",
    "    \n",
    "plot([V(r) for r in -N/2:N/2],color=\"black\") # Potential energy curve\n",
    "println(evals[1:10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.4.2",
   "language": "julia",
   "name": "julia-0.4"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.4.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}