{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dispersion Coefficients Near Surfaces"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this tutorial we reproduce the results depicted in Figure 5 from J. Block and S. Scheel \"van der Waals interaction potential between Rydberg atoms near surfaces\" [Phys. Rev. A 96, 062509 (2017)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.062509). We calculate the van der Waals $C_6$-coefficient between two Rubidium Rydberg atoms that are equidistantly placed in front of a perfect mirror (i.e. in horizontal alignment in front of a perfectly conducting plate). One finds that the relevant length scale is interatomic distance devided by distance from surface and that for decreasing surface distance the $C_6$ coefficient is significantly reduced."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As described in the [introduction](introduction.ipynb), we start our code with some preparations and load the necessary modules. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "# Arrays\n",
    "import numpy as np\n",
    "\n",
    "# Plotting\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Operating system interfaces\n",
    "import os\n",
    "\n",
    "# pairinteraction :-)\n",
    "from pairinteraction import pireal as pi\n",
    "\n",
    "# Create cache for matrix elements\n",
    "if not os.path.exists(\"./cache\"):\n",
    "    os.makedirs(\"./cache\")\n",
    "cache = pi.MatrixElementCache(\"./cache\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plate lies in the $xy$-plane with the surface at $z = 0$. The atoms lie in the $xz$-plane with $z>0$.\n",
    "\n",
    "We can set the angle between the interatomic axis and the z-axis `theta` and the center of mass distance from the surface `distance_surface`. `distance_atom` defines the interatomic distances for which the pair potential is plotted. The units of the respective quantities are given as comments.\n",
    "\n",
    "Be careful: `theta = np.pi/2` corresponds to horizontal alignment of the two atoms with respect to the surface. For different angles, large interatomic distances `distance_atom` might lead to one of the atoms being placed inside the plate. Make sure that `distance_surface` is larger than `distance_atom*np.cos(theta)/2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = np.pi / 2  # rad\n",
    "distance_atoms = 10  # µm\n",
    "distance_surface = np.linspace(\n",
    "    distance_atoms * np.abs(np.cos(theta)) / 2, 2 * distance_atoms, 30\n",
    ")  # µm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we define the state that we are interested in using pairinteraction's `StateOne` class . As shown in Figures 4 and 5 of [Phys. Rev. A 96, 062509 (2017)](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.96.062509) we expect changes of about 50% for the $C_6$ coefficient of the $|69s_{1/2},m_j=1/2;72s_{1/2},m_j=1/2\\rangle$ pair state of Rubidium, so this provides a good example. \n",
    "\n",
    "We set up the one-atom system using restrictions of energy, main quantum number n and angular momentum l. This is done by means of the `restrict...` functions in `SystemOne`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "state_one1 = pi.StateOne(\"Rb\", 69, 0, 0.5, 0.5)\n",
    "state_one2 = pi.StateOne(\"Rb\", 72, 0, 0.5, 0.5)\n",
    "\n",
    "# Set up one-atom system\n",
    "system_one = pi.SystemOne(state_one1.getSpecies(), cache)\n",
    "system_one.restrictEnergy(\n",
    "    min(state_one1.getEnergy(), state_one2.getEnergy()) - 30,\n",
    "    max(state_one1.getEnergy(), state_one2.getEnergy()) + 30,\n",
    ")\n",
    "system_one.restrictN(\n",
    "    min(state_one1.getN(), state_one2.getN()) - 3,\n",
    "    max(state_one1.getN(), state_one2.getN()) + 3,\n",
    ")\n",
    "system_one.restrictL(\n",
    "    min(state_one1.getL(), state_one2.getL()) - 1,\n",
    "    max(state_one1.getL(), state_one2.getL()) + 1,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The pair state `state_two` is created from the one atom states `state_one1` and `state_one2` using the `StateTwo` class.\n",
    "\n",
    "From the previously set up `system_one` we define `system_two` using `SystemTwo` class. This class also contains methods `set..` to set angle, distance, surface distance and to `enableGreenTensor` in order implement a surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up pair state\n",
    "state_two = pi.StateTwo(state_one1, state_one2)\n",
    "\n",
    "# Set up two-atom system\n",
    "system_two = pi.SystemTwo(system_one, system_one, cache)\n",
    "system_two.restrictEnergy(state_two.getEnergy() - 3, state_two.getEnergy() + 3)\n",
    "\n",
    "system_two.setAngle(theta)\n",
    "system_two.setDistance(distance_atoms)\n",
    "system_two.setSurfaceDistance(distance_surface[0])\n",
    "system_two.enableGreenTensor(True)\n",
    "system_two.buildInteraction()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We calculate the $C_6$ coefficients. The `energyshift` is given by the difference between the interaction energy at given `surface_distance` and the unperturbed energy of the two atom state `state_two.getEnergy()`. The $C_6$ coefficient is then given by the product of `energyshift` and `distance_atoms**6`.\n",
    "\n",
    "`idx` is the index of the two atom state. The command `getOverlap(state_two, 0, -theta, 0)` rotates the quantisation axis of `state_two` by `theta` around the y-axis. The rotation is given by the Euler angles `(0, -theta, 0)` in zyz convention. The negative sign of theta is needed because the Euler angles used by pairinteraction represent a rotation of the coordinate system. Thus, the quantisation axis has to be rotated by the inverse angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate C6 coefficients\n",
    "C6 = []\n",
    "for d in distance_surface:\n",
    "    system_two.setSurfaceDistance(d)\n",
    "    system_two.diagonalize()\n",
    "    idx = np.argmax(system_two.getOverlap(state_two, 0, -theta, 0))\n",
    "    energyshift = system_two.getHamiltonian().diagonal()[idx] - state_two.getEnergy()\n",
    "    C6.append(energyshift * distance_atoms**6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAEKCAYAAAAxXHOuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xl4FeX5//H3nYQkENawb2ENqKggxH1FrFutW9Val7q0RetSta2tdrF2+bXa1taq31qptWqrVqsVqSu4r6AgyI7sSwQCBEIgZL9/f8xEDjHASTg5k+Xzuq5zZeaZZ2buM2dy7jPb85i7IyIikigpUQcgIiItixKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgkVFrUATSGbt26+cCBA6MOQ0SkWZkxY8ZGd+++r8tpkYll4MCBTJ8+PeowRESaFTNbmYjl6FSYiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgklBKLiIgkVIt8jkVEpDWqrnZKK6vYUV5FaWU1pRVVlFVUU1ZZRVlldfCqiBmurJkeDCeKEouISJJUVlWzvayK7eWVbC+rZHt5FSVllZSUV1FSUUVpeRUl5ZUxw1W7DO+oCF6ln7+C5LEjTCDlVdVRv0VAiUVEZI/cnW1llWwrq6S4tOZVEZTVjJeFZeH454kjTCIl5VVsK6ukvDL+L34zaNsmlXbpqbRNT6Vtm1TapqfRtk0K2VnpZKalktkmhbbpqWSkpZLZJhjPbBPUzWyTEpYHfzPSUsiIHU5LDceD4fS0FNLuTMw2U2IRkVahvLKaLSXlFJaUU7i9nKKSCop27HxtCf9u3bFr+dYdFVT7npdtBu3T0+iQmUb7zDSyMtJon5FG9w4ZZKUH4+0yUmmfnka7jDSy0lPJykgjKyOVdulptEsPEkhmm53jGWkpmFlyNk6CRZJYzKwz8CBwIODAlcApwLeBDWG1H7v7i2H9W4FvAlXAd939laQHLSJNSlllFRu3lbOxuIwNxWVs3FZGYUk5m7eXU7i9gs0l5WzaHoxv3l5OcVnlbpeVlmJ0atuGTm3b0LFtG7q0S2dg16yYsjQ6ZLYJEkfGzuGa8az0NFJSmmcSaAxRHbH8GXjZ3c8zs3SgHUFi+ZO7/yG2opkdAFwIjAD6AK+a2TB3T9yVJhFpMraVVbKuaAdri0op2FrGhm1B4vj8FY4X7aioc/7MNil0zcqgS1ZNgmhHl3bpdM1Kp0tWOtlZ6XRuF0yrSRzt0lOb7dFBU5T0xGJmnYDjgMsB3L0cKN/Dh3oW8G93LwOWm9kS4DDgg8aPVkQSaUd5FZ8V7WDtllLWhsljbdEOPttSyrqiUj4r2kFx6RePLNq2SaVHxwy6t88gt0d7jhrSle7tM+jWISjr3iGDru3T6ZqVQdv01AjemcSK4ohlEMHprn+Y2UhgBnBDOO06M/sGMB34vrtvBvoCU2PmXxOWiUgTU13trC8uZdWmElYWlrC6sISVm0pYVRi8CreXf2Gebu3T6dUpk5yu7ThicDa9O7eld6dMenXMpGfHzOA6RYYuBzcnUXxaacBo4Hp3n2ZmfwZuAe4DfkVwzeVXwF0E117iYmbjgfEAOTk5iY5ZRELuzrqtpSxev40lBds+TxorN21n9eYdu9z5lJpi9OmcyYDsLE4Z0Yt+XdrSp3MmvTsFyaNnx0wy2+gIo6WJIrGsAda4+7Rw/GngFndfX1PBzP4GPB+O5gP9Y+bvF5btwt0nABMA8vLy9nIPh4jsTXW1k79lB4sLilm8fhuLC4LX0oJtbIu5EJ6VnkpO1yxye3Rg3P49yclux4Cu7cjJbkefzm1pk6oGPlqbpCcWd19nZqvNbLi7LwLGAfPNrLe7rw2rnQPMDYcnAY+b2R8JLt7nAh8mO26RlqyopII5+UXM/ayIReuKWVxQzJKCbZRW7Dz66N4huL7x1dF9GdqzA7k92jO0R3u6ZqXrwrfsIqoTl9cDj4V3hC0DrgDuMbNRBKfCVgBXAbj7PDN7CpgPVALX6o4wkYbbUlLOnPyiIJGEf1cX7vh8eu9OmeT27MDFh3f9PHkM7dGezu3SI4xamhNzb3lnjfLy8lx93otAcWkFs1Zv2W0S6Z/dloP6duLAvp2Cv3060SVLCaS1MrMZ7p63r8vRrRYiLUjRjgo+Wl7ItOWbmLa8kLn5RZ8/NZ6T3Y6D+3bmosMGhMmko45CpFEosYg0Y4Xby/kwTCLTlhWyYN1W3CE9NYVR/Ttz7dihHDowm4P7dVISkaRRYhFpRraVVfLu4g28t2QT05Zv4tP124DgafPROV24cdwwDh+czaj+nXUbr0RGiUWkicvfsoPXFqzn1QUFTF26ifKqarLSUxkzMJuzRvXliMHZHNS3M+lpuq1XmgYlFpEmprramftZEa/OD5LJ/LVbARjcLYvLjx7IuP16MHpAFz0fIk2WEotIE1BaUcX7SzcyZX4Bry9cz/qtZaQY5A3I5sen78e4/XsypHv7qMMUiYsSi0hEqqud95du4j8zVjN53np2VFSRlZ7K8cO7c9L+PTlheA+ydeuvNENKLCJJtnLTdp6esYZnZqzhs6JSOmamcc7ovpwyohdHDM4mI00X3aV5U2IRSYLtZZW8MGctT89Yw4fLC0kxODa3O7eevj9fOqCn7uCSFkWJRaSRuDvTlhfy9Iw1vDhnLSXlVQzulsXNpwznq6P70atTZtQhijQKJRaRBNtSUs5j01bx5EerWVVYQvuMNM4c2Yfz8/oxOqeLGmyUFk+JRSRB1m8t5cF3lvH4tFVsL6/iqCFduelLuZw6ord6NZRWRYlFZB+t2LidB95eyjMz8qly5ysH9+bqE4awX6+OUYcmEgklFpEGmvdZEfe/uZQX56wlLTWFCw7tx1XHDaF/druoQxOJlBKLSD19uLyQv7y5hDcXbaB9RhrjjxvClccMpEcHXYwXASUWkbi4O28sKuAvbyxl+srNdM1K5+ZThnPJEQPo1LZN1OGJNClKLCJ7MTe/iNuem8vHq7bQt3NbfnHmCC7I668L8iK7EUliMbPOwIPAgQRdEV8JnAt8BSgHlgJXuPsWMxsILAAWhbNPdferkx2ztD6bt5fzh8mLePzDVXTNSufOrx7EuaP7qfFHkb2I6ojlz8DL7n5e2O99O2AKcKu7V5rZncCtwI/C+kvdfVREsUorU1Xt/PujVfz+lUUUl1ZyxVGDuPFLuXTM1CkvkXgkPbGYWSfgOOByAHcvJzhKmRxTbSpwXrJjE5mxcjM/nzSXuflbOXxQNr8860CG9+oQdVgizUoURyyDgA3AP8xsJDADuMHdt8fUuRJ4MnYeM5sJbAV+6u7vJC1aaRU2FJdx58sLeXrGGnp1zOTerx/CGQf31lPyIg0QRWJJA0YD17v7NDP7M3AL8DMAM/sJUAk8FtZfC+S4+yYzGwNMNLMR7r41dqFmNh4YD5CTk5OcdyLNXmVVNY9+sJI/TfmU0soqvnPCEK4bO5SsDN3XItJQUfz3rAHWuPu0cPxpgsSCmV0OnAGMc3cHcPcyoCwcnmFmS4FhwPTYhbr7BGACQF5enjf+25DmbuqyTfz8uXksWl/MccO6c/tXDmCwOtMS2WdJTyzuvs7MVpvZcHdfBIwD5pvZqcAPgePdvaSmvpl1BwrdvcrMBgO5wLJkxy0tR0VVNb9/ZRET3l5Gvy5tmXDpGL50QE+d9hJJkKiO968HHgvvCFsGXAF8BGQAU8J/8Jrbio8DfmlmFUA1cLW7F0YTtjR3azaXcP0TM5m5aguXHJHDT798gPpCEUmwSBKLu88C8moVD91N3WeAZxo9KGnxJs9bx81Pz6aq2rnvokM44+A+UYck0iLpCqW0eOWV1dzx0kIeem85B/btyP9dNJoBXbOiDkukxVJikRZtdWEJ1z3+MZ+sKeLyowZy6+n7qU95kUamxCIt1stz13Lz07MB+Oslozn1wN4RRyTSOiixSItTVlnFb15YwCMfrGRkv07cd9Fo9ZEikkRKLNKirNy0nesen8mc/CKuPHoQt5y2H+lpajRSJJmUWKTFmDJ/Pd97chZmMOHSMZw8olfUIYm0Skos0iL8Z/pqfvTMbA7q24n/u3g0/bro1JdIVJRYpNl76N3l/PL5+Ryb240HLh1Du3Tt1iJR0n+gNFvuzp9fW8zdry7m1BG9+PPXR+lWYpEmQIlFmqXqaufXLyzgofeWc96Yftxx7kGkqWdHkSZBiUWancqqam797xz+M2MNVxw9kJ99+QBSUtSApEhTocQizUpZZRU3PDGLl+et48aTcrlhXK5aJRZpYpRYpNkoKa/kqn/O4J3FG7ntjAO48phBUYckInVQYpFmoaikgisf+YiZqzbz+/MO5vy8/lGHJCK7ocQiTd6G4jIu/fs0lm7Yxl8uVptfIk2dEos0aWs2l3Dp3z9kXVEpf7/sUI4b1j3qkERkL5RYpMlas7mEC/76AdvKKvnXtw5jzIDsqEMSkThEcuO/mXU2s6fNbKGZLTCzI80s28ymmNni8G+XsK6Z2T1mtsTMZpvZ6ChiluTaWlrBlQ9/RHFZJU+MP0JJRaQZieqJsj8DL7v7fsBIYAFwC/Cau+cCr4XjAKcBueFrPHB/8sOVZKqoqubaxz5m2YbtPHDJGEb06RR1SCJSD0lPLGbWCTgO+DuAu5e7+xbgLOCRsNojwNnh8FnAox6YCnQ2M129baHcnduem8s7izfym3MO4qih3aIOSUTqKYojlkHABuAfZjbTzB40syygp7uvDeusA3qGw32B1THzrwnLpAWa8PYynvhwNdeOHcIFh+qWYpHmKIrEkgaMBu5390OA7ew87QWAuzvg9VmomY03s+lmNn3Dhg0JC1aS58U5a/ntSws54+DefP9Lw6MOR0QaKIrEsgZY4+7TwvGnCRLN+ppTXOHfgnB6PhD707VfWLYLd5/g7nnunte9u25JbW5mrtrMTU/OYnROZ/5w/ki1/SXSjCU9sbj7OmC1mdX8JB0HzAcmAZeFZZcBz4XDk4BvhHeHHQEUxZwykxZgdWEJ3350Oj07ZvK3b+SR2UZN34s0Z1E9x3I98JiZpQPLgCsIktxTZvZNYCVwQVj3ReB0YAlQEtaVFqJoR3BbcXllNf8efyhd22dEHZKI7KMGJZbwYnupu1c1ZH53nwXk1TFpXB11Hbi2IeuRpq2iqpprHpvBik3befTKwxnao33UIYlIAsR1KszMUszsIjN7wcwKgIXAWjObb2a/N7OhjRumtDTuzk+fnct7Szbx23MP5sghXaMOSUQSJN5rLG8AQ4BbgV7u3t/dewDHAFOBO83skkaKUVqgv761jCenr+b6E4dy3ph+UYcjIgkU76mwk9y9onahuxcCzwDPmFmbhEYmLdYLs9dy58sLOXNkH773pWFRhyMiCRbXEUtdSaUhdUQ+XrWZm56aRd6ALvzuvIPV+6NIC7TXIxYzu5Tg2ZNrgErgbXdXe11Sb1tKyrnmXx/Tu1MmE3RbsUiLFc+psEOBke5+PoCZ/blxQ5KWyN35ycS5bNxWxrPXHE12VnrUIYlII4knsWwF+pnZt4HNQFbjhiQt0cRZ+bwwey03nzKcg/qptWKRliyeayw/AyYC2UAGwcONInFbs7mE2ybO49CBXbj6+CFRhyMijWyvRyzhA4oTa8Z195fUR1W1872nPsGBP14wilS1ASbS4tXryXsz+xtwhplVAp8Bs4HZ7n5vYwQnzd+Et5fx4fJC7jp/JP2z20UdjogkQX2bdDkO6OfuVWbWl6D3x4MTH5a0BHPzi/jjlEV8+aDenDtaXeiItBb1TSzTgK5AgbvnEzRf/2LCo5Jmr7SiihufnEV2Vjr/75wD9byKSCtS32bzHwDeMrMfmNmxYTfDIl9wx0sLWVKwjT+cP5LO7XRrsUhrUt/E8i/gUYIjnWuA981sacKjkmbtzUUFPPz+Cq48ehDH5qrTNZHWpr6nwta4+29jC8xMHWjI5wq3l3Pz07MZ1rM9PzxV3QuLtEb1PWKZZWY3xBa4e1kC45FmzN358X/nUFRSwd1fO0RNtoi0UvVNLD2Bq83sMzN73sz+n5md3xiBSfPznxlreHneOn5wyjAO6NMx6nBEJCL1OhXm7hfA56e/RgAHAYcB/6nPcsxsBVAMVAGV7p5nZk8CNedOOgNb3H2UmQ0EFgCLwmlT3f3q+qxPGt+qTSX8YtI8jhiczbeOGRx1OCISofo+IDmH8KHI8PUa0NBemsa6+8aaEXf/Wsx67gKKYuoudfdRDVyPNLLKqmpuemoWKSnGXReMIkVP14u0avU9FXY88DdgB3AhMBc4PZEBWfDAwwXAE4lcrjSe+99cyoyVm/n12QfSt3PbqMMRkYjVK7G4e6G7v+nu97j7ZQRN6i9pwHodmGxmM8xsfK1pxwLr3X1xTNkgM5tpZm+Z2bENWJ80krn5Rdz92mLOHNmHs0bp6XoRqf+psGHu/mnNuLsvNrOGNOlyjLvnm1kPYIqZLXT3t8NpX2fXo5W1QI67bzKzMcBEMxvh7ltrxTYeGA+Qk5PTgJCkvtydn0+aR+e2bfjVWQdGHY6INBH1fvLezFaZ2Qdm9oCZPQLMNbN6tS4YNgeDuxcAzxLcAICZpQHnAk/G1C1z903h8AxgKfCFjtLdfYK757l7XvfueigvGZ6dmc+MlZv50Wn70amdGr0WkUB9T4WNdfcc4GvA8wSnwdoSPN+yMJ5lmFmWmXWoGQZOJrhWA3ASsNDd18TU725mqeHwYCAXWFafuCXxiksr+M2LCxnVvzPnjW7o/Rsi0hLV98l7ANx9FbAK+F9NmZm1j3P2nsCzYaOEacDj7v5yOO1CvnjR/jjgl2ZWAVQDV7t7YUPilsS557XFbNpext8vy9NdYCKyiwYllrq4+7Y46y0jaG6/rmmX11H2DPDMPgUnCbWkoJh/vLeCr+X1Z2T/zlGHIyJNTH2vsUgr5+7cPmk+7dJTufkUtQUmIl+kxCL18sq8dby7ZCPfP3k4Xdur/VER+aJ6JRYze83MTq9VNiGxIUlTtaO8il89v4D9enXg4sN1S7eI1K2+RyyDgB+Z2c9jyvISGI80Yfe/tZT8LTv4xZkjSEvVwa6I1K2+3w5bgHFATzP7n3qQbD1WbSrhr28t5cyRfTh8cNeowxGRJqy+icXcvdLdryG4U+tdoEfiw5Km5lcvzCctxfjx6ftHHYqINHH1TSx/rRlw94eBy4DJiQxImp43FxUwZf56rj8xl16dMqMOR0SauLieYzGzewkajsTM7qk1Oa7nV6R5Kq+s5pf/m8+gbllceczAqMMRkWYg3gckp8cM/wL4+e4qSsvy0HvLWbZxOw9fcSgZaepqWET2Lq7E4u6P1Ayb2Y2x49JyrSsq5d7XFnPS/j05YbgupYlIfBpyz6gnPAppkn770gIqqp3bzjgg6lBEpBnRwwhSp2nLNvHcrM+4+rjB5HStV68IItLKxXvxvpidRyrtzKymky0D3N07NkZwEo3Kqmp+PmkefTu35TsnDI06HBFpZuK9xtKhsQORpuPxD1excF0x9188mrbpumAvIvWjU2Gyi+LSCv405VOOHtqVUw/sFXU4ItIMxZVYzOwsM7s2ZnyamS0LX+c3XniSbA++s5zNJRXccur+hJ2xiYjUS7xHLD8EJsWMZwCHAicAVyc4JolI4fZyHnxnGacf1IuD+qkZOBFpmHgTS7q7r44Zf9fdN4VdFGfVd6VmtsLM5pjZLDObHpbdbmb5Ydms2Ob5zexWM1tiZovM7JT6rk/ic/+bS9hRUcX3vjQs6lBEpBmL98n7LrEj7n5dzGj3Bq57rLtvrFX2J3f/Q2yBmR0AXAiMAPoAr5rZMHevauB6pQ5ri3bwyAcrOXd0P4b20L0aItJw8R6xTDOzb9cuNLOrgA8TG9IXnAX8293L3H05sAQ4rJHX2erc+/oS3J0bxuVGHYqINHPxHrHcBEw0s4uAj8OyMQTXWs5uwHodmGxmDjzg7jW9UF5nZt8gaJvs++6+GegLTI2Zd01YJgmyYuN2nvpoNRcfnkP/bD0MKSL7Jt7nWAqAo8zsRIJTUgAvuPvrDVzvMe6eb2Y9gClmthC4H/gVQdL5FXAXcGW8CzSz8cB4gJwcdZtbH3969VPapKZw7Yl6GFJE9l28T96bB14H6kwmNXXiWZ6754d/C8zsWeAwd387Zll/A54PR/OB/jGz9wvLai9zAjABIC8vT+2ZxWnB2q1M+uQzvnP8EHp0UF8rIrLv4r3G8oaZXW9muxwKmFm6mZ1oZo8QdPq1V2aWZWYdaoaBk4G5ZtY7pto5wNxweBJwoZllmNkgIJfGv67Tatw1+VPaZ6Rx1XFDog5FRFqIeK+xnEpwWuqJ8Mt9C9CWIDFNBu5295lxLqsn8Gz48F0a8Li7v2xm/zSzUQSnwlYAVwG4+zwzewqYD1QC1+qOsMT4eNVmXl2wnptPGU6ndm2iDkdEWgiL8+zVzhnM2gDdgB3uvqVRotpHeXl5Pn369L1XbOUu+ttUPl1fzFs3jyUrI97fGCLSUpnZDHfP29fl1LutMHevcPe1TTWpSHzeXbyR95du4tqxQ5VURCSh1AhlK+Tu/P6VhfTplMlFh+sOOhFJLCWWVmjy/PV8sqaIG08apn7sRSTh4r3dON6ftVvcfeveq0lUqqqduyYvYnD3LM4dredMRSTx4j25/gjB3Vp7akfdgYeBR/cxJmlEkz7J59P12/i/i0aTlqoDVhFJvHifvB/b2IFI4yuvrOZPUxYzok9HTlMnXiLSSPSTtRV5cvpqVhWW8INThpOSok68RKRxNPg+UzO7lKBByGsIHlx8293vT1Rgklg7yqu497XFHDqwCycMa2hPByIie7cvDzAcCox09/MBzOzPiQlJGsOjH6ygoLiM+y4arS6HRaRR7Uti2Qr0C/tp2UwDepKU5NhaWsH9by3l+GHdOWxQdtThiEgLty/XWH4GTASygXTg+oREJAn30LvL2VJSwQ9OHh51KCLSCsT7HMtQoKe7v1dT5u5uZhuBOe6+tLEClH1TtKOCv7+7nJMP6MlB/TpFHY6ItALxHrHcTXDqq7aicJo0UX9/dznFpZXceNKwqEMRkVYi3sTS093n1C4MywYmNCJJmKKSCv7x7nJOHdGLA/p0jDocEWkl4k0snfcwrW0iApHE+/u7yyguq+SGk3KjDkVEWpF4E8v08O6vXZjZt4AZiQ1JEmFLSTkPvbeC0w7sxf69dbQiIskT7+3GNxL0+ngxOxNJHsHdYOc0RmCybx58ZznbdLQiIhGIt62w9cBRZjYWODAsfsHdX2/ISs1sBVAMVAGV7p5nZr8HvgKUA0uBK9x9i5kNBBYAi8LZp7r71Q1Zb2uxeXs5/3hvOV8+qDf79dLRiogkV32bzV8avmqX16hPs/lj3X1jzPgU4FZ3rzSzO4FbgR/VrNfdR8W53Fbvb+8so6Siiu+O09GKiCRffZrN35t9ajbf3SfHjE4FzmvIclq7wu3lPPL+Ck4/qDfDe3WIOhwRaYWiajbfgclm5sAD7j6h1vQrgSdjxgeZ2UyCZ2l+6u7vJDieFqPmaOVGHa2ISET2pa2wfXGMu+ebWQ9gipktdPe3AczsJwStJT8W1l0L5Lj7JjMbA0w0sxG1T7mZ2XhgPEBOTuvsx33TtjIeeX8FZxzch9yeOloRkWhE0h+Lu+eHfwuAZ4HDAMzscuAM4GJ397BOmbtvCodnEFzj+cJj5O4+wd3z3D2ve/fW2Sz8hHeWsaOiihvGDY06FBFpxZKeWMwsy8w61AwDJwNzzexU4IfAme5eElO/u5mlhsODgVxgWbLjbuo2bSvj0fdXcubIPgztoaMVEYlOFKfCehI8E1Oz/sfd/WUzWwJkEJwag523FR8H/NLMKoBq4Gp3L4wg7iZtwtvLKKus4voTdW1FRKKV9MTi7suAkXWU13n+xt2fAZ5p7Lias43bynj0g5qjlfZRhyMirZz6vG8BHnhrKWWVem5FRJoGJZZmrqC4lH9OXcnZo/oyuLuOVkQkekoszdyEt5ZRXlnN9TpaEZEmQomlGSsoLuVf01Zy9iF9GdQtK+pwREQAJZZm7a9vLqOiyvmu7gQTkSZEiaWZKthaymPTVnLOIX0ZqKMVEWlClFiaqfvfWkpltXP9iXrKXkSaFiWWZmjN5hIem7aKcw/py4CuOloRkaZFiaUZumvypxhw05e+0GSaiEjklFiambn5RTw7M58rjxlEn85tow5HROQLlFiaEXfnty8toEu7NnznhCFRhyMiUicllmbk7cUbeW/JJq4/MZeOmW2iDkdEpE5KLM1EVbXz2xcXkJPdjkuOGBB1OCIiu6XE0kw8OzOfheuKufmU4aSn6WMTkaZL31DNQGlFFXdNXsTIfp348kG9ow5HRGSPlFiagX+8t4K1RaXcevr+pKRY1OGIiOyREksTV7i9nL+8sYST9u/BEYO7Rh2OiMheRZJYzGyFmc0xs1lmNj0syzazKWa2OPzbJSw3M7vHzJaY2WwzGx1FzFG59/XFbC+v5Een7hd1KCIicYnyiGWsu49y97xw/BbgNXfPBV4LxwFOA3LD13jg/qRHGpGVm7bzr6kr+dqh/cnt2SHqcERE4tKUToWdBTwSDj8CnB1T/qgHpgKdzaxVXMH+/SuLSEtJ4caT1HSLiDQfUSUWByab2QwzGx+W9XT3teHwOqBnONwXWB0z75qwrEWbtXoLz89ey7ePHUTPjplRhyMiEre0iNZ7jLvnm1kPYIqZLYyd6O5uZl6fBYYJajxATk5O4iKNgHvwMGTXrHTGH6+mW0SkeYnkiMXd88O/BcCzwGHA+ppTXOHfgrB6PtA/ZvZ+YVntZU5w9zx3z+vevXtjht/oXl9YwLTlhdx4Ui7tM6LK/SIiDZP0xGJmWWbWoWYYOBmYC0wCLgurXQY8Fw5PAr4R3h12BFAUc8qsxamsquaOlxYyuFsWFx7WvI+8RKR1iuLncE/gWTOrWf/j7v6ymX0EPGVm3wRWAheE9V8ETgeWACXAFckPOXn+M2MNiwu28ddLxtAmtSndWyEiEp+kJxZ3XwaMrKN8EzAQo7VgAAARwklEQVSujnIHrk1CaJErKa/kj1M+ZcyALpwyoufeZxARaYL0k7gJefCd5WwoLuPHp+9HeEQnItLsKLE0ERuKy3jgraWcOqIXYwZkRx2OiEiDKbE0Ae7O7ZPmUVZZzQ9PHR51OCIi+0SJpQl48qPVvDBnLd8/eTiDu7ePOhwRkX2ixBKxJQXF/OJ/8zl6aFeuOm5w1OGIiOwzJZYIlVZUcf0Ts2ibnsofLxilvlZEpEXQY90RuuOlhSxYu5WHLs9Te2Ai0mLoiCUiry1Yz8Pvr+CKowdy4n56ZkVEWg4llgis31rKzU/PZv/eHbnlNHXgJSItixJLklVVOzc9OYsd5VXc+/VDyEhLjTokEZGE0jWWJHvg7aW8v3QTd371IIb20K3FItLy6IgliT5etZm7Jn/Klw/uzQV5/fc+g4hIM6TEkiRbSyu44d8z6dUxk9+cc5DaAhORFkunwpLA3fnps3P5bEspT111BJ3atok6JBGRRqMjliR45uN8Jn3yGTeOy1UDkyLS4imxNLJlG7Zx23NzOXxQNteMHRp1OCIijU6JpRGVVVbx3X/PJD0thbsvHEWqmmwRkVYgssRiZqlmNtPMng/H3zGzWeHrMzObGJafYGZFMdNuiyrm+nB37nxpEXPzt/K7rx5M705tow5JRCQporx4fwOwAOgI4O7H1kwws2eA52LqvuPuZyQ3vIbbXlbJTyfO5dmZ+XzjyAGcPKJX1CGJiCRNJEcsZtYP+DLwYB3TOgInAhOTHVciLFpXzJn3vctzs/L53peG8fOvjIg6JBGRpIrqiOVu4IdAhzqmnQ285u5bY8qONLNPgM+AH7j7vNozmdl4YDxATk5O4iOOw1PTV3Pbc3Npn9GGf33rcI4a0i2SOEREopT0xGJmZwAF7j7DzE6oo8rX2fVI5mNggLtvM7PTCY5kcmvP5O4TgAkAeXl5nvDA96CkvJKfTZzHMx+v4aghXbn7wlH06KBm8EWkdYriiOVo4MwwSWQCHc3sX+5+iZl1Aw4DzqmpHHvk4u4vmtlfzKybu29MeuR1WLy+mGse+5glG7Zxw7hcvjsuV3d/iUirlvTE4u63ArdCcMcXwamtS8LJ5wHPu3tpTX0z6wWsd3c3s8MIrgttSm7UdXtmxhp+OnEuWRmp/PPKwzkmV6e+RESaWpMuFwJ31Co7D/iOmVUCO4AL3T2pp7pq21Fexc8nzeWp6Ws4fFA293z9EPUAKSISijSxuPubwJsx4yfUUec+4L6kBbUXSwq2ce1jH/NpQTHXjR3KjSflkpaq50xFRGo0tSOWJqlgaynTV27moxWFPPnRajLbpPLwFYdx/LDuUYcmItLkKLHUUl3tfFpQzPQVm5mxcjPTVxayunAHABlpKRyb241fn30QvTrp1JeISF1abWJxdyqqnNLKKubmFzFjxWamr9zMx6s2U1xaCUC39hnkDejCZUcOZMyALozo04n0NJ32EhHZkxaZWD5dX8xxv3uDqmqnoqr687+V1R68qqqpruPy/7Ce7Tnj4D7kDehC3sAu5GS3U4dcIiL11CITS2abVEbndCYtNYW0FCMt1UhLqRneWdYmHB7WswOjc7rQqZ064BIR2VctMrHkZLfj7gsPiToMEZFWSRcMREQkoZRYREQkoZRYREQkoZRYREQkoZRYREQkoZRYREQkoZRYREQkoZRYREQkoSzirk0ahZkVA4uijiMO3YAm0RPmXijOxGoOcTaHGEFxJtpwd++wrwtpkU/eA4vcPS/qIPbGzKYrzsRRnInTHGIExZloZjY9EcvRqTAREUkoJRYREUmolppYJkQdQJwUZ2IpzsRpDjGC4ky0hMTZIi/ei4hIdFrqEYuIiESk2SUWMzvVzBaZ2RIzu6WO6Rlm9mQ4fZqZDYyZdmtYvsjMTokwxu+Z2Xwzm21mr5nZgJhpVWY2K3xNaqwY44zzcjPbEBPPt2KmXWZmi8PXZRHH+aeYGD81sy0x05K5PR8yswIzm7ub6WZm94TvY7aZjY6ZlpTtGUeMF4exzTGz981sZMy0FWH5rETdPbQPcZ5gZkUxn+1tMdP2uL8kOc6bY2KcG+6P2eG0ZG7P/mb2Rvi9M8/MbqijTuL2T3dvNi8gFVgKDAbSgU+AA2rVuQb4azh8IfBkOHxAWD8DGBQuJzWiGMcC7cLh79TEGI5va0Lb8nLgvjrmzQaWhX+7hMNdooqzVv3rgYeSvT3DdR0HjAbm7mb66cBLgAFHANMi2J57i/GomnUDp9XEGI6vALo1kW15AvD8vu4vjR1nrbpfAV6PaHv2BkaHwx2AT+v4f0/Y/tncjlgOA5a4+zJ3Lwf+DZxVq85ZwCPh8NPAODOzsPzf7l7m7suBJeHykh6ju7/h7iXh6FSgXyPEsTfxbMvdOQWY4u6F7r4ZmAKc2kTi/DrwRCPFskfu/jZQuIcqZwGPemAq0NnMepPE7bm3GN39/TAGiG7fjGdb7s6+7Nf1Vs84o9w317r7x+FwMbAA6FurWsL2z+aWWPoCq2PG1/DFjfN5HXevBIqArnHOm6wYY32T4FdCjUwzm25mU83s7EaIr0a8cX41PCx+2sz613PeRIh7XeEpxUHA6zHFydqe8djde0nm9qyP2vumA5PNbIaZjY8oplhHmtknZvaSmY0Iy5rktjSzdgRfxs/EFEeyPS24PHAIMK3WpITtny31yftmwcwuAfKA42OKB7h7vpkNBl43sznuvjSaCPkf8IS7l5nZVQRHgidGFEs8LgSedveqmLKmtD2bDTMbS5BYjokpPibclj2AKWa2MPzFHoWPCT7bbWZ2OjARyI0olnh8BXjP3WOPbpK+Pc2sPUFyu9HdtzbWeprbEUs+0D9mvF9YVmcdM0sDOgGb4pw3WTFiZicBPwHOdPeymnJ3zw//LgPeJPhl0Rj2Gqe7b4qJ7UFgTLzzJjPOGBdS61RDErdnPHb3XpK5PffKzA4m+LzPcvdNNeUx27IAeJbGOZUcF3ff6u7bwuEXgTZm1o0mti1j7GnfTMr2NLM2BEnlMXf/bx1VErd/JuPCUQIvQKURXDgaxM4LcyNq1bmWXS/ePxUOj2DXi/fLaJyL9/HEeAjBBcbcWuVdgIxwuBuwmEa68BhnnL1jhs8BpvrOi3nLw3i7hMPZUcUZ1tuP4GKoRbE9Y9Y5kN1fcP4yu14c/TDZ2zOOGHMIrj8eVas8C+gQM/w+cGqE27JXzWdN8IW8Ktyuce0vyYoznN6J4DpMVlTbM9w2jwJ376FOwvbPRtvYjbiBTie4o2Ep8JOw7JcEv/wBMoH/hP8cHwKDY+b9STjfIuC0CGN8FVgPzApfk8Lyo4A54T/DHOCbEW/L3wLzwnjeAPaLmffKcBsvAa6IMs5w/HbgjlrzJXt7PgGsBSoIzkN/E7gauDqcbsD/he9jDpCX7O0ZR4wPAptj9s3pYfngcDt+Eu4TP4l4W14Xs29OJSYR1rW/RBVnWOdyghuHYudL9vY8huCazuyYz/b0xto/9eS9iIgkVHO7xiIiIk2cEouIiCSUEouIiCSUEouIiCSUEouIiCSUEksLZ2a3m9kPwuFfhg9m7q7u2WZ2QPKi22XdJ5jZUUle535hy7IzzWxII67nJTPrV6tsj59FWCeh28TMBprZRQlc3l7fwx7mfdjMzguHH9zTfmdBK9t9GhqnJJ8SSyvi7re5+6t7qHI2QSvQUTiB4LmTpDCzVIL3+7S7H+KN1MyLmbUFurr7mtjyOD4LaMA2CVub2J2BQMISS5zvIZ7lfMvd5++hyuWAEktz0pgP5egVzYvgQdBPgXcJHuD6QVj+MHBeOHwHMJ/ggak/EHyBFRI8VTsLGAJ8G/iI4CGuZ9jZ1P/DwD0ETwsvq1lmOO1H7Hwo8Y6wbAjwMjADeIeYBy3D6QOBdQTNRMwCjg3LXg/jew3IqeN9Hs/Oh71mEjQHfgIxzakD9wGXh8MrgDsJ2pm6KGadb4TTJ4YxzgPGxyzj1HCeT4DXwrIs4CGCh3BnEjR/UtdncRrwuzrKYz+LFcAvwnXMIWhFoK5t0j38HD4KX0eH898O/BN4L/y8B4bb+ePwdVRYbypBo6yzgJsIHib+R7jOmcDYsN7l4baYEsZ2HfC9sM5Uwqeua72HQwn2h0/CbdKh1vu18LNYRPCA8Isx875J0GZearjMuWFMNwHnAdvC+WYBbYHbwvc/l6ArXYtZzp3h+j8Fjg3LUwn28bkE+9P1YfkY4K3wM3+FmJYm9NrH76CoA9ArwR9o8M8yB2gHdCR4UnaXxELQ2vOimH/IzrHTY5bVNWb41zH/kA8TtG6QQnCEsyQsPy38cqlJQDVfQK8RNl8DHE5MnxQxy7+9Js5w/H/AZeHwlcDEOub5Hzu/XNsTNOdxAntOLD/cwzpr4m0bfgl1JfgyXw0MqlXnN8AlNdsv/CLLqiPGe4AT6yj/fFuHcdVs22uAB3cT3+MEDRdC0PTKgph6M4C24Xg7IDMczmXn0/O1t833CfuuIUhmqwiSzeUE+02H8P0XsfPp7D8RNGD4+XsgaDplGXBoWN4RSKv1fs8lSFSpBEcfW/hiYhlD0Dx7zTydY6fX/pzC4X8CX4mpd1c4fDrwajj8HYIuNNJq5gfaEOyr3cOyrxHTj49e+/ZS68Ytz7HAsx7292J195pYBJQCfzez54Hnd7OsA83s1wRfnO0JftXVmOju1cB8M+sZlp0E/KNm3e5eGLamehTwn6BbHCBor21vjiT4MoLgy+N3ddR5D/ijmT0G/Nfd18SsY3ee3MO075rZOeFwf4Iv5e7A2x704YPvbJ32ZODMmutXBF/IOQT9XMQ6GvgBe1fTKOAMdr7v2k4CDoh5jx3D7QtBs0A7wuE2wH1mNgqoAobtZnnHAPcCuPtCM1sZU/cND/rtKDazIoIkDsGPloNrLWc4sNbdPwqXVVeruccRtJRdBXxmZq/XUWcZMNjM7gVeACbvJu6xZvZDggSaTXCEWRNf7HYcGA6fRNB+YGUYX6GZHQgcSNCqMAQJb+1u1if1pMTSCrl7pZkdBowj+MV5HXU3h/8wcLa7f2JmlxP84q1RFjO8p2/zFGCLu4/al5jr4u53mNkLBL9O37Ogu+lKdr12mFlrtu11LcvMTiD4AjrS3UvM7M065t1lFuCr7r5otxWCpvpXe9Dh1N7UbM8qdv9/mQIc4e6ltdYDu76vmwjaohsZzrNL/TjFfr7VMePVe4hvn7j7Zgu6Qj6FoA2rCwiOVj9nZpnAXwiOYFab2e3s+jnFsx0h+PzmufuRCQpfYujifcvzNnC2mbU1sw4E/UDsIvyV28mD5sZvIvgCAigmOP1RowOwNmxu++I41j0FuCLs1Agzyw5/vS43s/PDMrOYftRj1F73+wStUxOu+5063scQd5/j7ncSnHPfD1hJ8Ks+w8w6EyTPeHQCNodJZT+C1l0huKZwnJkNqnlPYfkrwPVh76SYWV3N8Z9GcG2poWpvk8kEXS8TrnN3yboTwRFENXApwa/xupb3DuHnambDCI64dpso92AR0NvMDg2X1aGOmwjeBr5mZqkW9Eo4tvZCwmbvU9z9GeCnBF3+1o67JolsDPfj8+KIbwpwVU1M4We4COhuZkeGZW1iOguTfaTE0sJ40P3okwQXUV8i+MKtrQPwvJnNJrjA/72w/N/AzTG33/6MoJe594CFcaz7ZWASMN3MZrHzFNDFwDfNrKYl17q6iv0fcE54+++xBF+gV4QxXgrcUMc8N5rZ3LBOBfCSu68GniK4RvIUwQXneLwMpJnZAoIbG6aG72kDMB74bxh/zam0XxGccpptZvPC8dpOZd8SS+1t8l0gz4IePecT/Kqvy1+Ay8J492Pn0cxsoMqCXhdvCuulmNmc8H1d7jF9A8UrPCL7GnBvuM4pfPFo71mCbgvmEzTf/kEdi+oLvBnuO/8Cbg3LHwb+GpaXAX8j+Hxfoe79u7YHCa4fzQ7juyiM+TzgzrBsFkm8K7GlU+vGIo3AzDIIegzMizoWkWRTYhERkYTSqTAREUkoJRYREUkoJRYREUkoJRYREUkoJRYREUkoJRYREUkoJRYREUmo/w/qMnBBdSXM7gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10b7945c0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot results\n",
    "plt.plot(distance_surface / distance_atoms, np.abs(C6))\n",
    "plt.xlim(min(distance_surface / distance_atoms), max(distance_surface / distance_atoms))\n",
    "plt.xlabel(\"distance to surface / interatomic distance\")\n",
    "plt.ylabel(r\"|C$_6$| (GHz $\\mu m^6$)\")\n",
    "plt.show()"
   ]
  }
 ],
 "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.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}