kkr.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/alex/dev/pymatgen/pymatgen/io/vasp/outputs.py:818: UserWarning: No POTCAR file with matching TITEL fields was found in /home/alex/dev/kkr/kramers-kronig\n",
      "  \" was found in {}\".format(os.path.abspath(p)))\n"
     ]
    }
   ],
   "source": [
    "# First load the dielectric constant data from the VASP output\n",
    "import numpy as np\n",
    "from pymatgen.io.vasp.outputs import Vasprun\n",
    "\n",
    "def symmetrise(a):\n",
    "    \"\"\"Turn a XX, YY, ZZ, XY, YZ, XZ array into a symmetrical 3x3 matrix\"\"\"\n",
    "    return [[a[0], a[3], a[5]], [a[4], a[1], a[4]], [a[5], a[4], a[2]]]\n",
    "            \n",
    "vr = Vasprun('vasprun.xml')\n",
    "energies, eps_real, eps_imag = vr.dielectric\n",
    "eps_real = np.array([symmetrise(a) for a in eps_real])\n",
    "eps_imag = np.array([symmetrise(a) for a in eps_imag])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 19.3 s per loop\n"
     ]
    }
   ],
   "source": [
    "# python version based on code in linear_optics.F\n",
    "import math\n",
    "\n",
    "def kkr(de, eps_imag, cshift=1e-6):\n",
    "    eps_imag = np.array(eps_imag)\n",
    "    nedos = eps_imag.shape[0]\n",
    "    eps_real = []\n",
    "    for i_r in range(nedos):\n",
    "        w_r = de * i_r\n",
    "        total = np.zeros(eps_imag.shape[1:], dtype=np.complex_)\n",
    "        \n",
    "        for i_i in range(nedos):\n",
    "            w_i = de * i_i\n",
    "            val = eps_imag[i_i] * ((1. / (w_r - w_i - complex(0, cshift)))\n",
    "                                   + (1. / (- w_r - w_i + complex(0, cshift)))) * (-0.5)\n",
    "            total += val\n",
    "        eps_real.append(total * (2 / math.pi) * de + np.diag([1, 1, 1]))\n",
    "    return np.real(np.array(eps_real))\n",
    "\n",
    "%timeit kkr(energies[1] - energies[0], eps_imag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1 loop, best of 3: 282 ms per loop\n"
     ]
    }
   ],
   "source": [
    "# improved numpy version with Adam Jackon's help\n",
    "import math\n",
    "import scipy.integrate\n",
    "\n",
    "def kkr(de, eps_imag, cshift=1e-6):\n",
    "    eps_imag = np.array(eps_imag)\n",
    "    nedos = eps_imag.shape[0]\n",
    "    cshift = complex(0, cshift)\n",
    "    w_i = np.arange(0, nedos*de, de, dtype=np.complex_)\n",
    "    w_i = np.reshape(w_i, (nedos, 1, 1))\n",
    "    \n",
    "    def integration_element(w_r):\n",
    "        factor = w_i / (w_i**2 - w_r**2 + cshift)\n",
    "        total = np.sum(eps_imag * factor, axis=0)\n",
    "        return total * (2/math.pi) * de + np.diag([1, 1, 1])\n",
    "\n",
    "    return np.real([integration_element(w_r) for w_r in w_i[:,0,0]])\n",
    "\n",
    "%timeit kkr(energies[1] - energies[0], eps_imag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAD8CAYAAABjAo9vAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XeYXFd9//H3mTtte9Gu2qqsZFmWLTcJNVdswGBsQov5\nBfiZ0A1JIEAgxECeXzBJHggktAdiMGBsEgwYDIQYMNi4gG1Zliyr2JLV66ps7zs75Z7fHzM7W7Vt\nZvbO7n5ez6NHM/feOfera+njs2fOPddYaxERkZnD53UBIiKSXQp2EZEZRsEuIjLDKNhFRGYYBbuI\nyAyjYBcRmWEU7CIiM4yCXURkhlGwi4jMMH4vTlpVVWVra2u9OLWIyLT13HPPNVprq8c6zpNgr62t\nZdu2bV6cWkRk2jLGHBvPcRqKERGZYRTsIiIzjIJdRGSGUbCLiMwwCnYRkRlGwS4iMsMo2EVEZhgF\ne5Yc27+TbY/c73UZIiIK9mzpuv+DrHvy/TQ3nPa6FBGZ5RTsWXJRfA8AJ1582uNKRGS2U7BnQTwW\nS7/uObPfw0pERBTsWXH62L70a9tW52ElIiIK9qxoPLIr/drfqTF2EfGWgj0LIqeT4+uHnVqKImc8\nrkZEZjsFexb4m/ZTTyUtRSsoizd4XY6IzHIK9gy5iQSrWv9IfWgpseKFVLtNuImE12WJyCymYM/Q\nsz/9EiWmBwBTupCgidPUeMrjqkRkNlOwZ6q7CYDwjZ8jWLkYgNbTR72rR0RmPQV7hiwGgGWXXElJ\n9RIAOhqOe1mSiMxyWQl2Y0y5MeZnxpiXjDF7jTFXZKPdaSEeIWodHL+fioW1AESbTnpbk4jMatl6\nmPXXgIestbcYY4JAYZbazXsmHiFKkCBQUVVDzDrYthNelyUis1jGwW6MKQOuBd4FYK2NAtFM250u\nTKKXqAkA4HMcjvqXUtK40+OqRGQ2y8ZQzDKgAfi+MeZ5Y8x3jTFFQw8yxtxmjNlmjNnW0DBz5nr7\nEr1ECabfN5WsYm5MQzEi4p1sBLsfWAvcaa1dA3QBtw89yFp7l7V2nbV2XXV1dRZOmx98iV5ipj/Y\nXX+YALFRPiEiklvZCPaTwElr7ZbU+5+RDPpZwZfoJT4g2K0/TMjOmpEoEclDGQe7tfYMcMIYc0Fq\n0yuBPZm2O104bpSYb0iwE8W6rodVichslq1ZMR8GfpiaEXMYeHeW2s17jju4x44/hGMs0ViUYCjs\nXWEiMmtlJdittTuAddloa7rxu1GizoDZnb7kDJlEPAYKdhHxgO48zZDfRkk4ofR743MASCTiXpUk\nIrOcgj1Dfhsl4esPdvqCPa4VHkXEGwr2DAVsFHdAjx2TDHY3oSmPIuINBXuGwjaC6+8fS9dQjIh4\nTcGeoULbgw0Up9/3Bbt1NRQjIt5QsGcgFu0lbGLYUEn/Rl9yopF67CLiFQV7Bro7WgEwof4ee9+X\np25cwS4i3lCwZ6C7MxnsvgE99r6hGFc9dhHxiII9A5GudgCcgrL0Nn15KiJeU7BnoLerDYBAwYAe\nu5McY7eugl1EvKFgz0CsO9ljDxaWprf1D8VoVoyIeEPBnoFYd7LHHioeOBST7LFrjF1EvKJgz0A8\n0gFAqKg8vc2nHruIeEzBngE3FeyFA3rspMbYXVdLCoiINxTsGbC9qWAvGd5jRz12EfGIgj0DJtJO\njw0OeqBGeihGSwqIiEcU7BnwRVroMMWDNzrJB23oy1MR8YqCPQOBaBtdvtJB23xO3yJgCnYR8YaC\nPQPBeDs9/pJB29KrO2qMXUQ8omDPQGG8g17/0B576s5Tqx67iHhDwZ6BQreDWLBs0La+YNc8dhHx\nioI9AyW2Ezc0ONj77jxFY+wi4hEF+yT1RropNL1QUDloe9+Xp27C9aIsEREF+2R1NDcAYArLB233\n+ZKX1FoFu4h4Q8E+SW2NdQAEy+YN3mFSd54q2EXEIwr2SepsSgZ7QWXNoO0+n0m+ULCLiEcU7JMU\nbTkFQGn14kHbjXrsIuIxBfskxdtPA1A5b9Gg7ekblFwFu4h4Q8E+Sb7Os7RTRLigaNB2oy9PRcRj\nCvZJCvQ00uKrHLa9b1aMUbCLiEcU7JNUGG2k0z882NVjFxGvKdgnqTTeSE+4eth2n0kFu8bYRcQj\nCvZJiMeizHUbiZUsGrbPpIIda6e4KhGRJAX7JNTXHcFvXJzK2uE7fX2XVD12EfGGgn0Smk/uA6Bw\n3nnD9vV9eap57CLilVkT7Nsf/hEv/fN6djz+84zb6j57GIDKmpXD9vk0j11EPJa1YDfGOMaY540x\nD2arzWzpjfZS+9QnWZXYT9Gf/iXj9hLNR4lbH3MXLR+2r+8GJfXYRcQr2eyxfwTYm8X2smbXlseo\npJ3jzhLOix+mq7M9o/YC7cep91XjDwSH7TMaihERj2Ul2I0xi4Cbge9mo71sa9n9e1wMrZe+H5+x\nnNj3XEbtlfTU0RxcMOK+0cbYrbU88+t7iPR0Z3R+EZHRZKvH/lXgk+TpVJDqxs2cCJ1P1SWvAqDt\nyI6M2psTP0N3Yc2I+0Zb3XH3n37Fpq0fYfv3/y6j84uIjCbjYDfGvA6ot9aO2g02xtxmjNlmjNnW\n0NCQ6WnH7eSZBlYn9tG+4CrmL72AbhvCnt0z6fa6O9uoohW3rHbE/X2rO9oR5rHHOhsBCHefmvT5\nRUTGko0e+1XA640xR4EfA68wxvz30IOstXdZa9dZa9dVVw+/YzNXDm37HUGToPKSV+NzHE4GllLU\ntn/S7dWfOACAv6p2xP2jjrGnwt5iJn1+EZGxZBzs1tpPWWsXWWtrgbcCj1prb824sixJHHqcCEEW\nXno9AG0l57MwemTEHvV4tNYlg71k/ooR9/dNdzQjtG9JbTMKdhHJnRk9j91ay+KWLRwrvAQTKADA\nrVrFHNpoqq+bVJuR+kMAVC+5YMT9fV+eJiIjzLxRj11EpkBWg91a+7i19nXZbDMTh48e5nyO07vk\n2vS24iWXAXDmwPOTa7T1GF02TEXVuWbFJHvsV535r+E70714BbuI5M6M7rGf2v4QAHMvvzG9bf7K\nNQB0HN81qTbDnSeod+b3j6UPYUYbZtHCYCIyBWZ0sPuPPkEbxcxbuT69bc7cxbRSgq9hcvdSVURO\n0hpeeM795wr8JI2xi0juzdhgt67L8o6tHCld13+bP4Ax1AWXUdZxcMJtdra3sDhxksic1eM6fvPn\nXzu4JvXYRWQKzNhgP7bveebRTHzpy4ft6yo9jwXx4xNvc/fT+IylcPnGcR1/Re/T5wxzLRImIrky\nY4P9zNZfAlCz/vXD9tmypZTRRVtr04TabD/0DABLL7lm/J9paxlw4mTIr2t/hB1fuGFC5xYRGa8Z\nG+yVxx/ikHMeC5YMn28erK4FoOH4vgm1GT67nZNmPuVV88f9mdaz/T8ZDOy9r4lum9C5RUTGa0YG\n+/7tT7Ayvp/6824ZcX/fzUXtpw9NqN2arj2cKbl4Qp+J9Q5c8Etj7CKSezMy2Lsf+QKdtoCLb/7g\niPvnpm4u6m04Mu426+uOMJdm4gvWTqiWeG/PgHcKdhHJvRkX7Lv/+Asu736a3cveQ0lZ5YjHlFZU\n021D0D7+xbhOvvAkAOXnb5pQPfFof7CPtMyAiEi2zahgbzxbx/xHP8Zxs5C1b/3MuQ80hhZfBYGe\ns+Nuu/fYVmLWoXb1xII9EY2kX2u6o4hMhRkT7G2tTTR/542U2k5ib/4+oXDR6McHqiiIjH/54JKm\nnRzz1xIuGL3doRJRDcWIyNSaEcF+bN8OWr5+HbWxQ+y5+uucd8nYvepIqJrS+PimO7qJBEt699NU\nPrEvTgFsIjrgjYJdRHLP73UBmTh78iBHHvwya0//hC5TyL5XfZ8117xhXJ+NFc2nuv1Jor29BEOh\nUY89eWg3S+jGV/OyCdfoxmPp1xqKEZGpMK2CffcTP6fr4FOY7kYq2vZwXuwA84xla/lrqP2LL3LJ\nwtpxtxWs3Uj4zI/Yu/tpLlx3/ajHNhx8niVAxXnrJlzzoB67hmJEZApMq2DvfuHXbGr4GW0UcTK4\nnK1L3sfi69/H+uWrJtxW2aILk202HB3z2GhD8uEaC5aPb42YgeyAHruIyFSYVsF+2bu/RjxwJ2WB\nIGUZtlValVyhMdZWP+axTssRGimnqqR8wuexiQHBrqEYEZkC0+rL03BhMf5AMCttVVQtwLUG2zl6\nsCficTa0/obGwLmX6h29AQ3FiMjUmlbBnk2OP0C7KcJERp8Zc/ZkctmBzoKaSZ1HPXYRmWqzNtgB\nuk0RTqxr1GNa6pLrtofXv2NyJxmjx779sZ+z/bGfT65tEZERTKsx9mzr8RXij3WMekxX6uHVFTXD\nV4kcFzfe/3qEHvvaJ96dfHH9myfXvojIELO6x97rFBGM9/fY47EYT939DxzctTm9LdFwkLj1UV1z\n3rjbfbbi5v43A4ZiBg3LiIjkyKwO9phTRCjRH+zP3vXXXHX8WxT//Nb0tpKmnRz1LyMYCk/qHGbg\nUMyQYNdTlEQkF2Z3sAdKCLvJYK87vIcrG+4HYD6NNNefIh6LsSzyEk0Vl02s4YFDLgOGYgbfrAQN\nZ05OrnARkVHM6mBPBIoptMkHYZzc+iAAe0LJED91YDsnDuygyETwLV4/6XMYt7+Xbob02Ou///ZJ\ntysici6zOthtsIQim1x90Tbup9uGKHnLNwHoaThK08Hk4+uqV26YULtmwOwXX7x/dcehY+zWP7nh\nHRGR0czqYCdUQsjE6I10E+48wRlnAXMXLQfAv+9B4nU76bUBFq24dGLtDhiKKe06mn5thgzFdBck\nb3pqpXhy9YuIjGBWT3c04VIAOttbKO+to6VgCctT67iv6dmM2/0Mh/3LWTHhu137vxQtSHSmXzu9\nLYOOcmLJfY5NTKJ6EZGRzeoeu1OQDPbujlbK3RaiBfMG7fcZS0do3kgfHbdC2z/rpqTr+KB9gVh7\nsg40O0ZEsmdWB7sv1WPv7Wyl0PbgBocPiUQLqifcbt+zTdsposR2Ya0l2tvLqvjeQccVxlsB8KMe\nu4hkz6wO9kDqMXedTXUETQITLAFg380/o4NCANzCiQd739IBHaY4PYa/f+vvhx1VlEj22BXsIpJN\nszvYw8kgv/yP709uCCV77Besv4FDBckvTP0VSybcbt+smG4n+RNBS3MjkdYzw44rcZPB7jMWN6Fw\nF5HsmNVfngYLSwa99wUL0q8DieT89lDFJFZ1TA3FRPylEIcF37kUy/Cef4npnwoZj8cIOs7EzyUi\nMsSs7rGHCgaPqQ+8xb/4TV/hhdAaateM/ti8kSWDPRrofxzIQhpGPLKZZK8+Ho+OuF9EZKJmdY89\nVFQ66H3p4ovSr5deuA4ufHxyDad67PFQOYy+eCTtvnIq3XbiMS0QJiLZMat77AUDhmKeL7mOCze9\nNivt9o2xu+GxH6XXGZiTPFbPRhWRLFGwp8SXvzJr7dpUjz1QNXip360X3j7s2EgoGexxBbuIZMms\nDnafv38kyvhDWWvX9L0IhGihf7inYMGqYcfGU/PkE1qrXUSyJONgN8YsNsY8ZozZY4x50RjzkWwU\nNtV8gewFe9+Xp8b46DbJmTbHnKUsWX3V8EOLUsGuMXYRyZJs9NjjwMettRcBm4C/McZcNMZn8o4v\nkL2VFs2AYI/4kjc6tV/zTzDCdEZfYQWgMXYRyZ6Mg91ae9pauz31ugPYC0xi8re3nCz22I1NTps0\nPh9RX7LHbgIhjDHDjvWFk1MuE4n4sH0iIpOR1TF2Y0wtsAbYMsK+24wx24wx2xoaRp7T7SUni2Ps\n6UV7jYM1yV66Mb6Rg90JAJBIaB67iGRH1oLdGFMMPAB81FrbPnS/tfYua+06a+266urJrL+SW75g\nNr887V+P3fZ9lWpHXsHRpILd1Ri7iGRJVoLdGBMgGeo/tNb+PBttTjV/FsfY+3LdGB/WJC+xte6w\nHnvMOumZOa6GYkQkS7IxK8YA3wP2Wmu/nHlJ3ggEs/mYunSyp3vsA5cr6NNLEONL9dg13VFEsiQb\nPfargHcArzDG7Ej9uikL7U4pf1aHYlIhbgylN32OM84CVqx5+bAee9QE8DmpHrtmxYhIlmS8Voy1\n9kkG3JMzXTnZ7LHbvumOsPSya+GylwDo6Rq8cEwcf/rLU/XYRSRbZvWdpwMFQwVjHzRO/f+XG/3y\nuvjw+ZPBbjXGLiJZomBPKZ8zP4ut9Y+xDzR0KMZi0sGuL09FJFtm9bK9AIff+L+0HtzC2qw+5KL/\nztNBhgS7a3wDxtg1j11EsmPWB/vyy6+Fy6/NapvpO09HuCFpIIsPp28oxlWPXUSyQ0MxuTTGUEyy\nx54M9u62pikrS0RmNgV7Dgy883TQ9iFDMwO/PL167+fY/Iv/zHltIjLzKdhzaegY+xAWH06qxw4Q\nOPKHXFckIrOAgj0HzjXGPtKsGCfgH/Y5EZFMKNhzYuTpjsOOMj4cfzD93ol357IoEZklFOw5YNK/\nD768w748HToUk1Cwi0jmFOy5YM91g9Lgy21N/5enAAE3kvPSRGTmU7DnQN+sGN845rEHBoyx+13d\npCQimVOw50Sqx+4bex77wDH2sO3JeWUiMvMp2HPAjPAKRpoV48M/INgL7MwaY9+z+bds/u/P0t3Z\n5nUpIrOKgj0X+pbtHWM1Y2sMPl//f4Ji2z3iAzmmk7aWRgDOHD9A7UPv5IqDX2HPXe/1uCqR2UXB\nngPmXKs7+oZebt+gbQGTINIzfXvtLzz1v5R97Tx2PPJjjv3iDhxcniu+jrVtj3B8/w6vyxOZNRTs\nOdE3xj7Geuwj3Jna1d6ci4KmRPfW+wCIHHicy5sfYkfVTSy99Rsk8FH36Lc9rk5yIdobmfY/Zc5E\ns351x1zo67GPPRQzPNi7O5phwZKc1JVrZZ2HAVh35if4jcv8V/8dVfMXs6NoA0vPPuJxdZJtR/du\no/Anb+F4yRrWfXxaPsN+4qzFugncRJxEIk4iHifhJnBjMRJucpubSJCIx7CJBIlEDOvGSSQSJOJx\nrJtg3rLVlJbPyWmZCvYcSAf7GD32kX5gaqk7yJKVl+egqtzq7mxjWewAGPAbl72B1Vx4QfLP0VN9\nGQuPbaaro5WiknKPK5VzspZIpJvWlmbaW5uINB4n1nKcRNNRXNeFcBm+wnICJVXE9vyGS5t/R9jE\nmNvxBzbf+xlW3PB+qhfWDmqypfEM+/9wDxe88l2UV438MBubiNHT2UZPZyu9na1EerqI9kaIRntJ\nxKK48V7ceAw3HsXGo7iJ5O82EYNEFBJRbDwGbgyTiIIbxySiGDeGz41h3Dg+N4bPxvDbGD6bwG9j\nODaOY+P4ieO3MfzEcWwCHy4+Eji4ydfWxUm9d0zy8fRO6tdk7Lrue1x63S2T/PT4KNhzIb2448R7\n7NG2s9mvZwq89OQvWGsS6fdtNdelX4cWXgzHoO7ATlaufbkH1c1ciWiE7s5WIh2t9HS1Eu1sIxpp\nJ97dQSKS/GV7OzHRToh24cQ6ceJd+OPdBBPdBN1uQm4PBbaHAhshbBLMBwZGsGuTf499pn/V0m4b\nYlfla3BWXM/Ltn6cK458g55v38WzlTfghsoJr7gGu/0HXND1HBtNL+z9PLtDa4j5SwBLRc8JSt1W\nCm03BSZKIVCY4bWI4RDHTzKi/SSMn7jxkzABEsZPAj8JX4CEz0/CFKT2Jbe5JoA1Dq4vAMaH9fmx\nxofxOVjjgM+PNU7yvc/BGAd8yff4/JDaZ3wOOP7ka+NgnOR+n+NgfH6Mz8/iCzdl+Ccdm4I9J/r+\nAYwe7K4Z/v98Nzr95rK7iQQl277BWeYwj+S68gVL+n/qmLN0NWyG9pN7Ye3L6enqJFxQOI6faGa2\nRCxKZ2sDna0N9LQ30tveRLSziXhXK4lIO0TaobcTJ9aBP95NKNFJyO0m7HZRaHsost0ETZwSoGSM\nc/XaAF2E6TEF9PgKifoK6HSKiQXmEvcX4fqLsMFinIISggUlBIrLCVUspqBqMdU1KwiFQrR3tNLW\ncJLWU4dYtuZ6NpRW4CYSPP/iA0SqLyXYvJ8NLb9OnvDMD2mggt3VN3Ne46NU0UpRrAUn1ojBpTm8\nhPrwpSQCxRAsxoZKMKESfOFSAuFCAsEwwVAIxx/ECYTw9/0KBvEHQgRCIQKBMIFgEMcfAidAwBgC\nQPaeXjx9KdhzoH8oZvRgHyn4bWz6zYrZ/uu7WBc/wNY1n2fe858CoHLJhen9cxetACDefJz6uiOE\nvnMVu+bcwMYP3+tJvbnS1dlOy9mTdDSdoqflFNG2s9iOenzd9QR6WwjG2gjHOyhy2ymxnRTTQxlQ\ndo72IjZAlymkxxQS8RUSdYroCM2jxV9EPFCMGyyBYDGkAtEXLsFfUEqwsIxgYQmhojIKissoLC6j\nMBymcsy/j6MrLZ9DafkcFp9/WXqbz3FY88nfpt93trdw+Lk/kIh2sfq6v2BjKMyZEwd58cR+Vl95\nU/q46fkt0vShYM+BdLCPY3XHYdumWbD3dHWwePuXOOBfwcte9wFIBXtpxdz0MQXFpTRTimk/wZGn\nfsZGutjY9Et6I92Ewpn+AD41or29NJw8QOvpQ3TXHyXRcgKn4ySFPacpi56lwm2hyEQoGuGzLZTQ\n7iujxymhM1hNc2AFiXA5brgcU1CBU1RJoHgO4dI5FJRWUVw2h5KySsLhMOEp/5Nmpri0gkuvHzx+\nPH/xCuYvXuFRRbOTgj0nxjcUM/RBHHHrg9j0WQgs0tPFzh/9E5tooulV/4lvwAPBi0orBh3b5J9H\nQXcdkdP9U+OO7dmaF2Pu9ScP4wSCBEJh2hpO0d1WT9uxF3DP7iXccYQ5kRPMd89SYxLUpD7jWkOj\nqaDZP4+zRRdQVzgXWzwPp2QeofL5FFYuoKx6EeVVC6gIhqgYtQKR7FKw50DcKYIEBAKBUY/r67Hv\n3PAlCioWUPPQ+zCx6THG3tHWTNtXr2CTPcP2omtZu+nGQfv9geDg48MLqOo6RIcbp50iSulKj7l7\n4fj+HZx65gHKjz/Mqvje9PbSAcf02CCn/QtpKDqfurJX46taQWH1MioWnkdVTS1zQwXMHd60iOcU\n7DlQ/Z77eP7Ru1lz4fpRj7Op6Y6X3XQbAE2/C2ES+R/siXicQ3fdyuX2DPv8q1j89q+P+ZloUQ1z\nOzYTiEY5VLSGNV1PEms6mvNaX3r2YaKP/huFsWa6AxUUxlqoSDSyhDaWAAed83hm7l9Q0byD5vlX\n4a9eSbCkiorFK1m0/GKWO5Od1CbiHQV7DlTOX0rl2+8Y+8AhQzG9JoQvnt9DMdZanvn2X3FV11M8\nvfKTXPn2z4zvg+VLCJ+NsYAGjlW8nsauF3Daj+eszkO7nqbzt5/lsp4t6W0n3QU0hRfTHL6Qg/Mu\nZumVt7Biyflo9FdmGgW7h0qvft+g9zETwhfP7x77U/f9K1c33M+z8986Yqg33LaT3u5OFg3ZHq5e\nBvuSr03pfJr88ynsrstqbZ0drex+4IvMP/4g57nHaKeIzcs+xGV//knqDu5i2eqNLBoyRCQyEynY\nvfLZNlYO2RQ1Ifx5/BSl5x76AVfu/3d2Fl/N+vd/c8Rjht552Kd0/vL0a6eogo6CGhZ27h52XDwW\n5bkH/gO3u5nzb/4IVfMWseeZh1jw0Ps5VLKBl33sp8Pnv1vLtl9/l2Xb/pkraOOlwEVsrnk/F735\ndq6oqALg/MuvmdwfWmQaUrDnkbgvhJPo9bqMET3506+y6YU7OBhcyQV/82OMM7G/OlWL+gc8goUV\ndJUvY377o0R6uggX9E8S3PrD/8cVR+8EYM/3NlP16Sfx/eEOKmhnXccj7N36CBdufPWgtp+5+xNs\nOvFd9vkvoOm197DqZa/I4E8qMv3N7lv/8kzcCRN082soxrouT939Sa5+8Z/YX3g5Cz70EOHCse5z\nHG7gokfB4nICc1fiM5bTR/YMOtfyo/ezK7yeLSs+xkXR3ezd8jtWxfawuebddNkwHc/cM6jd/duf\nYNOJ7/Js2Y2suP1pVirURRTs+STuhAm4+dVjf/qeT3PV8W/zXPmNrPzYbykpq8y4zYKSSsqXrAag\n5dgL6e0Hdz7JPJqIXPAGll3/TlxrcB7+RwDKLrmR/cXr2dDya9pbm2htPMO+bY9Q+qv30Eg5F77n\nThy/fgAVAQ3F5JWEU0jI5s8Y+3O/+R5XHb+T58pezdq//VHW1nYpKp1DQUlyzL337L709sZnf8oy\n62PlNW+hvGo+BwIrWBnfT8w6LL/0ag7Go/DIn+j82hVUuY1cYBKcpprG193Nqiz8D0dkplCPPY8k\nAvkT7KcO7+GCLZ/hpcBFXPJXP8jqgl3F5XMoKimnnkr8LYcA6Opo5aLTD/BC0ab08q6N1VcAUOfU\nEC4s5uKrX089lSy0ZwmaBM8VXUPor59g1bpXZq02kZlAwZ5HrL+AgjwI9li0l4773oVrDGW33ksw\nnJ318nYWbCBq/en1YepDSyjrOgrAi7+7mzK6CF//8fTxodoNAHT5+9dwb3vjD9hS+Xq6Pn6Ml/39\ng1TOrUFEBtNQTD5xgviJe10F2+79B66I72Pb+i+zbunQSZmTt/rvfk1HaxMVqd5/V9ESLm7+Fc/c\n9y8UHn+ck2Y+Fwz48nPFxpt4fueVlLzmH9Pbzr/8GtDURZFRZSXYjTE3Al8j+VCR71prv5CNdmcd\nJ0jQJLCu69la5S9tfYSNJ+/h2Yqb2HDze7Patj8QpKJ6Qfq9W14LzbBp/5cAeGbe21g04M9dWj5n\n0JKwIjI+GaeHMcYBvgm8FrgIeJsx5qJM252VnBAA0Zg3M2Os68Lv/pFGU8Hq996Z8/MFKhcPer/o\nhg/l/Jwis0E2uoUbgIPW2sPW2ijwY+ANWWh39vEnV4OM9Xozzr7jkftYFd/L0Ys/PCXPJl2+4Wbq\nTP9D2BatuDjn5xSZDbIR7DXAiQHvT6a2yQQZJ7mOSTw69T32eCxK5TNf4LivhrVvmJqec+XcGuZ+\nahfdNsQz53987A+IyLhM2ZenxpjbgNsAlizRg7FGYvzeBftzv/w6G90TbL/iGyyZwoWyAsEQgTvq\nyf3jfUWRyCwWAAALFElEQVRmj2z02OuAgYOli1LbBrHW3mWtXWetXVddXZ2F085A/uQY++H//jBP\n3v+VKTvtrsd+xpoXPs+LwUtYc8P/nbLzikhuZCPYtwLnG2OWGWOCwFuBX2Wh3VnHlxqK2dD9BFfv\n+eyUnLO7s43aJ/6WOmcRiz7wgGezcUQkezL+V2ytjQMfAn4H7AXut9a+mGm7s5EvEMpJu1v/505O\n3bGSRHz4HPk9334PpXQRec0XKZszLyfnF5GplZUxdmvtb4DfZKOt2cznz02wr3/+dgA6utrTi3g9\n81//j5rDP2WdPcWu8Hou3fianJxbRKae7jzNIz5/br+07O3poqSsklNH97Hp0NfS28vf/O85Pa+I\nTC0Fex7x5Xg2St297+Xky97N+X/6KBECdJkiDldcxfqVl+f0vCIytRTseSRcXDHi9mcf/A7GF2D9\nTe+acJtuIpH+IuWyni3w5BbaTBEtt/6BRSsuZs6onxaR6UhTIPJISeX8Ebdv2PYJ1j/7kUm1eXDn\nk8O2nQit1F2eIjOYeux5pKCoLOttNu3qX0TrLHNo9VdR/EaNqYvMZAr2PBIqKhr7oAnYs/m36QdD\nAzQGa1j96T9l9Rwikn80FJNHgsHhD7SI9HSlXz/7H7fQ0nh23O2ZRz8HwD7/KgB6g7lf2EtEvKdg\nzyMj3fXZ3lyffr2h42EOPvXAuNrqaGtmefQAz1b+Gcv+/gm2VL2ZJbf+Z9ZqFZH8paGYPNfRXM/c\nAe/HO9d976M/ZIOJUbrxLwmGwmz80PdzU6CI5B312POYdV162gYPvSR62sf8XKSni3m7v8VxXw0X\nrH9VrsoTkTylYM9jsViM3vamQdvcyNjB/sKdf8lS9yT1az6iRb1EZiH9q89j8Xgv8Y6GQdtsb8eo\nn6k7vId17Y/wzLy3se7PPpDL8kQkTynY89gLf/wFbtfgYDdjBPup3Y8DMP+69+WqLBHJc/ryNM/s\nCq/n0shWADZs+dth+01P86ifjzcdBWDBMj1PXGS2Uo89z0QWrB91/8aOh3n+D/efc7/TfoJ6KgmF\nC7NdmohMEwr2PLP+Hf/CtgVvG/WYyKHh67/0Kequozkw8pozIjI7KNjzjPE5OEs2jn6Qc+657BWx\nM3QWLMxyVSIynSjY85AvNHxpgYGMMSNuj8eiVLuNxEoW5aIsEZkmFOx5yB8YPditGfk/W+PpYwRM\nAqdiaS7KEpFpQsGeh5xQ8ovPU2YuZ9+1mUbGt3hX86lDABRU1+aqNBGZBhTsecifGoqJmSDzai+i\nIbh40H5jnBE/F2k9A0DRnJrcFigieU3BnsfiJvklac0HHuDF67/HQf+K5I5zDMXEOhoBKKmcNyX1\niUh+UrDnIcdJ9sg7wsnZLaVz5rH65bdQ/cEHAbC4I37O7UoGe2nl3BH3i8jsoDtP89Cy1Zt4/sAd\nnP+KdwzaXlxenXyRiI/4OdPdTKctoFg3J4nMagr2fGQMa9700WGbHb8f1xpwYyN+zB9ppt1XQnGu\n6xORvKahmGkmjjNij926LpWdB+hysv9AbBGZXtRjn2LbL7+DgrkruHCSn4/jYEbosW978C7Wu0c5\n41ZlVqCITHsK9im29o3Dh1gmIm4ccIf32N3jWwAos2M/iENEZjYNxUwzCfwj9thdJwTAsdf9aKpL\nEpE8ox77NJMciunvsXd2tLL3rnczt+sIJ818VukZpyKznoJ9mkngxyR60+/3PHwvGzoeTb4OXuJV\nWSKSRzQUM80U207Wtz/M0/d+BgD/wd+n97nnWGpARGYXBfs0U2x6ALjyyDfo6eqktntXel/Pyjd5\nVZaI5BENxUxjO//3m2yinedKX0lw7dtZf90tXpckInlAPfZpLHT8cQB8q9/IJQp1EUlRsE9TMeuw\npvtpAAoq9Sg8EemnYJ+mAiaRfl0+b4mHlYhIvsko2I0xXzLGvGSM2WWM+YUxZnyP+pGs6bYh5tac\n53UZIpJHMu2xPwxcbK29FNgPfCrzkmQ8tlS/BYDT/oX4HE1zFJF+Gc2Ksdb+fsDbZwB9g5dj+157\nP4lYBE7ugQboDOppSSIyWDbH2N8D/PZcO40xtxljthljtjU0NGTxtLPLBRtfw0VXvwHrJsfYI4Xz\nPa5IRPLNmD12Y8wjwEjp8Rlr7f+kjvkMEAd+eK52rLV3AXcBrFu3zk6qWukX7Ur+Xl7raRkikn/G\nDHZr7airShlj3gW8DniltVaBPUUufvPfs+W+Vi778094XYqI5JmMxtiNMTcCnwRebq3tzk5JMh7F\npZVs/OCdXpchInko0zH2bwAlwMPGmB3GmG9loSYREclAprNiVmSrEBERyQ7deSoiMsMo2EVEZhgF\nu4jIDKNgFxGZYRTsIiIzjIJdRGSGMV7cLGqMaQCOTfLjVUBjFsvJNtWXuXyvUfVlRvVN3lJrbfVY\nB3kS7Jkwxmyz1q7zuo5zUX2Zy/caVV9mVF/uaShGRGSGUbCLiMww0zHY7/K6gDGovszle42qLzOq\nL8em3Ri7iIiMbjr22EVEZBTTKtiNMTcaY/YZYw4aY273up6hjDFHjTG7U0sYb8uDeu42xtQbY14Y\nsK3SGPOwMeZA6veKPKvvs8aYutQ13GGMucnD+hYbYx4zxuwxxrxojPlIanteXMNR6sunaxg2xjxr\njNmZqvGO1PZlxpgtqX/LPzHGBPOsvnuMMUcGXMPLvahv0qy10+IX4ACHgOVAENgJXOR1XUNqPApU\neV3HgHquBdYCLwzY9kXg9tTr24F/y7P6Pgt8wutrl6plAbA29boE2A9clC/XcJT68ukaGqA49ToA\nbAE2AfcDb01t/xbwV3lW3z3ALV5fv8n+mk499g3AQWvtYWttFPgx8AaPa8pr1to/As1DNr8BuDf1\n+l7gjVNa1ADnqC9vWGtPW2u3p153AHuBGvLkGo5SX96wSZ2pt4HULwu8AvhZaruX1/Bc9U1r0ynY\na4ATA96fJM/+EpP8C/F7Y8xzxpjbvC7mHOZZa0+nXp8B5nlZzDl8yBizKzVU49lQ0UDGmFpgDcke\nXd5dwyH1QR5dQ2OMY4zZAdQDD5P8ybvVWhtPHeLpv+Wh9Vlr+67hv6au4VeMMSGv6puM6RTs08HV\n1tq1wGuBvzHGXOt1QaOxyZ8/8613cidwHnA5cBr4D2/LAWNMMfAA8FFrbfvAfflwDUeoL6+uobU2\nYa29HFhE8ifvVV7WM9TQ+owxFwOfIlnneqAS+AcPS5yw6RTsdcDiAe8XpbblDWttXer3euAXJP8S\n55uzxpgFAKnf6z2uZxBr7dnUPzQX+A4eX0NjTIBkaP7QWvvz1Oa8uYYj1Zdv17CPtbYVeAy4Aig3\nxvQ9mjMv/i0PqO/G1DCXtdb2At8nT67heE2nYN8KnJ/6Nj0IvBX4lcc1pRljiowxJX2vgVcDL4z+\nKU/8Cnhn6vU7gf/xsJZh+gIz5U14eA2NMQb4HrDXWvvlAbvy4hqeq748u4bVxpjy1OsC4AaS3wU8\nBtySOszLazhSfS8N+B+3ITn+n4//ls9pWt2glJq29VWSM2Tuttb+q8clpRljlpPspUPyIeH3eV2f\nMeZHwHUkV6s7C/wT8EuSMxKWkFxh8/9Yaz35AvMc9V1HcgjBkpxl9IEB49lTXd/VwJ+A3YCb2vxp\nkuPYnl/DUep7G/lzDS8l+eWoQ7Ijeb+19nOpfy8/JjnM8Txwa6p3nC/1PQpUk5w1swP44IAvWfPe\ntAp2EREZ23QaihERkXFQsIuIzDAKdhGRGUbBLiIywyjYRURmGAW7iMgMo2AXEZlhFOwiIjPM/wcR\nfLvwDLLErgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fbf3a65d790>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# check that the real dielectric constants match\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "eps_real_calc = kkr(energies[1] - energies[0], eps_imag)\n",
    "\n",
    "plt.plot(energies, eps_real[:, 0, 0])\n",
    "plt.plot(energies, eps_real_calc[:, 0, 0])\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}