{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "### Stretched LDA Li_2" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from CADMium import Pssolver, Psgrid, Partition, Inverter\n", "import CADMium\n", "\n", "from copy import copy" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "----> Begin SCF calculation for *Isolated* Fragments\n", "\n", " Total Energy (a.u.) Inversion \n", "\n", " __________________ ____________________________________ \n", "\n", "Iteration A B iters optimality res \n", "\n", "___________________________________________________________________________________________ \n", "\n", " 1 -8.63010 -8.63010 1.000e+00 \n", " 2 -7.59689 -7.59689 1.538e-01 \n", " 3 -7.38753 -7.38753 5.774e-02 \n", " 4 -7.34922 -7.34922 2.769e-02 \n", " 5 -7.34315 -7.34315 1.305e-02 \n", " 6 -7.34281 -7.34281 6.081e-03 \n", " 7 -7.34335 -7.34335 2.619e-03 \n", " 8 -7.34343 -7.34343 1.427e-03 \n", " 9 -7.34358 -7.34358 6.655e-04 \n", " 10 -7.34365 -7.34365 3.113e-04 \n", " 11 -7.34368 -7.34368 1.462e-04 \n", " 12 -7.34369 -7.34369 6.884e-05 \n", " 13 -7.34369 -7.34369 3.256e-05 \n", " 14 -7.34369 -7.34369 1.547e-05 \n", " 15 -7.34369 -7.34369 7.381e-06 \n", "----> Begin SCF calculation for *Interacting* Fragments\n", "\n", " Total Energy (a.u.) Inversion \n", "\n", " __________________ ____________________________________ \n", "\n", "Iteration A B iters optimality res \n", "\n", "___________________________________________________________________________________________ \n", "\n", " 1 -7.32464 -7.32464 10 +4.650e-15 +1.000e+00\n", " 2 -7.32866 -7.32866 11 +5.560e-15 +4.569e-02\n", " 3 -7.34060 -7.34060 9 +6.817e-15 +2.899e-02\n", " 4 -7.34060 -7.34060 10 +6.045e-15 +2.343e-02\n", " 5 -7.33583 -7.33583 8 +5.524e-14 +9.256e-03\n", " 6 -7.33571 -7.33571 7 +4.544e-15 +1.038e-02\n", " 7 -7.33833 -7.33833 7 +3.567e-15 +1.726e-03\n", " 8 -7.33903 -7.33903 6 +5.195e-15 +4.218e-03\n", " 9 -7.33788 -7.33788 6 +1.915e-14 +1.499e-03\n", " 10 -7.33717 -7.33717 6 +4.731e-15 +1.894e-03\n", " 11 -7.33751 -7.33751 6 +4.848e-15 +7.625e-04\n", " 12 -7.33798 -7.33798 5 +6.510e-15 +6.076e-04\n", " 13 -7.33800 -7.33800 5 +4.045e-15 +5.146e-04\n", " 14 -7.33779 -7.33779 5 +6.872e-15 +1.337e-04\n", " 15 -7.33769 -7.33769 4 +3.152e-14 +2.090e-04\n", " 16 -7.33774 -7.33774 4 +4.111e-14 +4.090e-05\n", " 17 -7.33779 -7.33779 4 +4.975e-15 +7.323e-05\n", " 18 -7.33781 -7.33781 4 +7.882e-15 +2.652e-05\n", " 19 -7.33779 -7.33779 4 +7.553e-15 +2.497e-05\n", " 20 -7.33778 -7.33778 4 +5.466e-15 +1.356e-05\n", " 21 -7.33778 -7.33778 3 +4.740e-15 +8.645e-06\n", "Done with 5.122\n" ] } ], "source": [ "# dis_eq = np.linspace(1.0,5,30)\n", "# dis_st = np.linspace(5.1,10,10)\n", "# dis_eq = np.linspace(1.0,5,10)\n", "# dis_st = np.linspace(5.1,10,3)\n", "# distances = np.concatenate((dis_eq, dis_st))\n", "distances = [5.122] \n", "# distances = [2.0]\n", "energy = []\n", "\n", "for d in distances:\n", " a = d/2\n", " Za, Zb = 3,3\n", " pol = 2\n", "\n", " #Set up grid\n", " NP = 7\n", " NM = [6,6]\n", " L = np.arccosh(15/a)\n", " loc = np.array(range(-4,5)) #Stencil outline\n", " grid = Psgrid(NP, NM, a, L, loc)\n", " grid.initialize()\n", "\n", "\n", " #Fragment a electrons [alpha, beta]\n", "\n", " #Fragment a electrons [alpha, beta]\n", " Nmo_a = [[1,2]]; Nmo_A = [[2,1]] #Number of molecular orbitals to calculate\n", " N_a = [[1,2]]; N_A = [[2,1]]\n", " nu_a = 0.5\n", "\n", " #Fragment b electrons\n", " Nmo_b = [[1,2]]; Nmo_B = [[2,1]]\n", " N_b = [[1,2]]; N_B = [[2,1]] \n", " nu_b = 0.5\n", "\n", " #Molecular elctron configuration\n", " Nmo_m = [[3,3]]\n", " N_m = [[3,3]]\n", "\n", "\n", "\n", " part = Partition(grid, Za, Zb, pol, [Nmo_a, Nmo_A], [N_a, N_A], nu_a, \n", " [Nmo_b, Nmo_B], [N_b, N_B], nu_b, { \"AB_SYM\" : True,\n", " \"interaction_type\" : \"dft\", \n", " \"kinetic_part_type\" : \"inversion\",\n", " \"hxc_part_type\" : \"exact\",\n", "# \"k_family\" : \"gga\", \n", "# \"ke_func_id\" : 500,\n", " })\n", "\n", " #Setup inverter object\n", " mol_solver = Pssolver(grid, Nmo_m, N_m)\n", " part.inverter = Inverter(grid, mol_solver, { \"AB_SYM\" : True, \n", " \"use_iterative\" : False,\n", " \"invert_type\" : \"orbitalinvert\",\n", " \"DISP\" : False, \n", " })\n", "\n", " part.optPartition.isolated = True\n", " part.scf({\"disp\" : True,\n", " \"alpha\" : [0.6],\n", " \"e_tol\" : 1e-5})\n", " \n", "# atom = copy(part.E.Ea)\n", " \n", "# part.KSa.solver[0,0].optSolver.iter_lin_solver = False\n", "# part.KSa.solver[0,1].optSolver.iter_lin_solver = False\n", " \n", "# part.KSA.solver[0,0].optSolver.iter_lin_solver = False\n", "# part.KSA.solver[0,1].optSolver.iter_lin_solver = False\n", " \n", "# part.KSb.solver[0,0].optSolver.iter_lin_solver = False\n", "# part.KSb.solver[0,1].optSolver.iter_lin_solver = False\n", " \n", "# part.KSB.solver[0,0].optSolver.iter_lin_solver = False\n", "# part.KSB.solver[0,1].optSolver.iter_lin_solver = False\n", "\n", " part.optPartition.isolated = False\n", " part.scf({\"disp\" : True,\n", " \"alpha\" : [0.6],\n", " \"max_iter\" : 50,\n", " \"e_tol\" : 1e-5,\n", " \"iterative\" : False,\n", " \"continuing\" : True})\n", "\n", " energy.append(part.E.E)\n", " print(f\"Done with {d}\")\n", " \n", " \n", "# energy = np.array(energy)\n", "# np.save('h2plus_distance.npy', distances)\n", "# np.save('h2plus_overlap.npy', energy)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'Ea': -7.337782105762022,\n", " 'Eb': -7.337782105762022,\n", " 'Ef': -14.675564211524044,\n", " 'Tsf': 14.51597413441419,\n", " 'Eksf': array([[-3.74685652, -3.74685652]]),\n", " 'Enucf': -33.888205438333586,\n", " 'Exf': -3.041843099343745,\n", " 'Ecf': -0.30099395411366964,\n", " 'Ehf': 8.039504145852767,\n", " 'Vhxcf': 11.684638443813574,\n", " 'Ep': -1.806931800537365,\n", " 'Ep_pot': -3.6782187039973078,\n", " 'Ep_kin': 0.004774410724220246,\n", " 'Ep_hxc': 1.8665124927357226,\n", " 'Et': -16.48249601206141,\n", " 'Vnn': 1.757126122608356,\n", " 'E': -14.725369889453054,\n", " 'evals_a': array([-1.81050752, -1.81778993, -0.11855907, -1.81778993, -0.11855907,\n", " -1.81050752]),\n", " 'evals_b': array([-1.81050752, -1.81778993, -0.11855907, -1.81778993, -0.11855907,\n", " -1.81050752]),\n", " 'Ep_h': 1.8866167323214942,\n", " 'Ep_x': 0.0070897811354484475,\n", " 'Ep_c': -0.027194020721220125}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "vars(part.E)" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[6.63782316e+00, 6.63782316e+00],\n", " [5.94235124e+00, 5.94235124e+00],\n", " [4.77075119e+00, 4.77075119e+00],\n", " ...,\n", " [1.61342469e-08, 1.61342469e-08],\n", " [1.65561742e-08, 1.65561742e-08],\n", " [1.67729863e-08, 1.67729863e-08]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "part.nf" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "No handles with labels found to put in legend.\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAD4CAYAAAAZ1BptAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAy+ElEQVR4nO3deZTc5Xng++9T1bX03upVvWkBWhghgQAZsAFjB7ABOxaOgyMyY+MED2EG7swkOXeMx/fOzZybzOU6wZ7ruQSGTBjLJwlL7DjoBmEM2EBsbGsBoRUhIaTuVjdSq6VuSb1U1/LcP+r3qy61qnqrvfv5nNOnqn5L1VvSr+qp933eRVQVY4wxJhVPoQtgjDGmeFmQMMYYk5YFCWOMMWlZkDDGGJOWBQljjDFplRW6ANnU2NioK1asKHQxjDGmpOzYseOkqjal2reggsSKFSvYvn17oYthjDElRUSOpttnzU3GGGPSsiBhjDEmLQsSxhhj0lpQOQljjFnswuEwvb29jI+PX7AvGAzS0dGBz+eb9fNZkDDGmAWkt7eX6upqVqxYgYgktqsqg4OD9Pb2snLlylk/nzU3GWPMAjI+Pk5DQ8N5AQJARGhoaEhZw5iOBQljjFlgpgaImbZPx4KEWfQi0RjPbusmFIkWtBz9w2P8eE9/QctgzFQWJMyi98r+43z9h7v5xaGTBS3HM1t7eOBv3uLY0FhBy2FMMgsSZtH7yb7jAIyECluTGAlFAHjFKY8x85VuMbn5LDJnQcIUnSffeJ8vPv7mvC7ouYpEY/z03RMAjIULGyTc139lf36CxPEz49zy6Gv8zHn/ZmEIBoMMDg5e8PlxezcFg8E5PZ91gTVF5Zmt3fyXLe8CMDgyQWNVIKevt+PoaYZGwwCMF0mQ+NXhQc6Mh6kJzr4v+3zsOTbM+wMj/Ou/3cHT/+p6rlq2JKevZ/Kjo6OD3t5eBgYGLtjnjpOYCwsSpmi8su84//FHu2muDnDibIj+ofGcB4mX9x3H6xGiMS14kAiFY3g9QjiqvHZggM9f2ZbT1+sbjneFrAn6+P3vbeMH//rjXNxUldPXNLnn8/nmNA5iJtbcZIrCjqOnePDv3mJNey3fvecqgJwncFWVl/cf5+MXNwAwNhHL6evNZCwc5dKWahoq/XnJS/QNjVHmEZ65/3q8HuErf72V42fm1ofeLHwWJEzBnRkPc9+m7bTWBnnqqx+lqzn+a7Z/OLdB4uCJcxwdHOUzly/F7/UUPicxEaUy4OU3PtLMzw6cIBzNbdDqHxpjaW2Qi5qq+J9fvZbToxM88Dc7cvqapvRYkDAFd3hghKHRMN+48zIaqwLUV/oJlHnoy3FN4mXn1/qtl7UQ9HkK3tw0Fo4S9Hm5bXULZ8cjbP3gVE5fr29onLbacgDWdtTy1Y+vYGfPEJEcBydTWixImIIbPBcCoKUm3utCRGirK0+0mefKy/uOc0VHLUtrg5T7vYxNFDZIjIejlPu83NjVSKDMkwhiudI3PEZb3WRPl5aaIKpw2knkGwMWJEwRGDw3AUBDpT+xra0umNOaxHg4yju9Q9y8Kr5iY9DnZbzAI67HnZpEhb+Ma1fW8+sc1iSiMeXD4XHa6soT2xqq4v/+gyOhnL2uKT0WJEzBDY44QaJqMki01pbTP5S7mkTv6VFUSfTmKfcVviYx5tQkIF6u7sGRnI0VOXkuRCSmtCYHicp4T7JTTtA2BixImCIweC5EhT/+C9rVVlfO8bPjOUveHh0cBWBZQwUQr0kUQ+K63B8PEssbKhiZiCYCaLa5Pcfak5qbGp0gfTJHr2lKkwUJU3CDIxPn1SIA2mrj7eO56pLpBonl9fEgUe7zFjxxPR6OEfRNBgmYLGe2ubW01trk5qZ4TcLNERkDWQoSInK7iBwQkUMi8nCK/SIi33X27xKRq53tl4rIzqS/MyLy7519fyIix5L23ZmNspric/JcKNHU4XLbyvty1OTUfWqUqkAZ9U4epNzvZTxcuF490ZgyEY0R9MU/ksvqK51yjuTk9dx8T3JOoq7ch0cmc0TGQBZGXIuIF3gMuA3oBbaJyGZV3Zd02B1Al/N3HfA4cJ2qHgDWJT3PMeBHSed9R1X/ItMymuI2eG6C1trz55Nxe93kaqzEkcERljdUJObXD/oKO07CrcW4OYnO+nJE4MjJ3NQk+obHqPR7qQlOfgV4PEJ9ZcAS1+Y82ahJXAscUtXDqjoBPANsmHLMBuD7GvcroE5EWqcccwvwvqoezUKZTAkZHAld0NzkNoPkatR19+BookkHnJxEARPXboBycxKBMi9tteV0n8pRkBgao62u/IJFaBoq/Zy0moRJko0g0Q70JD3udbbN9ZiNwNNTtj3kNE89JSIpZx8TkftFZLuIbE81oZUpbqrKqZGJRHu4qzJQRm25Lyc9nKIxpef0aKJJBwqfk3ADlJuTAFhWX8HRwdw0N/UPj5/Xs8nVUOW3nIQ5TzaCRKr18Kb225v2GBHxA58H/j5p/+PAxcSbo/qBR1O9uKo+qarrVXV9U1PTHIptisGZ8QjhqJ43RsLVVleek7ESfUNjhKPKiqSaRKGDhLsqXnlSkFjRWJGzxHXf0Nh5PZtcDVUBTlnvJpMkG0GiF+hMetwB9M3xmDuAt1Q1McRUVY+ralRVY8BfEW/WMguM+6s11WyvbbXBnIy6dptwlk1tbgpH87KGRSru5ILn1yQqGRyZ4JyzGFG2jIejnDw3cV7PJldDpd8S1+Y82QgS24AuEVnp1Ag2ApunHLMZ+IrTy+l6YFhVkxfzvYcpTU1TchZfAPZkoaymyKQaSOfKVU0i0f21Iam5ye8lpjBRoHmLxsIX1iQmu8Fmt8npQyfwtqVobmqs8nM2FCl4d2BTPDIOEqoaAR4CXgL2A8+p6l4ReUBEHnAO2wIcBg4RrxX8G/d8Eakg3jPqH6Y89bdEZLeI7AI+BfxhpmU1xcetSUztAgvQWhdkeCycWNYzW44OjuD3elhaM9nc4v6CHy/QdOGTievJj+Sy+tyMlehzeoy11aZubgKsyckkZGXRIVXdQjwQJG97Ium+Ag+mOXcUaEix/cvZKJspbm5PmsYUNYl255du//AYlzRXZ+01jw6O0llfjtczmSpzf8GPR6LUktsV4VJxf7kHU9YkshwkhtLXJNzc0OC5iZT7zeJjI65NQbnt30vSJK4h+wPqjp4aPa+pCSZ/wReqG2yqIFEd9NFQ6c/6gLp+pwlv6TQ1iZM2VsI4LEiYghocCVFb7sPnvfBSdAfYZTMvoap0D44kmnJcwbL4l3OhBtS5wSk5JwHx5Houmpsaq/znBSSXW6Oz5LVxWZAwBTV47sJ5m1wtNUE8QlZ7OJ08N8HIRPS87q8AQX+Bg0SKxDXAiobKnDQ3pWtKqk80N1lNwsRZkDAFNTgSojFF0hrA5/XQXJ3ddSXcppsLmpsSiesCBwn/lJpEfQV9w2OJcRTZ0Dc0dsE0KK6qQBn+Mo8lrk2CBQlTUNPVJCA+h1M252+aOkW4KzlxXQju5IKBsvM/kssbKlCF3tPZ+TdQ1cSUHKmICI02NYdJYkHCFFSqacKTtdaVZzVx/e6HZ/GXeehcMqW5yQkSYwXqAhtflc5zwVxKlzTHF0Xa338mK69zZjzCyEQ0sbZ1Kg1VNsmfmWRBwhRMJBrj9OhEyjESrnZnQF22RkLvOHqaK9pr8U/5xe7WJAqZuJ6ajwC4rLWGoM/DW0eHsvI6bq1suu6t8fmbrCZh4ixImII5PRpGNfVoa1drbZBQJJaVNvJQJMru3mGuWX7hXJFBtwtsARPXqYKEz+vhyo46dnSfzsrruPmd1hTzNrkaKgOWuDYJFiRMwbhNGtPVJNxfvNlok99z7AwT0RhXpwgS7hd0qEBBIt7cdGGQALhm+RL2HhvOylQZ7r9j+zQ1icYqPydHJgo2j5UpLhYkTMG4TRrT1SQuW1oDwJ6+4Yxf762j8V/jVy9LUZNI5CSKM0hEYsqu3sz/DfYcG6a+0k9zdfrA3FDlZyISy/rEgqY0WZAwBXMyMQNs+iDRWV9OfaWfnd1DGb/ejqOnWd5QQVOKL0if10OZRwrb3ORPHSSucoLajqOZNznt7Bniyo7aCxLkydyaneUlDFiQMAXk5hmma24SEa7sqOWd3qGMXktV2dF9mmtS1CJc5c504YWQLnEN8QFuFzVVZhwkzoUiHDxxjis766Y9rt4ddW1jJQwWJEwBDZ6bwOsRasunn1BvXecSDp44x9nx8Lxfq/f0GANnQynzEa6g35sYr5Bv4+FY2uYmgGuWLeGt7tMZ5Ql29Q6hCutmCBKNiZqEJa+NBQlTQIMjIeor/Xg86Zs+AK7srEUVdh+bf5u8+ys8Vc8mVyFXp3PHSaRzzfIlnBqZ4EgGU3S80xP/97uyo27a4xqsJmGSWJAwBXPy3ETKZUuncn/57uwZmvdrvXFwgOpAGata0k85Xu7zFixxna4LrGv9inhw++eD81/HfWfPaVY0VKSccTeZzd9kklmQMAUzeC40bc8mV12FnxUNFbwzzyAxPBZmy+5+Pndl23lrSEwV9HmKMnENcHFTFZe31fD01p55Nzm90zM8Yz4C4j29qgNlNjWHASxImAIaHJl+tHWydZ11865J/OPbxxgPx/gX1y2b9rhggZubpqtJiAj3XLuM/f1neGceXWE/HB7nwzPjM+YjXA1VfmtuMkCWgoSI3C4iB0TkkIg8nGK/iMh3nf27ROTqpH1HnGVKd4rI9qTt9SLysogcdG7TNyabkjTT5H7J1nXWcfxMaM6T/akqT2/tZm17LWvaa6c9ttxfmCARiynj4RiBaYIEwIZ1bVT4vTz96+45v8bOnnhOZvZBwkZdm7iMg4SIeIHHgDuA1cA9IrJ6ymF3AF3O3/3A41P2f0pV16nq+qRtDwOvqmoX8Krz2CwQ4+Eo50IRGqtmV5Nwm0nm2uS09YNTvPvhWX53hloEFK4LbCgSS7z+dKqDPj5/ZRub3+mb8zQlO3uG8XmFy1prZnV8Q6Xfpgs3QHZqEtcCh1T1sKpOAM8AG6YcswH4vsb9CqgTkdYZnncDsMm5vwm4KwtlNUVicozE7GoSq9tq8HmFt+cQJGIx5c+27KelJsCGdW0zHh8sUJCYXHBo5o/jfTeuZCIa49GfHJjTa+zsOc3q1pppu9kma6gKWE7CANkJEu1AT9LjXmfbbI9R4CciskNE7k86pkVV+wGc2+YslNUUickpOWZXkwiUeVndWjOnmsQP3uplV+8w37jjMir8ZTMeH/R5CzJVeLoFh1Lpaqnmy9cv5+mt3bOePjwaU3b3zi5p7YrXJELEYjZ/02KXjSCRqrvI1CtrumNuUNWriTdJPSgin5jTi4vcLyLbRWT7wMD8uwea/DrpTu43y5wExNvTd/UOMzLLOYWe+vkHXN5WM6taBMSbewoxwZ+bB5ntr/w/vHUVQZ+XTW8emdXxe44NMzIRnXU+AuL/LzGFobH5D2A0C0M2gkQv0Jn0uAPom+0xqurengB+RLz5CuC42yTl3J5I9eKq+qSqrlfV9U1NTRm+FZMviZrELJubAD6/rp2xcJT/68X9KfeHIlGeeP19fv972/i/f/wu7354lt/5aOe08xQlK/cXpgusOzZjtkGitsLH7WuW8sKufr7143f52qZtPL/zWMqusZFojP+0eS91FT5uXjX7z4dbw7PktclGkNgGdInIShHxAxuBzVOO2Qx8xenldD0wrKr9IlIpItUAIlIJfBrYk3TOvc79e4Hns1BWUyTcL5/ZNjdBfNTx125cyd/8qvuCQWW7e4f57Hd/ziMvvstP3z3B46+9jwj85hWzq0VAvCYRiSnhaH6bnMYTOYnZBQmAu9a1czYU4S9fe59X9p/g3z2zk688tZUPh89fxe+/v3GYd3qG+NO71szp37rRCd6WlzAzN9TOQFUjIvIQ8BLgBZ5S1b0i8oCz/wlgC3AncAgYBX7POb0F+JHzS68M+DtV/bGz7xHgORG5D+gG7s60rKZ4DI5MECjzUDmLdvhkf/zpS/npuyf4Dz/YxWP/4mr8Xg9/++tuntveQ1NVgO/93kfpPjXKf3p+L6rMOLo4WTBpdTqfN39DiOaSk3C1JS0a9KN/83H2HBvmv2x5l9u+/Tpfu+kibl+zlJ5To/zXV97js1e08rk5BEtIqknYMqaLXsZBAkBVtxAPBMnbnki6r8CDKc47DFyZ5jkHgVuyUT5TfE6eC9FYFZh1U5Ar6PPy7S+t44uPv8lv/eWbAPi8wr+8bhl/dNul1Fb42PrBqXmVyQ0S4+EoNcHpJx3MJndSwbnUJAJlk8eubqvhqmVLuLGriUde3M93XnmP77zyHgBN1QH+zw1r5lymxPxNVpNY9LISJIyZq1Mjsx9IN9WVnXW8/Ec3c3jgHCMTUa5ZvuS8ldamrl89W+6X9HieeziNJRLXsy93IOlYv1PrWdlYyX//8nr29g1zeGCEqmAZa9pqE3MxzcWSCj8iNsmfsSBhCmTw3MS0iw3NZGVjJSsbK1Pu83nnVjtxuc09+U5ej88xcT312Km1scvbarm8bfrR5TPxeoT6Cr8lro3N3WQKIz653+wTqXMx33yC+0s+30FibB6J68A8a0tz0VDlt+YmY0HC5J+qcjKD5qaZzD9ITOYk8ml8Holrfx4S6/WVfktcGwsSJv/OhSJMRGJzGiMxF/NubvIVprkpkZMom32QmGvCfz7ik/xZTWKxsyBh8m5yIF1umpvm+yvb/SU/nueFh8bCUfxlnhlX6Mu3xko/Jy0nsehZkDB5NziPKTnmYr7NTYWqSYxPTL+WRKE0VAU4Mx6v9ZnFy4KEybt9/WcBWN6QundSpnzzTOoGC9jcVIxBYnlDBQDvfji7iQTNwmRBwuTd6wdOsKy+ghXOl1C2zTcnMZm4zve0HLE5jZHIlxsvaUQEXjtgE2cuZsV3ZZoFLRSJ8ub7g9y8qilnyVefJ8PBdAWoScxljES+NFQFWNtey2sHUs6taRYJCxImr7Z9cJrRiSifvDR3M/bONwHs8wpejyRmZc2X8XB0Tt1f8+mTq5rY2TPE0Kj1clqsbMS1yavXDpzA7/XwsYsbcvo6D37qYj516dzWqRKRgixhOjbPxPWjd19JNMeLAt18aTPf/ekh/vngSX7zyrlNEmgWBgsSJq9ee2+A6y6qn9VKcZn4Xz/zkXmdF/R58j+YLhKltnzuEwp+8ZqOHJTmfOs666ir8PHagQELEouUNTeZvOk9PcqhE+fmtPhNvhVineuxieLMSUB8Dqebupp4/b0BW8p0kbIgYfLm9ffivWRymY/IVLnPW4BpOWJFGyQAbl7VxMlzIfbNck1ts7BYkDB589qBAdrryrm4qarQRUmr3O/Ne+J6LByl3F+8H0W35ucGebO4FO+VaRaUiUiMNw+d5JOX5q7razYEy7wFGCdRnIPpXE3VAda011hX2EXKgoTJi+1HTjEyEeWTc+xxlG9Bf35zEqpatOMkkn1yVTNvdQ8xPBYudFFMnmUlSIjI7SJyQEQOicjDKfaLiHzX2b9LRK52tneKyM9EZL+I7BWRf5d0zp+IyDER2en83ZmNsprCeO29AXxeyXnX10yV57l3UygSQ3VuCw4Vws2XNhGNKT8/eLLQRTF5lnGQEBEv8BhwB7AauEdEVk857A6gy/m7H3jc2R4B/lhVLwOuBx6ccu53VHWd83feGtqmtLx+YICPrqinKlDcva7zPU5ifB4LDhXCVZ111ATLrMlpEcpGTeJa4JCqHlbVCeAZYMOUYzYA39e4XwF1ItKqqv2q+haAqp4F9gPtWSiTKSJ9Q2McOH62qHs1ucr9+e3d5OY/inXEtavM60l0hVW1rrCLSTaCRDvQk/S4lwu/6Gc8RkRWAFcBv07a/JDTPPWUiCxJ9eIicr+IbBeR7QMD1vuiGE12fS3ufARAoCy/vZsSCw4V4QR/U918aRMnzobY78ziaxaHbFyZqbqqTP2pMe0xIlIF/BD496rqdsZ+HLgYWAf0A4+menFVfVJV16vq+qam4v+lutioKq/uP05bbZCu5uLt+uqK1yTy17vJDUjF3twE8XmcAF7Zf7zAJTH5lI0G4l6gM+lxB9A322NExEc8QPytqv6De4CqJq5EEfkr4J+yUFaTJ6dHJvjhW708vbWb9wdG+P0bVhZ111dXuc/LRDRGJBqjLA/rSE/WJIo/SDTXBLlm+RK+/fJ7/PPBAX73umXcsaa1JMpu5i8bQWIb0CUiK4FjwEbgd6ccs5l409EzwHXAsKr2S/xb46+B/ar67eQT3JyF8/ALwJ4slNXkUCymbDtyiqe3drNlz4dMRGJcvayOP//tK/j8utKY9ycxXXgkRlUegkSpJK5df33vep7d1sPTW7v5w2ff4T//f/v44tUdbPxoJ10t1YUunsmBjIOEqkZE5CHgJcALPKWqe0XkAWf/E8AW4E7gEDAK/J5z+g3Al4HdIrLT2fYfnZ5M3xKRdcSbpY4Af5BpWU32RaIxtn5wihf3fMhLez/kxNkQ1YEyNn60k3uuXcZlrTWFLuKcuLmB8XA0Lz2xEkGiyBPXrroKP39w88X8q5su4peHB/m7X3ez6c0j/PXPP+CS5iruWLOUz1y+lMvbakqi5mhmlpVPgfOlvmXKtieS7ivwYIrzfk7qfAWq+uVslM1kXzSm/OLQSV7Y1c/L+49zamSCoM/DJ1c1c8fapdy2uiXns7zmSmIJ0zwlr0upuSmZxyPccEkjN1zSyMDZEC/u6efF3R/y2M8O8d9+eohl9RXcvmYpn7+yjTXttYUurslAaX6STUEcGxrj77f38Pfbezk2NEZVoIzf+Egzd6xZys2XNpVsYEjm/qLPVzfYUkpcp9NUHeArH1vBVz62gsFzIV7ed5wX93zI//zFBzz5xmEub6th40c7+fy69nlNiW4Kq/Q/1SanJiIxXt1/nGe29fDGwQFU4aauRr5x50e49bKWkvsFPBP3yzpfA+rGS7QmkU5DVYCN1y5j47XLGB4Ns/mdYzy9tYf//fm9/OkL+/ns2lZ+56OdXLuy3pqjSoQFCZPSweNneW57D//w1jEGRyZorQ3yv3zqEu5e30lnfUWhi5czk+tc56cbbKkMppuP2gofX/7YCr78sRXsOTbMM9u6ef7tPv7h7WOsbKzk7vUd/PbVHTTXBAtdVDMNCxIm4ex4mBd29fPs9h7e7h6izCPcclkzG69dxie6mvDOc+3oUhLIc00ikZMoK/7BdJlY017Ln7av5Zt3rmbL7n6e3dbDt358gEd/8h6furSJu9d38hsfacaXhx5lZm4sSCxyqsq2I6d5bnsPL+zqZywcpau5iv/ts5dx11XtNFYFCl3EvEo0N01E8vJ6Y+EoZR7Jy5iMYlDu9/LFazr44jUdHB44x9/v6OWHO3p5Zf8JGqv8/NbVHXxpfSeXlMDAy8XCgsQipKq80zvMlt39vLCrP5GEvuuqNu5e38lVnXWLtr3Y7QIbiuSruam415LIpYuaqvj67R/hj29bxevvDfDsth6e+nk82b26tYbPXtHKZ9e2sqKxstBFXdQsSCwSqsrOniG27O5ny+4POTY0hs8r3HhJI3902yruWLt0QfROypTb3BTKU04iFIkRKIF5m3KpzOvhlstauOWyFgbOhnh+5zFe2N3Pn790gD9/6UAiYNy5tpWVFjDyzr4VFrDxcJS3uk/z0/0neHHPZGC4qauJP7xtFbdd1kJthXVJTBYoc2sS+clJhMIxAmWLsyaRSlN1gK/ddBFfu+ki+obGnB81kwHjstYaPrt2KTevamZ1W82iyJMVmgWJBSQaU/b2DfOLQ4O8+f5Jtn5wilAklggMf3TbKm5d3WJ91afhdkXNV3NTKBJd9DWJdNrqys8LGC/u+ZAXdvXxFz95j7/4yXvUVfj42EUNfPySRm68pJEVDRWLtpk0lyxIlDBV5fDJEd58f5BfHDzJLw8PJpaXvLSlmt+9bhk3XNzIdRfVUx20wDAbkzWJPDY3WU1iRm115dx340ruu3ElJ86O8+ahQX5x6CS/OHSSF/d8GD+mNsjHL2nkhksa+PjFjbRY19qssCBRImIx5YPBEfYcG2Zv3xl29w6zt2+YM+PxXjhttUE+vbqFG7sa+djFDTRX2wdkPso8gkfyN+J6PBxNBCYzO83VQe66qp27rmpHVTkyOJoIGC/vO84PdvQC0FITYE1bLZe317KmrYa1HbUsrQlabWOOLEgUobGJKEcGR9jXd4Y9fcPsOTbMvr4zjDhTOPi9Hj7SWs3nrmxjbXst11/UYFXtLBERAmXePNckLEjMl4iwsrGSlY2V/MvrlxONKfv6zvDrDwbZ23eGPceG+dmBE8Sc1WsaKv2JoLGmvZaPLK2ms77CxmdMw4JEgQyPhTk6OMLRwdGk21GOnhrh+JlQ4rhyn5fVbTX89jUdzsVdS1dLlV3UORTweQjlqSYRisQsR5RFXo+wtqOWtR2TkwqOTkTY33+Wvc4Prt3HzvDkG4eJOJHD6xHa6oKsaKhkWX1F/LahIvF4IY6GnwsLElmmqpweDTNwNsSJs+OcOBNi4FyIE2fij3tOj9E9OMLp0fB55zVXB1jRUMlNXU2saKhgWUMlq1urWdlYZT048ixQ5slfTSIcJVC9uAYs5luFv4xrli/hmuWTKyCPh6O8d/wsBz48S/ep0cSPtRd29zM05bPZUhNgeX0l7UvKaa4O0OT8NVcHaa6J368OlC3YmrwFiWmoKmPhKMNj4fjfaHjy/liYM2NhBs5NMHA2xMDZ8fjtuRDh6IULxVf6vTTXBGmvK+fOta0sb6hgeUMlyxsqWFZfYWMUikjQl7/mpolIbMFM7ldKgj4vV3TUcUVH3QX7hkfDHD01wpHBUboH47dHB0fYduQUJ86GmEhxbQR9Hpqrg07wiP/VVwaoLS+jtsJHbfnkX41zWyodFuybCfjg5Aj/9ZX3LggAw2PhlF/4LpF4G2djVYDmmiBdLdVJF8nkBdNUHaAyDwvYmOyI1yTy19xkOYniUlvh44qK1AFEVTkzFmHg3LjTOhBKtBrEb0McPHGOXxw6mehUkk7Q5zkveNSW+6gJ+qgIeKn0l1HhL6My4D3/1u+lIjDl1l+GP4fXkH1zEf8193b3UOI/qq22PBHt66b8Ckj+NVAdKMNjTUELTqDMm8dZYK13UykRkXjNoMLHJc3TL9caicY4Mx4578fn1B+hya0TfUPjvDt+lrGJKCMTkTldgz6v8JtXtvHtL63L8B1eKCtBQkRuB/4f4suX/g9VfWTKfnH230l8+dKvqupb050rIvXAs8AK4suXfklVT2ejvFNdurSaN/7Dp3Lx1KYE5b8mURrNDmZuyrwe6iv91Ff653V+NKaMTkQYnYgyEppyOxFhNOTcOtu7WnIzKWLGQUJEvMBjwG1AL7BNRDar6r6kw+4Aupy/64DHgetmOPdh4FVVfUREHnYefz3T8hozk3jvJhtxbQrL6xGqg76CD4TNxtV5LXBIVQ+r6gTwDLBhyjEbgO9r3K+AOhFpneHcDcAm5/4m4K4slNWYGeVrnEQ0poSjas1Npqhl4+psB3qSHvc622ZzzHTntqhqP4Bz25zqxUXkfhHZLiLbBwYG5v0mjHEFfflpbnJ7yVjvJlPMshEkUmVup3YJSnfMbM6dlqo+qarrVXV9U1PTXE41JqV81STcQGQ1CVPMsnF19gKdSY87gL5ZHjPducedJimc2xNZKKsxMwqUefIyd5Pbe8US16aYZSNIbAO6RGSliPiBjcDmKcdsBr4icdcDw04T0nTnbgbude7fCzyfhbIaM6N8jbi2moQpBRn3blLViIg8BLxEvBvrU6q6V0QecPY/AWwh3v31EPEusL833bnOUz8CPCci9wHdwN2ZltWY2Qj4vHnp3eQGIuvdZIpZVsZJqOoW4oEgedsTSfcVeHC25zrbB4FbslE+Y+bCHSehqjmdjydkzU2mBNhPGGOmCPq8xJTELKG54jY3Ba0mYYqYXZ3GTJGv1ekSzU1WkzBFzIKEMVO4QSLXPZzc57fEtSlmdnUaM4X7yz5vNQlrbjJFzK5OY6Zwv7RzvTrdZBdYa24yxcuChDFT5C0nkejdZB9DU7zs6jRmioAvv81NNneTKWYWJIyZwhLXxkyyq9OYKfKeuLYgYYqYXZ3GTJHISeQhce31CGVe+xia4mVXpzFTuCOg85G4tlqEKXZ2hRozRT6bmyxImGJnV6gxUyTGSeR4dbpQJGo9m0zRsyBhzBRuTWI8x9OFj1tzkykBdoUaM8XkYLrc1yRstLUpdhYkjJlisndTHnISNm+TKXJ2hRozhYjgz8MSpta7yZQCu0KNScFdnS6XrLnJlIKMgoSI1IvIyyJy0Lldkua420XkgIgcEpGHk7b/uYi8KyK7RORHIlLnbF8hImMistP5eyLV8xqTK0GfNy9dYG1VOlPsMr1CHwZeVdUu4FXn8XlExAs8BtwBrAbuEZHVzu6XgTWqegXwHvCNpFPfV9V1zt8DGZbTmDkJlHnyMneT1SRMscs0SGwANjn3NwF3pTjmWuCQqh5W1QngGec8VPUnqhpxjvsV0JFheYzJikA+chI2mM6UgEyv0BZV7QdwbptTHNMO9CQ97nW2TfX7wItJj1eKyNsi8rqI3JSuACJyv4hsF5HtAwMDc38HxqQQKPNa7yZjgLKZDhCRV4ClKXZ9c5avISm26ZTX+CYQAf7W2dQPLFPVQRG5BvhHEblcVc9c8ESqTwJPAqxfv16n7jdmPgK+PCSurbnJlIAZg4Sq3ppun4gcF5FWVe0XkVbgRIrDeoHOpMcdQF/Sc9wLfA64RVXVec0QEHLu7xCR94FVwPaZ35IxmQuW5SdxbTUJU+wyvUI3A/c69+8Fnk9xzDagS0RWiogf2Oich4jcDnwd+LyqjroniEiTk/BGRC4CuoDDGZbVmFkL+Dw5nSpcVZ2chNUkTHHLNEg8AtwmIgeB25zHiEibiGwBcBLTDwEvAfuB51R1r3P+/wtUAy9P6er6CWCXiLwD/AB4QFVPZVhWY2Yt14lrW3DIlIoZm5umo6qDwC0ptvcBdyY93gJsSXHcJWme94fADzMpmzGZCOS4ucmChCkVdoUak0KgLLfNTW5SPGBThZsiZ0HCmBTivZtyWJMIW03ClAa7Qo1JIde9m9zntkWHTLGzIGFMCgFfbqflcJ/bahKm2NkVakwKgTIvkZgSieamNmGJa1Mq7Ao1JgX3y3siZ0HCrUlYc5MpbhYkjEkh16vTJWoSNuLaFDm7Qo1Jwe2amqvktfVuMqXCrlBjUnAXA8rVJH/u81rvJlPsLEgYk4KbKxjPVXOT1SRMibAr1JgUEjmJHNckLHFtip0FCWNScL+8c5aTsMS1KRF2hRqTgvvlnfPeTdbcZIqcXaHGpJDz5qZwFBHwe+0jaIqbXaHGpBDMdRfYSIxAmQeRVKv7GlM8LEgYk4Jbk8jV/E3jtr61KREWJIxJIR+Ja8tHmFKQ0VUqIvUi8rKIHHRul6Q57nYROSAih0Tk4aTtfyIix5ylS3eKyJ1J+77hHH9ARD6TSTmNmavJaTly1QU2Zj2bTEnI9Cp9GHhVVbuAV53H5xERL/AYcAewGrhHRFYnHfIdVV3n/G1xzlkNbAQuB24H/tJ5HmPyItG7KWc1CWtuMqUh0yCxAdjk3N8E3JXimGuBQ6p6WFUngGec82Z63mdUNaSqHwCHnOcxJi9y3twUtuYmUxoyvUpbVLUfwLltTnFMO9CT9LjX2eZ6SER2ichTSc1VM52TICL3i8h2Edk+MDAw3/dhzHm8HsHnldwlriNRm7fJlIQZg4SIvCIie1L8zVQbSDxFim3q3D4OXAysA/qBR2dxzvkbVZ9U1fWqur6pqWmWRTJmZoEcLmFqNQlTKspmOkBVb023T0SOi0irqvaLSCtwIsVhvUBn0uMOoM957uNJz/VXwD/NdI4x+RIo8+Rw7qYY1cEZP37GFFymP2U2A/c69+8Fnk9xzDagS0RWioifeEJ6M4ATWFxfAPYkPe9GEQmIyEqgC9iaYVmNmZNAmSeH03JY4tqUhkx/yjwCPCci9wHdwN0AItIG/A9VvVNVIyLyEPAS4AWeUtW9zvnfEpF1xJuSjgB/AKCqe0XkOWAfEAEeVNXcrUpvTAoBXw6bm6wLrCkRGQUJVR0EbkmxvQ+4M+nxFmBLiuO+PM1z/xnwZ5mUz5hM5LS5yXISpkTYVWpMGgGfN2eLDlnvJlMqLEgYk4bVJIyxIGFMWvEgkf2ahKpa4tqUDAsSxqQRKPPmpHdTJKbE1BYcMqXBrlJj0gj4ctPcZEuXmlJiV6kxaeSqucmdWdaam0wpsCBhTBrBHPVuGncCT9BqEqYE2FVqTBq56t1kNQlTSixIGJNGrib4S+QkLHFtSoBdpcakESjzMBGJoZpyAuJ5s8S1KSV2lRqTRq5Wp7PmJlNKLEgYk0auVqez5iZTSuwqNSYNt/dRKMur07mr3dncTaYUWJAwJg2rSRhjQcKYtNwv8Wx3g50MElaTMMXPgoQxabhBItsD6tygY72bTCmwq9SYNAK+HDU3ha25yZSOjK5SEakXkZdF5KBzuyTNcbeLyAEROSQiDydtf1ZEdjp/R0Rkp7N9hYiMJe17IpNyGjMfieambCeuI9YF1pSOTNe4fhh4VVUfcb78Hwa+nnyAiHiBx4DbgF5gm4hsVtV9qvo7Scc9Cgwnnfq+qq7LsHzGzFvQahLGZNzctAHY5NzfBNyV4phrgUOqelhVJ4BnnPMSRESALwFPZ1geY7Iml4lrv9eDxyNZfV5jciHTINGiqv0Azm1zimPagZ6kx73OtmQ3AcdV9WDStpUi8raIvC4iN6UrgIjcLyLbRWT7wMDA/N6FMSlMBonsJ66tFmFKxYzNTSLyCrA0xa5vzvI1Uv1cmjoZzj2cX4voB5ap6qCIXAP8o4hcrqpnLngi1SeBJwHWr1+f3Ul2zKKWSFxnvXdTzHo2mZIxY5BQ1VvT7ROR4yLSqqr9ItIKnEhxWC/QmfS4A+hLeo4y4LeAa5JeMwSEnPs7ROR9YBWwfabyGpMtOWtuCscsaW1KRqY/ZzYD9zr37wWeT3HMNqBLRFaKiB/Y6JznuhV4V1V73Q0i0uQkvBGRi4Au4HCGZTVmTnI1TmLcmptMCcn0Sn0EuE1EDhLvvfQIgIi0icgWAFWNAA8BLwH7gedUdW/Sc2zkwoT1J4BdIvIO8APgAVU9lWFZjZmTyd5NOahJ2LxNpkRk1AVWVQeBW1Js7wPuTHq8BdiS5jm+mmLbD4EfZlI2YzJV5hE8Yolrs7jZlWpMGiKSk9XpQpGYBQlTMuxKNWYaAZ8nMbV3toTCUWtuMiXDgoQx01haE+TI4GjWnk9VOXpqlNaaYNae05hcsiBhzDTWddbxTs9Q1ta5PjI4ytBomHXL6rLyfMbkmgUJY6axrrOO4bEwH5wcycrz7ew5nXheY0qBBQljpuH+4t/ZM5SV59vZPUSF38uqluqsPJ8xuWZBwphpdDVXU+n3Zi9I9Ayxtr0Wr03uZ0qEBQljpuH1CGs7arMSJMbDUfb1n7F8hCkpFiSMmcG6ziXs7z+TcVfYff1nCEeVqywfYUqIBQljZrCus45wVNnbd8EkxHOys3vIeb6UCzgaU5QsSBgzg6uylLze2TPE0pogS2ttjIQpHRYkjJlBS02Q1tpgVoKEdX01pcaChDGzsK6zjre7T8/7/MFzIbpPjVrS2pQcCxLGzMK6zjp6T49x8lxoXue7tRCrSZhSY0HCmFlwv9zd5PNc7ewZwiOwtr02e4UyJg8yWk/CmMVibUd8ANzOniFuXd1CLKb88vAgp0Ym8Jd5CJR5KPN4CEdjhCIxJqIxLmqsZI0TFHb2DLGqpZrKgH3kTGmxK9aYWajwl7GqpZpffzDIUz//gE2/PMLRWcwOe/WyOr56w0re6Rnis1e05qGkxmSXBQljZmldZx1Pb+1m25HTXLN8CX/86Uu5bGk1oUi89hCJxvCXeeJ/Xg8/P3SSTW8e4d8+/XbifGNKTUZBQkTqgWeBFcAR4EuqekEXEBF5CvgccEJV18zmfBH5BnAfEAX+raq+lElZjcnUvR9fTrnPy11XtXFFR92Mx3e1VHPvx1bw+nsDvHbgBLdfbjUJU3okk3nyReRbwClVfUREHgaWqOrXUxz3CeAc8P0pQSLl+SKyGngauBZoA14BVqnqtPMirF+/Xrdv3z7v92OMMYuRiOxQ1fWp9mXau2kDsMm5vwm4K9VBqvoGcGoO528AnlHVkKp+ABwiHjCMMcbkUaZBokVV+wGc2+Ysnd8O9CQd1+tsu4CI3C8i20Vk+8DAwBxf3hhjzHRmzEmIyCvA0hS7vpn94ky+bIptKdvFVPVJ4EmINzflsEzGGLPozBgkVPXWdPtE5LiItKpqv4i0Aifm+Prpzu8FOpOO6wD65vjcxhhjMpRpc9Nm4F7n/r3A81k6fzOwUUQCIrIS6AK2ZlhWY4wxc5RpkHgEuE1EDgK3OY8RkTYR2eIeJCJPA78ELhWRXhG5b7rzVXUv8BywD/gx8OBMPZuMMcZkX0ZdYIuNdYE1xpi5y2UXWGOMMQvYgqpJiMgAcLTQ5chQI3Cy0IXIkYX83sDeX6lbzO9vuao2pdqxoILEQiAi29NV+0rdQn5vYO+v1Nn7S82am4wxxqRlQcIYY0xaFiSKz5OFLkAOLeT3Bvb+Sp29vxQsJ2GMMSYtq0kYY4xJy4KEMcaYtCxIFAERuVtE9opITETWT9n3DRE5JCIHROQzhSpjtojIn4jIMRHZ6fzdWegyZYOI3O78Hx1yFtBaUETkiIjsdv7PSn5aAxF5SkROiMiepG31IvKyiBx0bpcUsoyZSPP+5vXZsyBRHPYAvwW8kbzRWaFvI3A5cDvwlyLizX/xsu47qrrO+dsy8+HFzfk/eQy4A1gN3OP83y00n3L+zxbCWILvEf9MJXsYeFVVu4BXncel6ntc+P5gHp89CxJFQFX3q+qBFLtshb7ScC1wSFUPq+oE8Azx/ztTpNKsljmrlTZLwTSrgc6ZBYniNusV+krMQyKyy6kSl2yVPslC/X9KpsBPRGSHiNxf6MLkSKYrbZaCOX/2LEjkiYi8IiJ7UvxN94tz1iv0FZMZ3uvjwMXAOqAfeLSQZc2Skvx/mqMbVPVq4k1qD4rIJwpdIDNn8/rszbgyncmO6Vb4m0ZJrtA32/cqIn8F/FOOi5MPJfn/NBeq2ufcnhCRHxFvYntj+rNKTqYrbRY1VT3u3p/LZ89qEsVtwa3Q53z4XF8gnrQvdduALhFZKSJ+4p0NNhe4TFkjIpUiUu3eBz7Nwvh/myrTlTaL2nw/e1aTKAIi8gXgvwFNwAsislNVP6Oqe0XEXaEvwsJYoe9bIrKOeHPMEeAPClqaLFDViIg8BLwEeIGnnNUVF4oW4EciAvHvjL9T1R8XtkiZcVbL/CTQKCK9wP9BfGXM55yVM7uBuwtXwsykeX+fnM9nz6blMMYYk5Y1NxljjEnLgoQxxpi0LEgYY4xJy4KEMcaYtCxIGGOMScuChDHGmLQsSBhjjEnr/weZoDiXkgXpawAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x,d1 = grid.axis_plot(part.nf[:,0])\n", "x,d2 = grid.axis_plot(part.nf[:,1])\n", "\n", "x,pot = grid.axis_plot(part.V.vp_pot[:,0])\n", "x,hxc = grid.axis_plot(part.V.vp_hxc[:,0])\n", "x,har = grid.axis_plot(part.V.vp_h[:,0])\n", "x, xc = grid.axis_plot(part.V.vp_x[:,0] + part.V.vp_x[:,1])\n", "x, vp = grid.axis_plot(part.V.vp[:,0])\n", "\n", "x, xca = grid.axis_plot(part.KSA.V.vx[:,0] + part.KSA.V.vc[:,0])\n", "x, xcb = grid.axis_plot(part.KSB.V.vx[:,0] + part.KSB.V.vc[:,0])\n", "\n", "# plt.plot(x, pot)\n", "# plt.plot(x, har)\n", "plt.plot(x, vp)\n", "\n", "plt.legend()" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 30, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuU0lEQVR4nO3de3xcdZ3/8ddnZpJM7kmba9MbvVEKtKUU5KKIIlJQFnAFwVVQVwGV9e7q/h7u/tTHug9/rruCLoroouKKyKpguYlcVESLvdDSe2na0iRtbk2aTG4zmcv398fMSaZpLpPMOZnb5/l49JFk5syc7zSZ857vXYwxKKWUyl2uVBdAKaVUamkQKKVUjtMgUEqpHKdBoJRSOU6DQCmlcpwn1QWYiaqqKrN48eJUF0MppTLKtm3bThhjqsfenpFBsHjxYrZu3ZrqYiilVEYRkaPj3a5NQ0opleM0CJRSKsdpECilVI7TIFBKqRynQaCUUjlOg0AppXKcBoFSSuU4DQKlYna29LDt6EnHz2OM4dHtLfQOBh0/l1KJ0CBQKuZfn9jHF3610/HzHOrs59O/eJVfb29x/FxKJUKDQKmY471DHOrspz8QcvQ8O5p7AWjt9Tt6HqUSpUGgFBCJGNp9foyBXS29jp5rV0sPoEGg0ocGgVJA9+AwwXB029adsQu1U16NBU27BoFKExoESgFtcRflncecqxEMhyLsbfVFz+nTIFDpQYNAKUaDoKGi0NEawWvtfQyHIjRUFNLm82OMcexcSiVKg0ApRj+dX7aimubuIYLhiCPnOXxiAIA3La9iOBThpA4hVWlAg0ApoN3nxyWwtLoYgD6/MyOHfEPRC//KulLg1CYppVJFg0ApoiN4akq9zCnOB0Yv2Hbz+aPPu6I2FgS+IUfOo9R0aBAoRbRGUFvupcybB4xesO3mGwqR5xYWV0VrHm29AUfOo9R02BIEIrJBRA6ISKOIfHGc+1eKyCYRCYjI58bc97qI7BKRHSKi+0+qlGjr9VNXVkBZYSwIhhxqGvIHKfPmUVNagEt05JBKD0nvWSwibuBe4EqgBdgiIhuNMXvjDusGPgFcP8HTvMUYcyLZsig1Ux19AS5aMpeywuhbwrkaQZCywjw8bhdzSwro7NMgUKlnR43gQqDRGHPYGDMMPAxcF3+AMabDGLMF0CESKu0YY+gPhCj1ekabhhzrIwhR5o2GTanX41intFLTYUcQNADNcT+3xG5LlAF+JyLbROT2iQ4SkdtFZKuIbO3s7JxhUZU6XSAUIRwxlHg9o01DDtcIAEoLPI6va6RUIuwIAhnntunMkrnUGLMOuBr4uIhcNt5Bxpj7jTHrjTHrq6urZ1JOpcZlXYxLCjwU57txifN9BAAlXg/9WiNQacCOIGgBFsT9PB84nuiDjTHHY187gEeJNjUpNWsGYkFQnO9BRCgrzHN01JDVD1GiNQKVJuwIgi3AchE5Q0TygZuBjYk8UESKRaTU+h54O7DbhjIplTDrYlxcEL1Al3nzHGu774uvERQ4dx6lpiPpUUPGmJCI3AU8A7iBB4wxe0Tkztj994lIHbAVKAMiIvIpYBVQBTwqIlZZHjLG/DbZMik1HVbzTIkVBIUeRzqL/cEwgVBktI/AqzUClR6SDgIAY8xTwFNjbrsv7vs2ok1GY/mANXaUQamZGhi2agRuIFojcKJpyPr0b40aspqGjDHEPgwplRI6s1jlvP5AGIh+QodYEDjQWWyFi1UjKPF6CEcM/qAzC9wplSgNApXzBsb2ERR6HKkRWM1NVh+Bdb6+gE6vUamlQaBy3mlB4M1zpI/AZzUNxUYNlcbOp0NIVappEKic1x83fBSiTTcDw2FCNu9JMLZGYHVOa4exSjUNApXzBgIhCvPcuF3RDlurr8DuoZ3j9RGA1ghU6mkQqJzXHwiNNAsBji1FbXVAj60R9GmNQKWYBoHKef2BMCWxoaOAY0tR+/xB8tyCNy/6trNqHgMaBCrFNAhUzhsIhEaaaWB0nL/9NYLorGJrzoD2Eah0oUGgcl5/IDTSUQzxNQKbg8AfGnluGO0j0GUmVKppEKicNxAIjXw6BxxbijpaIxg9T4HHTZ5btEagUk6DQOW8gdM6i2NNQw70EcTXCCC2zITWCFSKaRConNcfCJ8SBMX5nuieBA71EcQr0YXnVBrQIFA5rz8QPGXUkMsllDowuzjaR3DqOo+6FLVKBxoEKqeFwhH8wQglBad+Uo+uN2Rz09A4NYLSAg99Dm2Co1SiNAhUThsYjq48WhxXIwD71xsauxeBpbjAzWCsDEqligaBymkDgVM3pbHYvSfB2L0ILEUFHp1QplJOg0DltLErj1qiu5TZd4Eeu86QpSRfO4tV6mkQqJzWP0s1grErj1qKtUag0oAGgcppYzeut5QV2ttHYHU8l3rHjhpyMxgME4kY286l1HRpEKicNto0dHpnsZ17EozUCE7rLPZgDAwFtcNYpY4GgcppI/sVjzN8FOxbB2ikj2BM01BRga5AqlJPg0DltMlqBGDf7OKRvQhOm1AWPa92GKtU0iBQOW2iPgK7dynr8wfxuITCvFMDx1r1VOcSqFTSIFA5bSAQwuMSCjynvhXsXoraWnDO2ovAonsSqHSgQaBymrXy6NgLtBNNQ2Mnk4H2Eaj0oEGgclrfmL0ILFZbvl2TysZbghq0j0ClBw0CldPGbkpjsXtzmvEWnIPRvomBgPYRqNTRIFA5bSAQPm3EEESXfhCxs4/g9CWoYTQIBoe1RqBSR4NA5bT+MbuTWVwuobIonxMDw7acp6s/QEVR/mm3W6OGtGlIpZIGgcppEzUNAdSVeWnr9Sd9Dn8wzMnBIPPKvafd53YJ3jyXdharlNIgUDlt7H7F8eZVeDneM5T0OVpjYVJfXjju/SUFnpEZzkqlggaBymn9k9QI6ssLRy7iyWjtHYo93+k1Aoj2E2gfgUolDQKVdoLhCDd9fxM/39zk6HmMMbE+gtM7iwHqyr30DgWTvki39sRqBBXj1wiK82dnKeo7frqV/3rhoOPnUZlHg0Clnce2H2PzkW427jju6Hn8wQgRw2n7FVvmVUQ/wSdbK2jzRR9fVzZRjcDteGdxz+Awz+xp579+30hnX8DRc6nMo0Gg0kooHOHe3zcCsKO5h6BNy0CPZ3RTmglqBGXRT/DWJ/qxAqEwf248wUN/bWJ/mw9jxt9T4HjPEJVFeRTmj3+e6OY0zvYRbDt6EoiG3w//dNjRc6nMM37jqFIp8lLjCV7vGuQdq+t5cmcr+1p9rJ5f4ci5Jtqm0jJaIzi9w3h/m4/bH9xGU/fgyG3XnFvHN29cQ1H+qc/X2uufsKPYOn/88zhhy+snyXMLb11Zw/+8fJR/3LASt0umfqDKCVojUGllX2sfAJ+8YjkQvYA5ZaKVRy21ZeM3De1r9fG33/0L/mCY+963jhc++2Y+e+UKfru7jff+4K8Mh06txUSDYPxmIYhOXnO6j2Db0W7OaSjnipW1DAyHaTnpbPCozKJBoNJKY0c/tWUFrKgtpaGikG1Hux0718AE+xVbvHlu5hbnnxIE/mCYTz28g6ICDxvveiMbzqlnSXUJ/3DFcr59y3nsaO7h7udeO+V5WnuHqK+YOAicbhoKhMK82tLL+kWVLK0pAaL/z0pZNAhUWmns6GN5TSkAaxdWsPuYz7FzDQxPXiMAqK/wntI09O3nD3KgvY9/f/dq6sZ8yn/n6nncfMECvvfHQ+xo7gFgaDhMz2BwiqYhNwPDoQn7GJLV2NHPcCjCmgUVLIsFwUENAhXHliAQkQ0ickBEGkXki+Pcv1JENolIQEQ+N53HqtxhjKGxo3/kYjW/opC2Xr9jF0hr05mJagQQnUvQHGu/b+4e5IcvHeFd6xq4/MyacY//0jtXMbe4gK89uRdjDMd6BmPPM3mNwMl9izt8gVgZCikvzKOmtICD7RoEalTSQSAibuBe4GpgFXCLiKwac1g38AngmzN4rMoRrb1+BobDI0FQV+5lOByh26b1fsaymmMmC4KLlszlUOcArzSd5BvPHECAz1915oTHlxR4+MyVK9jy+kme3t3GI1tbcAlcsHjOhI8pdnhzmo6+aNNWTWkBAMtrS2js1CBQo+yoEVwINBpjDhtjhoGHgeviDzDGdBhjtgBjl3Kc8rEqd1jNFctjQWB9irbG4dttov2K4918wQLKC/P4+M9e4fFXj3Pnm5dO2swDcNP6+ayqL+MLv9rJz14+yjtXz2PBnKIJj7eGrzrVT2DVCKqtIKgp5VBHv2M1LZV57AiCBqA57ueW2G22PlZEbheRrSKytbOzc0YFVenN6sC0agTWqB07Fn4bz8ioofyJawTFBR5uu2Qxrb1+Npxdxydio5km43G7+MFt6ynMczMwHOajly+d9Hjr/E6NHOroC1BemIc3tl/y0poS+gMhxwJWZR475hGMNxg50Y8aCT/WGHM/cD/A+vXr9aNMFmrs6GNOcT5zS6KfXK1P3nas9zOegUCIonw3rinG09/55iU0VHi5bm1DwmPvGyoKeeSOi9nf5uOs+rJJjy12eLvKjj7/SLMQjNa4Drb3T1m7UbnBjhpBC7Ag7uf5QKJrAyTzWJVlGjv6WVZdMvJzdWkBbpc4ViMYGJ545dF4Rfke3nPBwpFP1IlaXFXMhnPqpzxuJAgcWniuoy9ATdloEOjIITWWHUGwBVguImeISD5wM7BxFh6rsogxhoMd/SyrHQ0Ct0uoKS1wrAmjPxCetKN4tozuW+xcH0FN6eiopbnF+VQW5elcAjUi6XeBMSYkIncBzwBu4AFjzB4RuTN2/30iUgdsBcqAiIh8ClhljPGN99hky6QyT9fAMD2DwVNqBBDtJ3Csj8AfTIsgcLJpyBhDZ1/glKYhEWFZTQmNHX22n09lJlveBcaYp4Cnxtx2X9z3bUSbfRJ6rMo91rj25bWnBkF9uZfX2p25YE20X/FsK3Kws7h3KMhwODIyYsiyrKaUp3e3YoxBRNccynU6s1ilBWtcu9V+bakr99Luc2bZ5Mk2pZlNxfnODR/tiC05XTNmCexlNSX0DAbpcmiOhsosGgQqLTS291FS4Dltzf66Mi/9gRB9/rFTUJKXaGex0zxuV3TfYgc6i605BDVjagTLdc0hFUeDQKWFxs5+ltaUnNZMYY12cWIzlcn2K55t0X2LHQiCMbOKLTpySMXTIFBp4WB7/8in1HjVJdEaghNBkC5NQ2CtQOpEEIzfNFRf7qU4302jQ/0vKrNoEKiU6x0K0tEXOK1/AEaXRTjRb29bdigcwR+MpE0QFOU7sxR1Z1+Aonz3aa9zZOSQrjmk0CBQaeBg7FPpeDWCqpJ8ADr77B1Cal1006dpyO1YjWBss5BlWU0pB9p0zSGlQaDSgLV2/7kN5afdV1mUj9sldPbb2zTUPzz5fsWzrbjA41Bnsf+UyWTxVs8v50R/gOMOzdNQmUODQKXcjuYeGioKT2vHBnC5hKqSfNv7CKbar3i2FTvUWdzZF6C6bPwawdoFFQDsaOqx/bwqs2gQqJTb3tQzclEaT3Vpge19BFPtVzzbSvI9DDo0j2CipqGz6svI97jY0ezcvtAqM2gQqJTq7AtwrGdo0iCoKimwvUbQn8DuZLPJiRrB4HCI/kDotFnFlnyPi3PmlbFdawQ5T4NApZTVP7B2YcWEx1Q7EATWNpXlhXm2Pu9MlRVGgyAcsa/jdnQy2cTbZK5dUMmuY70EwxHbzqsyjwaBSqntTSfxuIRz5p3eUWyJNg0FiNh4kfTFZiqXedMkCGLlsGoqdhiZQzBBjQDgvIUVBEIRDrTpfIJcpkGgUmp7Uw8r60spzJ949E5VSQGhiKF3yL5lJnyx5yorTI+mobJYzcRn41IaI7OKJ+gshtEO41eatJ8gl2kQqJTpHQyy5fVuLl1aNelxVhu3nUNIff4gHpdQOM3NZpxS5o0Gkp1hl0jT0PzKQhbOKeL5fR22nVdlHg0ClTLP7WsnFDFsOKdu0uNGgsDGfgLfUIiywry0WYJ5pEZgZxD0BchzC5VFEzd/iQhXn1PHXw6dsDWEVGbRIFAp8/TuNurLvayZXzHpcY4EgT848ik8HVh9BHY3DVWXFEwZdhvOqSMYNjy/r922c6vMokGgUqI/EOLFg51cdXbdlJvHV5VY6w3ZWSMIjnwKTwdWX4VvyL7O4uhksombhSxr5ldQX+7lqV1ttp1bZRYNApUSL+zvYDgU4Zpzp97cvczrId/jsrlGEEqbEUPgUGexb+LJZPFcLuGqs+t48WCnI/s+qPSnQaBmnTGGH/35CPMrCzl/UeWUx4uI7XMJojWC9GkaKsn3IGJ3H4E/oSAAeNe6BoZDER7e3Gzb+VXm0CBQs27zkW62N/Vw+2VLcE/RLGSpLi2wfdRQOtUIXC6htMCDz6Z5BMOhCCcHg5OOGIq3en4FFy+Zy3+/dIThkE4uyzUaBGrW3ffHQ8wpzufG8xck/Bi7l5mwRg2lk7LCPNtqBFZoTjaHYKw7L19Km8/PYzuO2VIGlTk0CNSs2tXSy+8PdPKBSxZPOolsLGt2sR2GQxGGguG0GjUE0ZFDdvURdPjG36JyMpctr2JVfRnf+8MhrRXkGA0CNWsiEcO/bNxNVUk+t12yeFqPrS4toGtgmJANa+JYHaLpVyPw2DZqaHR5icSahiDaF/O5q1Zw5MQAP/7LEVvKoTKDBoGaNb/efoztTT18YcPKaS/2Vl1agDHQPZj8ctRWO3w69RFAdAE822oEfdNvGgJ468parlhZwz3PHRypVajsp0GgZkX3wDBff3of5y2s4G/XzZ/246tHtqxMvnko3dYZspR5bewj8PkRgbnF+dN+7L9cu4pg2PDVJ/baUhaV/jQIlOOMMXzpsV30DgX52vXnTjmBbDzW7OIOG4LAWkqhNM1qBGWFebYt89DRF2BucQEe9/Tf4ovmFvOJK5bxxM5WHn/1uC3lUelNg0A57jc7jvPUrjY+feUKVs0rm9FzWG3ddjRXnIw1L1UWTf/TspMqi/IYGA4TCCW/U1mbz09d+fSaheLd+ealrF1QwZce202b7mmc9TQIlKOOdg3wz7/ZzfmLKrnjsqUzfp7aMi8i0GrDRenkQDQI5syg2cRJlbHy9AwmXyto6/VTV1Y448d73C7+86Y1BEJhPv2LHbZumKPSjwaBcow/GOaj//MKLhHuuXltwpPHxpPvcVFVUmDLp9PugWFE0md3MsucWA2leyD5DvHWXj/zKhIfMTSeJdUlfPW6c9h0uIu7n3st6TKp9KVBoBzzlcf3sLfVx7fes4b5lUVJP199udeWGkH34DAVhXlJBZMTrBpBskEwOByidyhIXXlyQQBw0/oF3LR+Pt95oZHfH9A9C7KVBoFyxIObXufnm5v52OVLeevKWlues67Ma0uN4ORAMO2ahWB0hE+yQWCFZb0NQQDw1evO4az6Mj7x8+00dvTb8pwqvWgQKNu9dPAEX3l8L1esrOGzbz/TtueN1giGkn6eroFAWgaBVSM4meRcCSssk+kjiOfNc/ODW88n3+3iwz/ZQo8NczlUetEgULZ6rb2Pj/1sG8uqS7jnlvNsbX6pKy/E5w8xEEhu9u3JgWDajRgCqIj1WdhVI0i2jyDe/Moivv/+8znWM8Sd/7PNlpFNKn1oECjbtPYOcdsDmynIc/PD29ZTUmDvhC2rqaMtySGk3YPDzC1JvyDwuF1UFOUlHQRtsVpTbQKb0kzH+sVz+Ma7V/Py4W4+88irRHQkUdbQIFC26Bkc5gMPbKHPH+LHH7yABXOS7xwey+r8TKafwBjDyYHhtKwRQHTkkB01gjnF+XjzEl/UL1E3nDefL169kid3tvLVJ/ZijIZBNkivOfYqI/n8QW57YDNHTgzwow9ewNnzyh05T0NFtM372MmZ9xP4/CFCEZOWfQQQ7SdIto/gWM+QbR3F47njsiV0+AI88OcjFOW7+fxVZ065L7JKbxoEKikDgRAf/NEW9hz3cd/7zufSZVWOnau+3IvbJTR1D874OazJZOlaI6gsyudYT3Id4k3dg5xZW2pTiU4nIvzzO89iKBjmu384hDfPzSeuWO7Y+ZTzNAjUjPUOBfnQj7ewo7mH/7rlPN62yp5hohPxuF00VBQmFQRd1qziNOwjgOgQ0l3Hemb8+EjE0NI9xJUO/y5EhK9dfw7DoQj/+exrhCKGT79tudYMMpQGgZqRrv4Atz6wmdfa+/jOLedxdQKb0Nth4ZyipILA2tymqnjm6/A4aW5JPl39w0QiZkaL87X3+RkOR1joQB/NWC6X8I13r8btgm8/f5CBQIgvveMsDYMMZEtnsYhsEJEDItIoIl8c534RkW/H7t8pIuvi7ntdRHaJyA4R2WpHeZSzmrsHuen7m2js6Of+W9dzzSyFAMCCJIOgPTbiqDaJBdmcVFfuJRQxIzWX6TraFf2/mY0gAHC7hK+/azUfuGQx//3SET7/y50Ebdg8SM2upGsEIuIG7gWuBFqALSKy0RgTv5j51cDy2L83AN+LfbW8xRhzItmyKOftbOnhQz/eQjBsePBDF/KGJXNn9fyL5hbRPTBMnz84o2Wk231+3C5hbprWCKxVVtt9/pGlt6fDCsnZCgKI1gz+77WrqCjK4+7nDtLW6+e771uXdhv/qInZUSO4EGg0xhw2xgwDDwPXjTnmOuBBE/UyUCEis/cxUtni6V2tvOf7L+PNc/Orj14y6yEAoxe45u6Zdai29QaoKS1Iu3WGLNYQ2fYZzpVo7h7EJTCvwp5ZxYkSET71thX8+7tX8/LhLm783iaOdg3MahnUzNkRBA1Ac9zPLbHbEj3GAL8TkW0icvtEJxGR20Vkq4hs7ezstKHYKlHhiOHrT+/noz97hZX1pTz6sUtZVlOSkrJYQTDTi0xHn58amyda2ak2trVku29mG/C83jXIvIpC8mawIY0dbly/gJ986ELafH6u/c5L/EEXqssIdvy1jPfRauwsk8mOudQYs45o89HHReSy8U5ijLnfGLPeGLO+urp65qVV09Lh83PbA5u574+H+Ls3LOTh2y+aUZOFXZZWlyACB2e4+Fm7z09tCss/laqSAkRmXiM42N7H8hSFtOXSZVU8ftcbaags4oM/3sLdz71GSPsN0podQdACLIj7eT4wdn+7CY8xxlhfO4BHiTY1qTTw/L52NtzzJ7Ye7eYbf7uar91wLgUe+2erTkdhvpv5lYUzDoK2Xr8tyzM7Jc8d3XdhJkEQCkc4fGKA5Q7OIUjUwrlF/Pqjl3DD2gbufu4gN9//Ms1JdPIrZ9kRBFuA5SJyhojkAzcDG8ccsxG4NTZ66CKg1xjTKiLFIlIKICLFwNuB3TaUSSXB5w/yfx7dxd//ZCu1ZV6e+Ic3ctMFC6Z+4CxZUVPKwfa+aT9uaDiMzx+yfQ0eu9WWzSwImroHGQ5FUl4jsBTmu/nP96zl7vesZX9bH9fc8yce2dKsy1KkoaRHDRljQiJyF/AM4AYeMMbsEZE7Y/ffBzwFXAM0AoPAB2MPrwUejY079gAPGWN+m2yZ1Mw9u7edf35sNx19fj7ypjP43FVnprwWMNay2hL+dPAEoXBkWpuzWxfXmjRuGgKoLfVyfAbrKVm1pHSoEcS7/rwGzl9UyWcfeZV//NVOfvPqMf7thnNZNLc41UVTMbZMKDPGPEX0Yh9/231x3xvg4+M87jCwxo4yqOQcOTHAvz21j2f3trOyrpT73n8+axdUpLpY41peU8pwOMLR7kGWVif+6bcltkZRQ+XsjqiZrobKQja/3o0xZlqTs6xaUqo68iezYE4RD99+EQ9tbuLrT+/nqrtf5I7LlnLHm5dQlK/zWlNNfwM5rncoyL2/b+RHfz5CvtvFP244kw+/cQn5nvRdmHZlXfQT797jvmkFgTXGPt0/iS6cU0SfP7rdZMU01kTa2+pjwZxC25f/tovLJbzvokVccVYN//rkPu55/iC/2NLMF64+k+vWNMxoJrWyR/q+25WjeoeCfOvZ13jj/3uBH/zpMDec18DvP3c5H7t8WVqHAMCZdaXke1zsbOmZ1uOOdg+Q5xbq0ryPYHSI7PQ6V19t7mX1/AoHSmSv+vJC7n3vOv73zoupKSvg0794lQ33vMgTO4/rHgcpkp4fHZRj2n1+frrpKD/Z9Dp9/hBXnV3LJ69Ywap5ZakuWsLy3C5W1ZfxakvvtB7X1DXIgsqitJ1MZrFqLEe7B1mTYPNcV3+AYz1D3HbJIgdLZq8LFs/hsY9dypO7Wrnn+YPc9dB2VtQe5CNvWsK1a+Y5sp+CGp8GQQ4wxvBKUw8PbnqdJ3e2EjaGt6+q5RNXLHds7wCnrZlfzi+3tRCOmIQv7E3dgyycO3tLL8zU6OzpxGsEO49FQzETagTxXC7h2jXzuObcep7c1cq9LzTy+V/u5OtP7+fvLlrELRcuoL48vft0soEGQRZr6/Xz6+0t/HJbC4c7Bygp8HDrxYv5wCWLM+KCOJm1Cyv4yaaj7G/zJRRmxhiaugZZv6hyFkqXnMJ8N9WlBdOaPb3t9ZO4XcI5DZkZ7G6X8Ddr5nHt6nr+cqiLB146wrefP8h3XjjIG5dV8e7z53PV2XVaS3CIBkGWaeoa5Jk9bTyzp41tTScxBi5YXMkdly3hHavnpW1H4nRdFFvnaNOhroSC4ET/MH2BUNp3FFvOmFvM4c7Eg2DT4S7ObSjP+N+viHDpsiouXVbF0a4BfrWthV+9coxPPryD4nw3l6+s4aqz63jLmdUzWnRQjS+z/2oUff4gfz3czUuNJ/jLoRO81h4dS35WfRmfvGI5169tYHFVZlz8pqO+vJDFc4t4+XAXH37TkimP39/mA2BlfXqNsZ/ImXWlPLb9WEJDSAeHQ7za3MNHLpv6/yGTLJpbzGfefiafetsKXj7cxeM7W3l2bztP7mwl3+3i/EWVXLpsLpcuq+LchvJpzSlRp9IgyCChcIQjJwbY0dzD9uYedjT1cKC9j3DE4M1zccHiOdy0fgFvX1WX8U0/ibhkWRW/2X4MfzA8ZZPB/tboGPuVdZnRKb6yvpS+l0Mc6xlifuXkv8uXD3cRihguTsFqsLPB5RIuWVbFJcuq+Nfrz+GVppM8u7edlw6e4Ju/e41v/u41Sgo8rJ5fztoFFaxdUMGaBRXUlBboJjkJ0iBIQ71DQZq7B2nqHuT1rgEOtvezv62PQx39DMcW7yot8LBmQQUfu3wplyytYt2iirSbAey0DWfX8dBfm/jja51cdXbdpMfua/NRW1aQtpvWj2UF1v7WvimD4Nm97ZQUeHjDkjmzUbSUcruECxbP4YLF0dfa1R9g0+EuXj7cxY7mHu5/8TCh2BDUiqI8VtSWcmZtKctqSlg4t4iFc4qYX1mYc++VqWgQOCwQCtPnD8X+Benzh/ANxb76g5wcHKbdF6Dd56fDF6DN56d3KHjKc8wr97KirpTLllexoraUNQvKWVJVkvMTcC5eOpfKojye2tU6dRC09nFmhtQGINo0BNFJYpPtBR2OGJ7d28GbV1Tn5MVtbkkB71w9j3eungeAPxhm97Fe9hz3caC9j9fa+nhsxzH6/KGRx4hElxmpK/NSXeqltqyAmlIvc4rzqCjKp6Ioj4rCfEq9Hgrz3Xjz3BTmuclzS9bWMHIqCO774yGe3tUKIgjRP4jo1/ifozcK4BKJ3ha73fobiBhDMGQIRiKEwoZgOEIoEvsa9/NAIEQgNPnyu26XUFNaQE1pAQvnFnHBGZUsnFPEwjnF0a9zizK+A9ApeW4XG86pY+OO45M2D/X5gxxo83HlquWzXMKZKynwsLymhG1HT0563EuNJzjRH+Cdq3WfJwBvnpv1i+ewfvFo7cgYQ2d/gObuQY52RWvaLSeH6OgL0HJykG1Huzk5GJzkWaNcAoV5bgrz3RR43LGQcOFxuXBJ9DoS/9UlEncNiV5jIsZgDBgMkUjs51gZIybu/tjPpxwX+/lrN5w7UiOyS05dYYry3VQU5Y/8xwMj/7nR//y474l+2or/RZnY8R6X4HELJXke8twuPC6JfnVHv+a5BY/LRVGBmzJvHqVeT/RfgfV99GtZYR6lBZ6c/2SfjHecO4+fb27md3vb+Zs188Y95pWmHiIGLrT5zeO0C86Yw+M7jk86V+KRLc2UF+bx1rNqZrl0mUNEqCn1UlPq5fxF4/8NDIci9AwN0zsY5ORgkJ7BYfr8IfyhMEPDYfzBMP5ghKFgmKGg9XP0vvDIhTt67YjELuLhSOS0i7vLJSMfMl2xT6IuAZfLdUpouE75MGoFTPTnQgeG0OZUENx68WJuvXhxqouhbHTx0rmcUVXMD148zLWr68etum8+0oXbJZy3sGL2C5iECxfP4aG/NrGv1Tfu/ICjXQM8vbuV2y9bmpPNQnbK97hGwiIX6XgrldHcLuEjb1rCrmO9bDrUNe4xL+zvZN3CCoozrIntkmVzEYEX9o+/3eMP/3QEj8vFBy9dPLsFU1lHg0BlvHeta6CmtIBvPHPgtEXLmrsH2dfqm7IzOR3VlHo5b0EFv93ddtp9R04M8IstzbxrXUPab7Sj0p8Ggcp43jw3X9iwkh3NPfxyW8sp9/3v1mZEYMM5mRcEANeumcfeVh/bm0Y7jY0xfOXxPeR7XHzmyhUpLJ3KFhoEKiu8a10D6xdV8vXf7qc3NgJkIBDioc1NXL6iesqx+OnqxvULKPV6uP/FwyO3Pb+vgz8c6ORTb1tOjdYGlA00CFRWEBG+ct3Z9AwO863nXgOiw4VP9A/zD1dkzrDRsUoKPLz/okX8dk8bB9r68AfDfPWJvSyvKeG2SxanungqS2gQqKxx9rxy3vuGhTy46XU2Heri/hcPc+2aeaxbmP4rjk7mw29aQpk3jy9v3MPDm5to6h7ky39zNnm6to6yif4lqazyscuXETFwyw9eRgT+6eqVqS5S0uYU5/O5q85k0+Euvvz4XlbWlXLpsqpUF0tlEQ0ClVXmVRTSUBHdyOT6tQ3Mq8iOTU3ee+FCKouiyy6/UUNA2UyDQGWdwvzo5KqzM2j7zam4XUJVSQFATqwsq2aXBoHKOv5gGCDrRtSUFUZrBMX5mTUxTqU/DQKVdWLLSFGaYTOJp+LNi75dPW5dm0rZS4NAZR3rQpnnya4/7+hyZUy4AJ1SM5Vd7xSl4niy7IJprafncenbVtlL/6JU1gnH1hvKtnH21sqq2RZwKvWy652iFIwsPJd1QRD7qk1Dym7Z9U5RCgjHeouzrVPVahqKbpGklH00CFTWCcd2B83LsrZ0K9aM5oCyWXa9U5Ri9JNztu0zHo4FgDYNKbtpEKisk63XyXAkWtXRUUPKbvoXpbKOK9uqAjGhWJVAawTKbhoEKutkaxBYw2KzrRNcpZ4Ggco6I6NrsqxTNRQLgmwNOpU6GgRKZQgr1zQHlN00CFTWyrYLZpa9HJVGNAiUUirHaRAopVSOsyUIRGSDiBwQkUYR+eI494uIfDt2/04RWZfoY5VSSjkr6SAQETdwL3A1sAq4RURWjTnsamB57N/twPem8VilpiXbRgsp5TQ7agQXAo3GmMPGmGHgYeC6McdcBzxool4GKkSkPsHHKqXiaNApu9kRBA1Ac9zPLbHbEjkmkccCICK3i8hWEdna2dmZdKGVyjTZNgpKpQ87gmC8P8+xn1kmOiaRx0ZvNOZ+Y8x6Y8z66urqaRZRqcynNQHlFDt2924BFsT9PB84nuAx+Qk8VimllIPsqBFsAZaLyBkikg/cDGwcc8xG4NbY6KGLgF5jTGuCj1VKoU1DyjlJ1wiMMSERuQt4BnADDxhj9ojInbH77wOeAq4BGoFB4IOTPTbZMimllEqcHU1DGGOeInqxj7/tvrjvDfDxRB+rlFJq9ujMYqWUynEaBEopleM0CFTWyt7OVR1HquylQaBUhsjaXFMpp0GgVIbQeoByigaBUkrlOA0ClXVMlq7FoE1DyikaBCprSfb2FitlKw0ClbWytWaglN00CJRSKsdpEKispU1DSiVGg0CpDKMtXspuGgRKKZXjNAiUyjDa4qXspkGgVIbRpiFlNw0ClXWy9Tqpnd/KKRoEKmvpZVOpxGgQqKyVrTUDpeymQaCUUjlOg0BlLW0aUioxGgQq62T7qJosf3kqBTQIlFIqx2kQqKyT7aMss/zlqRTQIFBZR5uGlJoeDQKVtbKtZpBlL0elEQ0CpZTKcRoESimV4zQIlFIqx2kQKJVhsr0zXM0+DQKVdUyWjqvJts5vlT40CJTKEFoTUE7RIFBZR3SgpVLTokGgso42DSk1PRoEKmtpzUCpxGgQqKxTmOcGsu8TtDf2ulxZ9rpU6nlSXQCl7Pbgh97A4zuPU1NakOqi2Oo/blzDT18+yrqFlakuisoyYjJwKML69evN1q1bU10MpZTKKCKyzRizfuzt2jSklFI5ToNAKaVyXFJBICJzRORZETkY+zpu46WIbBCRAyLSKCJfjLv9yyJyTER2xP5dk0x5lFJKTV+yNYIvAs8bY5YDz8d+PoWIuIF7gauBVcAtIrIq7pBvGWPWxv49lWR5lFJKTVOyQXAd8JPY9z8Brh/nmAuBRmPMYWPMMPBw7HFKKaXSQLJBUGuMaQWIfa0Z55gGoDnu55bYbZa7RGSniDwwUdMSgIjcLiJbRWRrZ2dnksVWSillmTIIROQ5Edk9zr9EP9WPN/3FGrP6PWApsBZoBf5joicxxtxvjFlvjFlfXV2d4KmVUkpNZcoJZcaYt010n4i0i0i9MaZVROqBjnEOawEWxP08Hzgee+72uOf6AfBEogVXSillj2RnFm8EbgO+Hvv6m3GO2QIsF5EzgGPAzcB7AawQiR13A7A7kZNu27bthIgcHeeuKuDEtF5B+tPXlP6y7fWAvqZMMd3XtGi8G5OaWSwic4FHgIVAE3CjMaZbROYBPzTGXBM77hrgbsANPGCM+Vrs9p8SbRYywOvAHXHBMJPybB1v1lwm09eU/rLt9YC+pkxh12tKqkZgjOkCrhjn9uPANXE/PwWcNjTUGPP+ZM6vlFIqeTqzWCmlcly2BcH9qS6AA/Q1pb9sez2grylT2PKaMnL1UaWUUvbJthqBUkqpadIgUEqpHJcVQSAiN4rIHhGJiMj6uNsXi8hQ3Oqm96WynIma6PXE7vun2CquB0TkqlSVMRnZtOrsRCvrZjIReV1EdsV+Nxm5A1RsyZoOEdkdd1tCqyWnowlej23vo6wIAqIT0d4FvDjOfYfiVje9c5bLNVPjvp7Yqq03A2cDG4DvxlZ3zUQZv+psAivrZrK3xH43mTru/sdE3yPxplwtOY39mNNfD9j0PsqKIDDG7DPGHEh1Oewyyeu5DnjYGBMwxhwBGomu7qpSQ1fWTVPGmBeB7jE3J7Jaclqa4PXYJiuCYApniMh2EfmjiLwp1YVJ0lQruWaShFadTXPZ9PuIZ4Dficg2Ebk91YWxUSKrJWcaW95HGRMEM1wFtRVYaIw5D/gM8JCIlM1OiSc3w9cz2UquaWWK15fwqrNpLmN+H9N0qTFmHdEmr4+LyGWpLpAal23vo2QXnZs1k62COsljAkAg9v02ETkErABS3gE2k9fDJCu5pptEX1+GrzqbMb+P6YgtEYMxpkNEHiXaBDZe/1umSWS15Ixh5+rNGVMjmAkRqbY6U0VkCbAcOJzaUiVlI3CziBTEVnNdDmxOcZmmLfYmtCS86mwaGllZV0TyiXbkb0xxmZIiIsUiUmp9D7ydzP39jGWtlgwTr5acMex8H2VMjWAyInID8B2gGnhSRHYYY64CLgO+KiIhIAzcaYxxrMPFLhO9HmPMHhF5BNgLhICPG2PCqSzrDH1DRNYSt+psSkszQ8aYkIjcBTzD6Mq6e1JcrGTVAo+KCESvDw8ZY36b2iJNn4j8HLgcqBKRFuD/El0u/xER+XtiqyWnroTTM8Hrudyu95EuMaGUUjkuq5uGlFJKTU2DQCmlcpwGgVJK5TgNAqWUynEaBEopleM0CJRSKsdpECilVI77/177RZ5C/pP7AAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "h_energy = part.E.Ea\n", "plt.scatter(distances, part.E.E - 2 * h_energy)\n", "plt.axhline(y=0, alpha=0.5, c=\"grey\", ls=\":\")\n", "# plt.ylim(-.2,.1)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }