diff --git a/Wavefunctions/GaussianOrbitals.ipynb b/Wavefunctions/GaussianOrbitals.ipynb new file mode 100644 index 0000000..3ecfc63 --- /dev/null +++ b/Wavefunctions/GaussianOrbitals.ipynb @@ -0,0 +1,815 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "import numpy as np\n", + "import math\n", + "from IPython.display import display\n", + "from sympy import *\n", + "init_printing()\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "\n", + "def display_eq(a,e):\n", + " if isinstance(a, str):\n", + " a = Symbol(a)\n", + " display(Eq(a,e))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Gaussian Type Orbitals\n", + "\n", + "Gaussian-type orbitals are frequently used for single particle orbitals in QMC." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "x,y,z = symbols('x y z')\n", + "alpha = Symbol('alpha', positive=True, real=True)\n", + "r = Symbol('r',real=True,nonnegative=True)\n", + "i,j,k = symbols('i j k',integer=True)\n", + "N = Symbol('N')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The form for the primitive orbital is" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAH8AAAAbBAMAAACn9k1sAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vEM1mIomZ\nMqu7iC+qAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACX0lEQVQ4EaWUT2jTUBzHv02aZU2aNN3BOayu\nDkHYYWbgLiK0IKKIh4Eg7rQpAxEchuHJzbbYk16sgoq3HIZehrSnoRYtO02nuIMO9JRdd6pOEcZk\nvrw/GUs6J+s7vN/39/l9fz/yXpoCbSzjdX8b3aRVwlB7AzpRdKITMgEqc7WvGqCwOB0GJP8WsDxX\nih2gkFDvhUDrtNVzMmehxQlazNjxptT3U1H7oMeZmslzVTcOccVDrQ4n8e7BqrL5a3vBz+Z+c5aM\nl7i6e2r7IRTbtKqJnH2B17s2bcRHr7DjaDeHOb6le1xtLHHBwzwSZft51xa8nwUUcR2SJQpShSnj\n8hGO4n1kHcU45OvAivBBfrgO6CIt3hYq5TIVyy4bgtG4DoOccyJg+swG8FmkZ4NRyxxJ1peAUVSG\n3nS0puiAiVGgW6SFV0KJQ3a65kHBaJztqSzm5WzAkuitoAGo6dXBmk/VuoVnDg4TKV+0v/qIwFrG\npSK6HUdxGFVgCo+sS35Zj4+h1zxxlcjJ9AfPR8CQqxFPyzWH5Jpsk18A3jrnfMdMsoQ+841LpIMX\nPgG0cUwvMRnZG9B+6I5vvsNqToeF79xmeEyY19JPOIqEBnBsllLRlXOMP9x2kseUx0U0yHkgt+hz\n0qXS+gJiY1Qg3mARqZ0ugFwZeXpljfjOxJr8bQ5A8ljjeYCpDjLAYiy8KwSQVsg/paZWodV55Jgw\nRl7utymKlTDtUhXeJm+MEOSRF73SfeApqxZqn1yqkngs/tcyPWwmc+y6L+zq+JfBLKn+nex9KdbW\nN72nKXL64//3/QXCBXzOM63xGwAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$N x^{i} y^{j} z^{k} e^{- \\alpha r^{2}}$$" + ], + "text/plain": [ + " 2\n", + " i j k -α⋅r \n", + "N⋅x ⋅y ⋅z ⋅ℯ " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "gto_sym = N * x**i * y**j * z**k * exp(-alpha *r**2)\n", + "display(gto_sym)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The normalization ($N$) is given by" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAANkAAAA/BAMAAACfjXDaAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHWklEQVRYCbVYaagbVRT+8pLJ+l5eVBBUbEcr\nKqg0bqggNQoKitYRxVq3xA3FNdal7g2CiAu8J25V/0wVXH5In4IbD226WEGqDfpPi0ZxQdyeGy61\nPr9zl0lmJnkvYnogc875zrnnJHfu3LlfgAFk1exwZGaAXoA/UNawkqaGVWiQOil3kKxh5RTbw6o0\nSJ2RyiBZ/yVn6eaukiEHKPQptCHACw1jjj4eYHMZR491LYWQA+T7DDyeeFbFxuTakMvzcplXvEy9\nkxNygN07kZiV7CBlms5vHX8OK3VDVzDkANu6QlEz0i3rRxN6+9mumUTIwS29RzxYI666nX+IZMhv\ny7in1sSeR5r4qZMRcoAjO5Fu6yaZNdXN/Vhw6Vb4/p4B5tKZTE3KCCUhh4gv6FK5BPIFkLjApyvd\nnNLNEpBu4/uJNYfsuugIRrd1PwEhB46USZThbFnRsHXSFTZq4qrpN6an28ASpKenn5p+CydvbKmU\nS5Cv0Rh7VLRzjE/7WX7g1LCqIUZfydYZKngYAa4JkqaA6ln09CpRd0G+1MbUZEVyrsXCGapiTemC\ny687K3jaQ24S43ylSFhJyAFGW0TfAT4B9jIZyv84Q09103dBuh2b8us6p+p39DgrOH+Jn/OR/lUH\n+lzTJQYeAviMnBukLAfOvp6e6jYqjdR9ewlXV8TmHZzq6AmP9nbxC7/pbklBrIQc2UrkOZo4EnfZ\nBBRa2lTdeAcpZY2Ya66mDdG3i/mY9pGXRfuNdMscfKDClGOiapvMuJzz2a/bQYoAIryZuGIP9U1d\n8ZU4K5rIu8Cyjdz2qF/D54/gaR3DRJ1f/hSOSDyP5aW9jWOCwAm0km1eFm/3gpTiVBDH5Rd1bGVl\nspNIsqiba4t2/lh2zmocZZKeob4wz265OpLu28YJKvxAK9/iAttv3eogZcwP4nFj75EppNsolrg8\nqMd+bTLnWp03UqZuSrfFXOZTdePoIK/ywOYawFco/unZlNRMEI8b3ngTmRYcTJREj2znaNyp8zZT\npVbeeiPwMG/L455xdJDXKX6k26PA6S2Tgjm74VMPoxUOu43PTwX59VfR3sQPn3cXu1DlPDi/s9sO\nweh0xKfJmXT4a9Jf2hRZpP3lVQcpKfEyf4eHaqXaKuFElc5FsJ4P36p9gH/YjTW1Y2slarRklfCR\nS5ZsSvcqsZmBdnbIg8866qvjR1QbdR1MPLdyf1+b2B/4fEZ/B4NQjTV5ydSABzx8F6SMuET7SXZm\nAZZPcSgfLepLMF6pO09JdoF7lm+GZe76AQc0jWNVpkFLJi51MXdlmxKcRGxaSB/bVLvyCEdxV16C\n4sFmVw5l9XL0vrImEtoj4vd05cn+jyIbF2c4MuqOiN/DTXkLKz3guaEzVThZCmU5vCvzycL6YfOl\nxOP3KihRDkUizUMx66Q3tK05kE5I1ns6dUFoxGchbziOanTRcGrNX0Uds/3584aSUfyFZQZZD0Pp\nduELLJNyh1LLFOnLqIH2dcwptnmRc9H/F9Oxn9pa4bbPz86SMHP7tM19dGe1Yt0wc1voclfdid3C\nzG2cG8hJO7FbmLkVJsFXmpFU21ooVgJzaWBxQfWG8UVXTtjMdm+3ab4N9QGTSbvxBLjv4SyYLuNS\n+s/KQG6fGhSCEYc135CFxtdcXMLMbfQP4BWbRONsZIik27J0NIsgGdEgiiwXhwuuDOcvUCRE7C6J\nMLcsTyG+CY/WgcuAJ5RLR7OId7pAnq9isPANRVGqvhoYvkSY288yVVpyJeBJ4LqKuM6kYRE8GQVg\nL3jCk3RSFEtCxO0n2zFWN7H7qdd5phve1Cwi63eDiMOKbwhFydVMnS4VYW5YJ0ddLSu0Wust44kG\n9J6myrgKXespgtEDVnxD0nofUfSZR9XgZbEXvKSXKCz7O+r30foWikUk24ISVAQjDmu+AR4+k3aO\nZEAg36iJtu6PjWAreVFh+XKidRCtjzSLEDLC1V3WBCMOa74BUhQu5C7Zyi3+5yiXwkTzPJuzWhmL\nSCX+pnW6ZhG5hqCLuD6EYMRhzTeENAR3RAZkbty2YD3zFbESQEu1drc1VbfRmj4AS7dNjKhuAiqC\nobqFYM03pJsiIZInchquV74iVhqSa85/3zpqJj9AolSocbK5Qk9kRM2kgIpgxGHNN+SkrUiIrQWU\n+edClEshOXOpzZBVwvf4aKnaqqvlIAFZJQo0BIOLRySADd+IUZRUTa+HHHt2ZOQn3zryd8+eK1dc\niYlEG9gXikVkagZUBCMOK75Bcu4qEmJLEWiPSx9FrDpocYdvHc4d1s3O/oXkFlqv612ZU6RBCMGI\nw5pvoNCI7MpVb7wh+WFJ/VKzQL5kLerEpHHWWDDv9oR1OEZR3hX+HpfZpsWKrrWoM03jGDJiCEYU\ntiNiFIX/i/XqtrZhR+DQwAKCp9BsNYZgRGEzYuAj6XFe0INv00CCzuYVkd7QllgUNvnmKxlvDnVr\nJyZL0UixZC0sCCyeFHrD+KwrZ07zwzmjww6eMdyC/wKu0Imh24CfHwAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\frac{2^{\\frac{3}{4}} \\alpha^{\\frac{3}{4}}}{\\pi^{\\frac{3}{4}}} \\sqrt{\\frac{\\left(8 \\alpha\\right)^{i + j + k} i! j! k!}{\\left(2 i\\right)! \\left(2 j\\right)! \\left(2 k\\right)!}}$$" + ], + "text/plain": [ + " _________________________\n", + " ╱ i + j + k \n", + " 3/4 3/4 ╱ (8⋅α) ⋅i!⋅j!⋅k! \n", + "2 ⋅α ⋅ ╱ ─────────────────────── \n", + " ╲╱ (2⋅i)!⋅(2⋅j)!⋅(2⋅k)! \n", + "────────────────────────────────────────\n", + " 3/4 \n", + " π " + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n1 = factorial(i)*factorial(j)*factorial(k)\n", + "n2 = factorial(2*i)*factorial(2*j)*factorial(2*k)\n", + "norm_sym = (2*alpha/pi)**(3/S(4)) * sqrt((8*alpha)**(i+j+k)*n1/n2)\n", + "norm_sym" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Sometimes the normalization is written in terms of a double factorial. The expression can be converted between notations using identities like $\\frac{(2i)!}{!i!2^i} = (2i-1)!!$.\n", + "\n", + "For more information, see the [Double Factorial entry at MathWorld](http://mathworld.wolfram.com/DoubleFactorial.html)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "QMCPACK splits the normalization between the radial and angular parts. The normalization can be derived." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Check normalization\n", + "\n", + "Derive the normalization for the simplest case." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHIAAAAbBAMAAABSCMbcAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIokQdkQymVRmzbur\n3e+SS/cOAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABg0lEQVQ4EWNgIAskaQuQpY+BYQPrBTJ1LmD+\nTKZOBg5y7WSwItefzIfJday5sAJpWpmcCpgZps/0Ebj/H+xPRqK1X2U4wsQw9StvAFQHE7E62Tcw\nzAplNNkAV0+Ezu7dQGDA18AQV8DA4wDXaQJnEWDwCzDkGzDwJcCVhcFZBBhcAQz6zAzyCFUKCCZ+\nFmPJTHYNhmiYIraEwkYYGytt5FKIVZyhgnEClwCLigIDDIOUaa0CgeUgJutFBiEQjQG4J3AKMB2Y\nxNDDAMNoSuQamMvQhCBcXwYmBvYL1xj4DWAYTdm5KZ4GaEIQ7nWgTq7NHxn4mqG4AU3ZHzQ+nPuT\nQZJBvu0zA+80KJ4AktK9CwJ3QMxPIAIbWMIwkWE7G1CXGxQ7oKl6xMDAgiYE4cYyKHAXsAFd6gbF\n6Dq9GJixxwpb5U41BuZvDPzBUAzLBzB7WEpmwphoNCPIY1eBYQvDaPK4uZwCQDl3YHzCMG6laDLg\nPMZUVMwAw2jyuLnYAw63epgMAJXSX3BU4RMrAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$\\phi_{000} = e^{- \\alpha r^{2}}$$" + ], + "text/plain": [ + " 2\n", + " -α⋅r \n", + "φ₀₀₀ = ℯ " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Derive the normalization for the simplest case (i=j=k=0)\n", + "gto000 = gto_sym.subs({N:1,i:0,j:0,k:0})\n", + "display_eq('phi_000',gto000)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAADoAAAAzBAMAAAA0mgbTAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAInarRM2ZVBDdiWbv\nuzJCz3LGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACKklEQVQ4Ea1UPUwUQRT+doFzbnWXtSIWBrAw\nJhhjQQFo9BoT7SCxsDBh9VQSfwINJFZsoMJgpLCwc2NBcSbmQmEhxVWGSCCRHhISCg2xoJB4YMz5\n3u7O7MzCXWF8xfv5vjf7vZ2dWSC17gabrPKxPw/otTWsV/m8MGEgiY4SmjNI5HTGYD3a2pYtOR2x\ngFk4PyWb03GqGADWJJvTofILcF0+eky2JbEPGA0VSzrS2kuUJeXbMAFJx+adC4CLtA92D8Pth+zJ\nSKfyrXfpNGDN0KTFGCx+jQPQB2/C8S9xJf4An2L4GXml4/JLkI37eMPRvkxO6eyks3QHJyNmp2D5\nmc4T1BlEW70YUvAC2L7SEXUkO+cdTnHTmcriIIVUx17AKqPA3TvsRxuNAwqpjlvFEKPAdBCHxKU6\nGdLhZzlSHQ0x0kTHgLTC0NHwJDV0jrD/DeCv3NT+RcVbutZimVuabMEKVFqw2JVH7PimG8fDMWqj\nU/9uuc55/NAQPhzOuw8KOXFOe6OzdPzEbTxUrJ68LBPbEWFLgby1qrCJrTUdg9lN1UuJV9nISmb3\nd8+XaLT7Dxh2cctgxX4VK7BeY97voiHwKlQ0rRWNEE/DUxHc4ALjzxUJfjLdg1ppPERhOCJCjJjs\nZ2Kr66S8GhIxyy61dObaTfo1OL8IFPdeSI7uNT15mnTxm9g9wovxhdPWupFYQRn4uPc9W8ZZ4erB\nFWCxvA3n/TIe95hs8+ovPnWwlHVUvZAAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\frac{\\sqrt{2} \\sqrt{\\pi}}{16 \\alpha^{\\frac{3}{2}}}$$" + ], + "text/plain": [ + " √2⋅√π \n", + "───────\n", + " 3/2\n", + "16⋅α " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Integrate the radial part\n", + "val_int = integrate(r*r*gto000*gto000, (r,0,oo))\n", + "display(val_int)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJoAAAA2BAMAAADND52xAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vEM1mIomZ\nMqu7iC+qAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADRUlEQVRYCa2WS2gTURSG/0mamTaZSVLxgaWk\n2UlBaRQFN9p014VgEEREpEUEwSJOxSriI12Ij5WjglJ0EUEXSqERoYWmpaO48IXNwgraSiOuhdIs\nArVSz9yZSWamaRdzPYvc/5x7z5+ZmztfBuAI5e5vjm5vayDV7S1x5BIKHN2rWmeTq0o8hb08zZ5e\nAbGcp8SRXsBXjm6jdfjNjqpDQzPnCZFKmEpV7XhFVIXYa5usUNja1yhmEP1b69xWk35UY8Xp1nTK\nj4erJ1QBJlrZ5u074prxk3Tm8SwlnQ6nIWkON58IOAT0AdsFFQIcbv4QICch072+iwCzhb5i9e78\nISABiElgdMzwcf4KfhDQoONFNg9Mpcns8nKOPu3wgYA5oCWmkZvtYY9+EBA+VvickUtomlZU28cc\n/SCgkR6mDApDmvzSbQYLARtWUmjoOeH5Ks/adVMnV3AnCYQ4zCTGFQsBwXtLAJ0e32FzhSEg8nQZ\n+OLbi46fyRUTAVH0AFs43CyumAiQ0aZB53CjVuKKhYDdyGaQd7pNF9HKcqnD+HEit1joziUeTVyx\nEDAGuRxMOaezwHcz3+Msr6OJKxYCdIQXI64DUnu7mPA6GOfAHQu0hLhCwRCgAx3DRlaNmsftam1d\nQVyxEBBMA50fgcdDiZz0OlHElebr2FkC03QSad8OsygZum4YXLEmjJsMlYH75VBqoDivSseVCsYq\nMLTSW7fbW5wjrli1EI3CAsL9SeATZiDrSm/4fIZpgWprxRNz+2jbGFfMZefOHiVRAqJx4CQSENNC\nCYE40zTUjYAK/KrOmFyppkyIRUhLaL+YzQUeqtmrTIs/3Gvs7IGK8KCd1B1HqHpGKWticf4buiNM\nx4S6S4MzKuR03Sm7uIvE/o2tWvBn13tcG2dafmXPusZLxK8BV4UniZPbWx4DZ69SGJ2UMs4KnxbV\niAa00RlZ5jMyuqWpD7Rtkcmu5y05fjdy2GrgelwocpsxerFtS9I/Am9sHqRHMU4uSsl4Kvki3E+I\nYDayFuO/NmkReGRcUVaNpfiuzOiezplvgZsg5vndRvQm5tKO6H9wi1Z43jO8d6P8OeAt8eQ3bvJ0\ne3sP6t4KT96YW6P7HwYt8fXjSQTPAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$N_{radial} = \\frac{2 \\alpha^{\\frac{3}{4}}}{\\sqrt[4]{\\pi}} 2^{\\frac{3}{4}}$$" + ], + "text/plain": [ + " 3/4 3/4\n", + " 2⋅2 ⋅α \n", + "N_radial = ───────────\n", + " 4 ___ \n", + " ╲╱ π " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Normalization for just the radial part\n", + "radial_norm = 1/sqrt(val_int)\n", + "display_eq('N_radial',radial_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAHsAAAA0BAMAAABY51r1AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAdrtEVN3vEM1mIomZ\nMqu7iC+qAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACiklEQVRIDaWWMWgTYRTH/5c0PZq7M9VFkAaD\nSxfBGwRdxGw6iARExEECggp28HBwEezqlqubiJChXcShXexwBYOLSgU7FJeCZJEOIhaDBGolvu9y\nd9/3JV/P+y4PXvP/v/d+l4+70HdA9nCWfvBhyfByiiq4l3lXMrycokwEvCsZXk5TOzWhKxmhniLP\niz3JiI1DtIFKO2lJJqmmiUf4wtuS4WWlKn86tQVMHRUenGSUEC8+gdXjTlvtAneBAUWCSiapqsUd\n4CedHvNiWzJiY0y3vBCfuSd0JCPU1bLpARdusN7GnMs+hoapDFH8A5g+w1dc8365HpkMZDhSqgEG\nGL4AnDa8yGTF39HgTrCwBbsPfLAik5U2uuEk3brpGvBmnTnxPobdw/9cQbkNPD5oY3EVeFunSWYy\nhtOBEQ1XfMIzYvHYSrDxLdJ2FzOfHXqKGtEaDH7H48Fz316LDX0eG7iYat7Wu6DAt2pAKTddfLYP\n0JPMGdbyAbCdEwaOoAkcz43bOOmjkxs/i8UGVmGs8aRrWU/D+P9l12H3ii5e4GqSOkfpoPzL8nAL\nFS/OEbzC/rGpYo8G6XxnXsPcx/SlKOsjeJot0vDFTTh9lF5G6dO8dT2M7jgqr2Y6N0q9EK8SznJ2\nnBEq8mouUcfYg0MHr0aZjkur+eGDm8R3Yf5F5VyUrvBdCqlczdfozsepYMSSajXP0XOPU5wd0+rV\nbG2+R5xjiFjQWc0iF2md1azAJyyxX++El5gIn89FG+zYHUBvsydfFWyfWHpFTmuzJ7RTt9pVcsM1\nn5SziwJ2aXi45rNTyeRXNEiHaz6paYiPoHcDCo3NLlzd7CN8OdTY7AINo8FebHJHwcf3EfgfeQPK\nAD6Tv5oAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$N_{000} = \\frac{2^{\\frac{3}{4}} \\alpha^{\\frac{3}{4}}}{\\pi^{\\frac{3}{4}}}$$" + ], + "text/plain": [ + " 3/4 3/4\n", + " 2 ⋅α \n", + "N₀₀₀ = ─────────\n", + " 3/4 \n", + " π " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# The full norm includes the factor of 4*pi from the angular integration\n", + "full_norm = 1/sqrt(4*pi*val_int)\n", + "display_eq('N_000',full_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAADIAAAA0BAMAAAA6SHafAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABiklEQVQ4Ed2Tv0oDQRCHf0nucvljSLRXIkK0\nEExhYSVBsBGUa200j3BiINiYvIEBC9uAhZZ5hKRIChtfIIXY2UgUg5XG2ezu3YxcYe3CZeebb+cu\nd/AD9Doa1UxFm4CdQi8yAnwviIwA5zwSEIAMu5uAPt6iIQ5ux+mExsLS2jb1xvxfa3DraLXD47xI\n+8h2UJzNZhPbNpDtIv1BvZRvhYX8VJsXbizkpkDmQBlvc0NNaqDiKgBOcmQS9zguLRtQR27p6iuT\nDZAqDw1QM1kFnOZFA9jykewFBsiM6KLzNHNNT7qhfQ5AoYxFwG2twv0k80VCAa09YKB2Wt9kJrqk\n38Rds9I1VAGeJ/tW5em7WONdvmK9b81/2+lF49efX1TkQ0yJfAgj8iGMzIdQIizc8HzwPmw+RFOD\nCEuMj2mJsMT4+JaI0fzII33wd6pscsI5rzFeGZRYckJziLOFGtE8RmFXF1WKcpQcJp16jkjHiLWp\nTD4V1YxJDnenfrFNbJLDzYMKaezaRfqX+QGRc3/QBcDTygAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\frac{2^{\\frac{3}{4}} \\alpha^{\\frac{3}{4}}}{\\pi^{\\frac{3}{4}}}$$" + ], + "text/plain": [ + " 3/4 3/4\n", + "2 ⋅α \n", + "─────────\n", + " 3/4 \n", + " π " + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Compare with the previous expression for N\n", + "norm_sym.subs({i:0,j:0,k:0})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Higher angular momentum\n", + "Now check the normalization for higher values of $L=i+j+k$." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAI0AAAAbBAMAAACzTxHIAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIokQdkQymVRmzbur\n3e+SS/cOAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACIElEQVQ4Ea1UT2jTYBx9+ZLFpG3WwBDx5IcM\ni4IQpLDLoKUMxnBiTt7KcvHczosOPdSrIFb0bgWHsIvTHYZ4MII3ETxMxIPsMrTeOspQEdHvT2PM\nt6RG5g9ef+/33suPL00IcPAyjz8++BK2oYjmf9lTQM3Lv8h69CQzfDPTSTHupWhS0u9nWinGqxRN\nSrP/cFvATtYefflMlpWi63uqOF2BZ6/cfeH8/CEsTQ2kztZAkR2/VKf2kv820klExnbSU+yrsFf9\nqbOxmm9PoQ3TFxdZz1htoI+JT8DLeE81pmNYuYNSwt6FOQQ+xtqFmI5hLReLCXsVZOBpf/xrNGFn\nDGT78uvkEzuy0L3UmViP8mawfD3i6b26eLhRH1n69LybnrqidQsubnynQKMhwYJkTVSbUZ3aW06P\nEV5NV6O876ti13ZJCHsA2NQKOZQM8cjQGh1C62MmUHw5ngPBoTcorQNOgCGHktNRYKas0ue5OxFP\n9ndsD8tN1oGaiz5HMsCmMjNllcOIqf0bjuJYBzUfaHl4z8Ej5IOoUMRbATxBUKay7/99gNvYBG4x\nh+3Y4lBClncexmjPJAXqii/HJmixDZxiE3vTdjiU3FLnOTuzLKOHGVfx5WiuPD3B2EPgYi3AHoeS\nm52rnuxG2vzCbxpJo64J4ysQOj57XgxKIOdo83OauzBCo20J5LxQiYmvxukv17Z9VCqegJLINxr5\nYn9J/QJmLHbBJIgx0gAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\phi_{L00} = x^{L} e^{- \\alpha r^{2}}$$" + ], + "text/plain": [ + " 2\n", + " L -α⋅r \n", + "φ_L00 = x ⋅ℯ " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# L = i+j+k acts like angular momentum state\n", + "L = Symbol('L', integer=True, positive=True)\n", + "# Use just powers of x, general enough for now\n", + "gtoL00 = gto_sym.subs({N:1,i:L,j:0,k:0})\n", + "display_eq('phi_L00',gtoL00)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKkAAAAVBAMAAADP89MSAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEHarIkSJZt3NVLsy\nme8Q6PJIAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAC1klEQVQ4EcVUS2gTURQ9k04nmWSahBYERW2o\nUHVl1K3Y2YgLpc2qiAodWrUgClkURPATVESEYhQUdKGDuHBnKEJxoR0EUVw0sRQFUWlREQShtZRC\nQfDcl89M6Ih15YG8Offc+868d9/LAP8bVh7QTiHqhS/kpcg6cCY8HVQnM360j/TQugp2+VKA6SUJ\nDKDFC6jhtOz4+nXAmsU8TF8KMDMtAV1RDKh/pREHiOWxDCEr0a2kdo5VtrIiVEnawDlgCdpMWP6L\nEg9z3BiWbtL2e43wANkorGXgTkPzSWtBcYej6QLtm/Z0Z8gD6N/e38Xw+f2t+kTWuvuha4ckTwKt\ni4gsKgZsOHYc2mB3TsqYNSocLG/6gaKaE70Qk/4atxUk6Txxn/ESlRBDIgtzDH021XEe78LO3XPA\naUb6d0zar2yMqTIgkaf4Xs+bNqJFGGljPi6TfFjuRbA/rbM8UbomZpCSKVd5WEWkisBDRmYRj7zL\nQPmFlDHlinPUNrKIF6DBLFEMQtO4S+LamK1caZVh+BVIOSi7dGLUWeGG2eOUI2VAMgdsoX9bAdYc\nYzWFTx9xkYE1934p11LDtdPFJybEdSADxPn2VEnKqq6jdOUalWvZQ5qycUUhS9pW5ADdxvmcdKDm\nyg505vCTGbkNjbVOS1m1A0vYiz4XvA7x9DZExDWApCNB3EFLPuDK00rmjAIzclrsK9wfXPNBKaue\n1g2M4CnktAbcCb6hGX18N11nkbC50Ppah7iJTK/NzEf+9G8w7EkPt1QZt5qVPTiJCqsq6OhqH+ar\nmiA75B0Y3nzU6llY37PweuIzY6ra1BvJ3JShY+oI9LdDOSljKP8C6934CVJTLYpkNYjZtSp9Jqz8\nEkVdra83LP0nLcrNKRiZGml6nGUUtUUalGHVeFyrXBs6I2GzuZLRiqH5P4ny1RbU3atRfZSvdkQC\nsf8HVLtW2+bKefWbNAL8BnADrTi2KbmMAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$x = r \\sin{\\left (\\theta \\right )} \\cos{\\left (\\phi \\right )}$$" + ], + "text/plain": [ + "x = r⋅sin(θ)⋅cos(φ)" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEwAAAAZBAMAAABp+tVcAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARImrEHa7zVTvMt2Z\nImbh7FZmAAAACXBIWXMAAA7EAAAOxAGVKw4bAAABsElEQVQoFY2SP0gbYRjGfzljNHfeJbipg4JI\nBzMIOjnUc2mXokEROigVFEtBQRfJGBykpYJRcLNUkKLtoG5OLVmkS4dM3ZRQEBcHo7bVwV7fL3ff\nmSFK3uF7n3ueH9+/+6BcqVnXF4+ORqFu/lHAD62CeVcDFk/apRowaLyuCbPSVbA98YaKOFthlgqU\nuRFa2JMQmV6Bfu05+UCZl9qCBhc+uwmXBu0106rlfR8R+Y1ElvoO3zQzb6fvY60+gJMjUSSa8624\n5134yv6iISIHYOXZz0KGSOdgV3BK8+XY/M8rdjY/jSrHyEN7kldJSGHEMrGcP8MzOGAKFtO2+iex\nIgz09n4Uecau1RPZEiXV4tKhsDfwTz6tLBzCX5HvcYdlUr+sPxMo7DX8FicuwSXmjch9WA0giD7x\n8gqbCbHoDU1qR4J9D7FjnFIFJotGSyTSkrfBeogtwdMKTB1hjW4Vn1Hxgt65zKkVg0WNHjiy5Y7h\nR/kSlZJqPv2V3Pb6tr3+5VtB1PW+GHdVsqCGh+qrDuycVtX6uTaNtFbVeqwQuM+rpaGnnmW55oL+\nQDvxfafAfxPMYQlZFON7AAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$r^{2} \\sin{\\left (\\theta \\right )}$$" + ], + "text/plain": [ + " 2 \n", + "r ⋅sin(θ)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Will need some values for the spherical integral\n", + "theta = Symbol('theta')\n", + "phi = Symbol('phi')\n", + "x_in_spherical_coords = r*sin(theta)*cos(phi)\n", + "surface_element = r*r*sin(theta)\n", + "display_eq('x',x_in_spherical_coords)\n", + "surface_element" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAAcBAMAAAD2EOVAAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMARImrEHa7zVTvMt2Z\nImbh7FZmAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAEtElEQVRYCcVWXWgcVRT+djLZzc5mNmsehNoH\nx0qt2BQjjYhFdCzUIhRdUoKCfwuWFqGV9UGbBx8GH2JFJangi1i6WqTqg4mI+GLsokaltDKICAFt\nVqEUJA1JqrGW2PWce+fn7uydEBOkB87ec8733fPdnTt3ZoBrZ+/e4l8TceObr4GXl5wVxSuZiRXx\n/wvchruA4vzK7avmHysTArTvoLsq3mpJkxj30VFpo7fqdK3qyhl+R7Wt0XoKRzFTQr6cbJHQ2e0n\nCbq84Fv/6OrrqI26eLiUnN+qYw4kcW2eL9kLWmDtxU+B0bbZrTp7fnTaGNpC1+ruTe1cXdGoA1s0\nAOkYg3UTs3N7/VPNCQ1DUyrI++MjDSRKDDzUQLGmJ4h5NvA7weyYJX8ReIdj1UjnY9xp4NeFQtum\nB7xASrSRpT4x2JWAwIP1ZpwwkNn3GrAjrimRnGcAHTXpGSdTw2WgXyGJsA/dFVy4zt5VURGdFLeS\nVqyLMefKVPxai3HCwAduj4tcXFMiOY8Wh6r0k0dGXDpjndHiRvaTuSCdvIcNdXSUlenQSnErYb24\ngcdBmbX/MvAFejx0Ou1gOG8XQUwkP95sYuvyV6dKCTbp9PgYd5GvJZAojaQ4YLOGj+zj8Rj/6IyA\n4hh6GjDHUmAq8/97MHAdS+rkSrjdxIyeQNVIilux5ZvNeRpa33X2hwLjHwYKdYx7wHBUjINgnkOV\nnCc9BtWIdeypue5bcb1ahlaKWmVuvn9zOWAadaB3oHeI/uKjQ9WfLuHk2+/tZZCBG0t4qgTwyTG/\n3ORh59BNghbAsGpnX5dMZqu2c+ismibiVClqY2SHs2MBP9sAnAPeM8ADwASeBl4o2/zqYOC+7duP\nU3iBfNqzjlmTyJYFTcL4zm7kfBSrdAOQK5aZxA9KmgxTpajN+4X+TC2YUPBgeZ/wvbPBhcOLo3Ve\nJZAAfAb8ReGr5IfRdaXgwLoiaBLONoq+0Q9+l7MrNu2Zp5U0GaZKURtXeQPmSzBN8R4rLD0BXtwB\n4E/qRgAWYdGDC+O0q8s8elQSNAk/BgPdE7Dm6Z4gV+zu80+6SpoMU6W4zWjM5jV0ic7mlmadF7c/\nWpx5Gd1jRKXFZXi9o8RdEjS5uM9pcbmKbnGXYgFdlCrFi/s2nsG7113l/BcUF5TFEUBXtKdMyEbl\nyi0LmtzWqziHGTpfE/Jsc5PAFsNAP6ZKcas34jl83+cdzl8C7lUWx8BRbGWED8RhusD5OjILgiYP\nxPO4iIPyMBTF/2OyMBLoDGPNmCpFbdRvJbqh+SOR7BUXh3hPg21lYMqmxyNwhny6gW3WW8g5gkbP\nD4I3wsnW6cIHjtgegbnSoyRViloFD1DRi2PaNrLe87+VTjQHTjR3jPxNygzsedxl5Dly855BHz9v\nukPQKGfY+v7Z2yik579wCkPrnJoLQ92YKsWtVJtUEzWOAHtMLYcxw3aDs+nAOV6bRVLcSrVZNVHj\nCDDKajmMGS76nG0OnOO1WSTFrVTLCgG1EsQRsFsD0omgefzBBLMqXctaXTGU4lYt1vKxqSIRcEit\nRjHD4kRy57B7hP63IJRqb3MurVEAyL1rJ4XzLhLEvh4LelGbfwHLsmutHmzYGwAAAABJRU5ErkJg\ngg==\n", + "text/latex": [ + "$$r^{2} \\left(r \\sin{\\left (\\theta \\right )} \\cos{\\left (\\phi \\right )}\\right)^{2 L} e^{- 2 \\alpha r^{2}} \\sin{\\left (\\theta \\right )}$$" + ], + "text/plain": [ + " 2 \n", + " 2 2⋅L -2⋅α⋅r \n", + "r ⋅(r⋅sin(θ)⋅cos(φ)) ⋅ℯ ⋅sin(θ)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# substitute x to get the integrand in spherical coordinates\n", + "e_sph = gtoL00.subs(Symbol('x'), x_in_spherical_coords)\n", + "#display(e_sph)\n", + "e_int = e_sph*e_sph*surface_element\n", + "e_int" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAAyBAMAAABv8PuQAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMrvvq4mZVCLdRM0Q\ndmZ2luhKAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAHoUlEQVRoBe1YaYwUVRD+eobZuXoPFWOi4k7w\nCIkKi1cQMQy6URFlRyNGJOKIgEQFR4xXjKFljUdi3P2FG43sqOEHv5io0ehqtr3wB1EmATyi4BgC\nITHCgsuCoIxVr/v1tX0sRmVIrGT61fHVe1X96h09QBQltt0UBTlJ7Cuw6CSJ1ArzTItzMUuwD2h2\nqRpcUHJBAT4AZPJBxgbUZ/WAoNSpZFgYYGwg9akzhzvPPUgBLadfS9fRud7Y1B+qpErrXn3jye05\n4DYg9QeH1lrxCZBXfHPRx9BgqvY8EAMSItQNBRmdqkkO7SVKdMgSG5bhTJqAOIUL3F3mJ1OyarS0\nSnpqxK41xEZ+ciZJoF3nIAf5IUhmgj7s1ElzEhwqnAnRTvE8Kp78sDI5fdsNLG/hR2OTmQkdGnRs\n8CZmkJWJKW/QTWZszYm44ZiZTOQIE0esOL2ZtLdZpjEw66aNAfRPQ8xM3uV+s7xD6czZ1SUkoKds\nMkFNcvB8h6nPwf9XrJnJGh6vNUeno0ZzMzDwwVUDA0VW1Ymo7RF7GzGBRNczm05gJmJONtDWSzsy\nk7e6IucE1xiOxvN4M0ksehB46kDB2YcvHwJ0VtfdbcAOowMrE+FKu7Tm27FDebWDp737+OgOzKFy\nGIp2CgHybYVIvNF+Kixxa3HMiXAFIvcu1d4tqLvjzaSPX9W4HEcSTsHAZ/sPLGNfOi/UOfXOOfW8\n0ZM1J8JVnjeGzfeZdL7QCSMlX1CgcjpvKfFqoN0yRAO979zKRLiO4YxPVKzR/hazSsfm8lg8o4Bx\nTy+ZmtUruQIfm+I5lpqZn+m3ciKRjmwRStnkecNzTpH0UWo083vR0iYVjvZ1YJUUeZCva/44eIFI\nGXFI56aC5Ea15ArVXD6pHFJLJ1qRjLM4tJYwznI9/VZN8qdVJQesJvbNL4uYJVXK1ZJLFIEFpkCD\nIHPRY7BxEkWtF0gawBEHlA4H2sWyK1ryho6+LD93XgUqhpqetGK3WkIu0yv59oLkgI8AZQb+QFaq\nFOt6tJtU7wHizwT+fD1Nb6V5ljhH6wWKTGDHAVzvQLtYdkW6zdBtA+hKdakFINmgRP/S2fbuVZFz\nKM2ibaak0jUcAjNuyhQybWRAntXc6Vqe5NE4jAICX5ku7Mq0yWhGPYUrdpn6acCVQJeZF/DZKLyh\nSPb6GOKa2AOPQe3wWs96fKWuHEZzng00SEsHWms+OHiBBH+Dfs44kkVS+BC7giuXiYtmUHdkki0Z\nBu9zkyY1z7RJDuuJWwOFXv37pi51isnMpCvR/JHF/WWSeZB0UXyqSpyJosYLJFWBfq44HrfhTo5d\n0aQbKrFmgH4dieX3s8pUGFb7qRq3+glP/5CaklfeuWf7zWy7naI8gmYqQuKgfLu9snoYZ10xYWvV\n9iSO+9xQxk9lA6cumVfCF9vvFA6jgFDalj0RGIcL7hJiNRbpMyb1Kj7TnqcqqLjsUtj4Er8oNYc0\nYnlkL0aPRvJk2mSGO68dAn4h6VOgF+cBb1dTh0m0iQfp6uycSRrGdZeUaUofYlXhYMPo7wcC4r5U\nLasFxeFEu/l0ieV0DvECYsUf6VLT6waYUn+d9ZkZtEPG8oiJogdeIdeK+P/mSbK+pVNdUCbnAsdI\ntIkHeQ04ShrGDSC5P12Asl842DDqjICxWouWyHMcdIaNkYaok3iZe1oEPKSjqVKgEmFtIPVfrCFG\nmVD4VQJNp4+Egqj/dpLSB74DZ3Ih8CeJNvEgB8VyAuHUEbLw3zsHhYMNM6L5nl5WU29EHE4nkxeZ\n8AZzCZXm5XpUJl8MjiBGmeSsTDaXcC8oMvqpC+pFzmSiTybqITR1GLgMp7mqDBwQDiRYxNGsoUyy\nuePPRFTXY0i9RCda4hB1aZ+AVv82k9KwrxyzM6Hq2lwWE8C72AtoOeKbCQ1CF2sxi4Sz5mREONi9\nG9V1DC/zV2BoHE4fyfMaU4po1qiyE1xYASvewCcLGFeL2ZnQio+XE71k5JX8KDDXNxMeZDrmE0Tg\nBmg1xou06wkH1kpi4MfYgwsi4pB4Z0trCy/Ou302FgO7h86m7bnoNHv45AzENCotWV23Eb7arRFq\nKf0e0TGJS2tUdfEgy1N0PBq47hruUKYiWxAOrJXEwPUoxIoRcUi8s+VJHKzXf0di4V4sqdKJVHaa\nPXxm4Y47la7hT7qGz5hyGdm4VpY/yKAP6bfu113l8fWbxtdnrfwtz0pJPMjGb3QWGadet03Dc/Nu\nFA6stIiByorJt5DCNw6qnmDq85i6PXKomNZMc6ojFGcNEoETX6apGvflGwedFsG022O6yyOHii1F\n05yohuKsQSJwYGCLxn35xrElx6YAimkug1pxiVHCJBOwKRxoDRKBo1VIGw/35RtH6pcc2wIolXMZ\nrDFd2kBhtWmRGQUArUEicOJi28yd+MaRyOYCBhDql13GPS4pUjBq2iyIELQ5iFE4ITg6SwzyjWMP\nZ2LedMM6aXibWqVM5E234aMNC1DcYuRNNwzY8LaH+T72kHnTbfhoQwJUC5yJvOmGABvelOzs7Lrs\nx4PmTbfhww0PMJ7jbzhx0w0HNry1NWfddBs+1tAAE1OGa/KmGwr83/gvvYG/ABH+GL/r8Pp3AAAA\nAElFTkSuQmCC\n", + "text/latex": [ + "$$\\frac{\\Gamma{\\left(L + \\frac{3}{2} \\right)}}{4 \\alpha} \\left(2 \\alpha\\right)^{- L - \\frac{1}{2}} \\sin{\\left (\\theta \\right )} \\sin^{2 L}{\\left (\\theta \\right )} \\cos^{2 L}{\\left (\\phi \\right )}$$" + ], + "text/plain": [ + " -L - 1/2 2⋅L 2⋅L \n", + "(2⋅α) ⋅sin(θ)⋅sin (θ)⋅cos (φ)⋅Γ(L + 3/2)\n", + "───────────────────────────────────────────────────\n", + " 4⋅α " + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Perform radial integral\n", + "val_L00 = integrate(e_int,(r,0,oo))\n", + "val_L00" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAP0AAAAyBAMAAABlpcCuAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAMrvvq4mZVCLdRBBm\nds2z0OeMAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFT0lEQVRYCa1YTYhbVRT+MsnLzySZDG5cKCS0\nOCAKFnWhiGPQ2RRLZwSdiqgTSqnianChq2JQN2JxBqticTFx4UakfThd6GoiqHQ3IygoKJOFtIKQ\nxDLW1tLGc3/z7su9Ny8yB+bdc8733fu9d99999wMELXnokHML9ZjCTPMD7j1zexkUdDw8c/5QGRa\nXjgR+LST9SchpY4TJuA9H5gMyz3i4r15gyF3uWCWf8kHJsM2Ww7elye4/uaqA6d0cd2NJUXmnMQy\n188dchKQ9dybu5eB5LiIkVKB0MeHKh5tvx5NTZopNJw9pP4bTSfjeWCXPsCrTsJ4YNs9utR33yG9\n/uCpU6d/dg8xXv99N0Xql5xvqNzG7Xi2POseYiyS+tdNkfrFmzbK55QsMKCWCW14wlzWs7qlPrqW\nscrzlPyN/lJ1fhdRyvJciFzPu29oerqu3RFH6e82RyC8u9sC1imfbVdizx80Mx1gm7AEVm27SUp/\nbcfC2aqhtET5rbDSMuFCM3UNeNxMuiLbsymu0k/XVCbSZq+hwB78F6TZbUQs3crTono0kvG4B91Y\n9onrRzg6dchCyl/Fdyx9LzIxfYD2tKJe13canb8R0W2P7S0c/Jv8+0VcWrwltAyyDAKazlFbDD+J\nJPMnD8zqsLCDXF9G+RoEJiQyxCrUgWoNoMJbVFtXpSP5lqaknyUKVo81IuEFBEPWMhB0JDgTQmJC\ngvJrdG/VOjAFTKtv+/yS5OMr5eg2pUg6w5yZ+TASPwk8gFcPkIUoNQhtINVi+E8AxwAhQfHWutDP\n0pv6h3HIVlq8ocsp5ehWT5LOMCenpo5nHwYWZyX+GT5FZRUZHj4ESExIbAJsItjz54DyX7JPT7Y2\nfexpMOr8EQ16odZPXX79HlSb+JHh0/SsEhMSM6tgew7TJ6OvSNgt2Vr1r2jQ52zI95EeDPrBxsnD\nfD3IwyXDhAQlWEGL6U+zb0HY6PxjniNBl1tD8uJNZIgINNVmAcMkXupgaknrZ24IbmTxWvTpBY63\nQs3GKayyLMOkRG4dmfZQv89gWq6sDblr0T/DgTGX41Y83WJphkmJVB8B3ZOc/4zUr9TopN0EXux2\nz3S7l4iOCv9lwTyM6guQM9iFRinXOTV+4fockxKkzw4MUj8r5//8Dq1F0dfy/Gc5EpzlVo9LiPg1\n5JsWhM8/x6QEfRCpUOurrXWFbulX0d2in+D9pxoo2/TZ+hOYlKD1x4ztv2Tq+9+g9SGnwqIv1j/v\n4Lq888zyYRsW1AGBSYlsA5sdvL1x5WVGF1t78ehg4eiAmMws+gm+/95gcF30N698/yFMS8y0eP2R\nrNywaKh+o/r2/Vfxx7Qfx/C3jFjXn2H2C+XqkmqvP8C3g9/bwAeKb2+/j6VfMGLfo+mSqouU0ZUC\ntqPmB/GsGU81jbjYMUK5tZo5GcmySStTFYkYje+oRetrHzLp/BG12O3gviho+vTViZKqNikTVjsq\nmwSfnTbAi0ZEP99oL3CYLqnpmp3Bd1Q8aAcTZqurPqIoqWs7dk6lw/If2cGEWd/vD1Uyd5v2wcSB\nKtnvHPsItOn3XQjlCzUO3u2giAPVEQeaLD09ugENOx4X7s1hxvDYgaqIY0Zu4sDz+mRJdf7+Ztvy\nRfwwsaTRYdv9AciSOlM3OuiA750N3KET/8tJN1zdRNmkc3vTzgjo7Fxu27HkWff/n1RJnXMMVqCl\nuxIWExRnxwAifcmFypLq+v/khd7eQm9vbP1xja7za7PatTqbLWt635J0IPaaa/q9nSYBT3jJpZoX\n3gcws+Mb5JwP3B/sFc8wRToF7L/9B35pZBo4+USHAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$\\Gamma{\\left(L + \\frac{1}{2} \\right)} = \\frac{\\sqrt{\\pi}}{L!} 2^{- 2 L} \\left(2 L\\right)!$$" + ], + "text/plain": [ + " -2⋅L \n", + " 2 ⋅√π⋅(2⋅L)!\n", + "Γ(L + 1/2) = ───────────────\n", + " L! " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Replace the Gamma function, from https://en.wikipedia.org/wiki/Gamma_function\n", + "gamma_half = 2**(-2*L) *sqrt(pi) * factorial(2*L)/factorial(L)\n", + "display_eq(gamma(L+S.Half),gamma_half)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAAA0BAMAAADGTneNAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKmklEQVRoBeVafYxcVRU/M/PezufODNR/Skw7\ntFJDBLulIjGSMiFCJHzsREwbpboTAxuwhJ1WZIGq+xJNjPyzLViooum0RvxDsKOkBFzNDoSUiLVd\nIdGAmJ2ikZgauu2229JKx9+59933sfPum5nuh3ycZOeej9+5573z7r3vnvuWyEPxlZcQ/W264lF9\nSNnNtIkoO/UhvXvvbe+m/jz19HlUH20KmvfkxC1PkNns+tmKxZd30nKLknVP4Dtd3mi4PEZO0St1\nzW8J80jUwqwLZ7tw5afczidLNGC5YsRzTUuIJPRfzWcaRLcD9YSL7I6Ll3X4fWw4orMuqN4s04jl\nRHiSaNIRiKITrrCPFPQsK9NEkaZr7Y57raTBLznJhlhFY15QdU+JkttUhHiBaKsS0H7Z5TMVsqGJ\nU6yFbJ5xzV1x5k4N/NZVIhPmTzX2BVUnq9Rzonct6EqiNxHqN0SfVxEvVwxRMk8CShQ/zVoT6ROD\nw4V0zMUKOmhGZIJe09kXUp8+iUzYARKVRJHeJSrbslm1GTQPYUJIaHpKaH9P9Khr7oDDkLIJq5GG\n7EykKhr7AqtT8kEQ3fj3V0vGMcqU7YDxBtEhvEaPQ5arPUNzNWGGYo9gOvxJPJBXSL2fnYneYwq6\nuO2oegK7mk168PiqvZYdH8tE/N43lj3Pt7BO6Bhqw7HAX2XjOmtSKhNZ/X3amegux52F7wT1Mx3o\nINHNtDlTZPtTAsRQ+1XzOtE9Qhfy49u5O5lI9WldVCaOlrSQBTREgy4sYSFilaP2yfktlnsBHYfS\nJOon+jbbw8i3c3cy0W9pfVQmBupayAIa9gf1/YcaXpR8OUY5JewiEwI6DfmwyMSLwhLy49u5O5m4\nVO+hMpFG+K6o4+IxBNhboAtagkY24lUpUhBt5MRI5dkhoIn/giuIl8l1LX6zFL6du5OJkFeOykRG\nv5TMCmGLHRePIcBriZ5v7d3Ek3+L1QOlnMUtr5gCGsc7N9PobEs8YMHLJuOQvTsxzilVa6syEeF0\nd0OtxaPGWw+M/GJ4VTXAa3Wefs7qlylZ4fY+7K4FNDWFRRPDZAWZj7MhjLw7dwcX9rxVJugRB94Z\n01o8avz0wDT2C9UAr+WFrMjAOuoRLe+sGLp+fGZwfAYOzzoVmNHw+meLrrTVZV0uts3lZ3HRu87Y\n76PPzLJIkSOZbwcXwt7i8RbAv9oIxtFsIEWI/hkYTip7Tqbw4F1yZrmtirj341Sp0oQqVe3eeeeu\neKRRVmzJPoHLDp0NeQuP5FXoG+qKI0Ik+u5XCqIQFlrjJ47xSbd4jCBC4spXZMHsABQzC4g9E1FP\nUVlbW+MUevJQtuARwMbrjuxWqUKVdgzenTsZw2ukYaAg21zNBbZwRy2l6q8ojmgf3mU76CQXwpIM\nUQ0y7y0e0yWiG0pY6R2cDQ8CikxQ2LU8h7Besu9DqW5TDFelvoIWsk3enTvFuHRhGq2LRu1VIZh5\nqfH89tc9gmK551QD+XVDKJOveHwJ2l9RbiII1wIk2gg4u+hosuC38Mj0kJuYpFOlSjNXqZK8O3c8\nn4GSUE9OSOuIJVtMprrinFYNHEfBDCLRv4nOiULYZyFf8bgdZ2rbKNcIwLUA0c338Pclf3c+KY2w\nXoo0vFLWtfJa6iloiVClBhL250yHiqIh3qtKCshErqyMHykqTgyqPWSgWlYhIp+zjd7iMVHF0ClQ\nPzKucE4XviqTgaAK/tIAz5lQlIK4St3ygsWsVDDnI/PXUlxtJ/GsYw3KRB9bb3rjj5FrysZzHztw\nN4voOHGaMtj0cQjjTwdqS2boxodv2l+nXZ7ikZeMUYsOWRJnfnx4gjYcuF84tADJKF78Z7nMoM+5\n0Trhjgr1IsvcHkFujwT3t74k9XYm5NGXUAVkIlmDxezDPjdWpvTVtDwPGZF6ZgbvmpIhbuXC6NNE\nT9cjx2B0KdYgGhoc3AUNX8rSCWO7sZtideHgwnBKCCBdHmlg+Gc53lzpKdEBqtTHiK7IlIjsWTCr\nX/PiDVIzBAhIHn0JNigTVVgSO/CGi5UJOxBMelEPp2rifIRD/KCEcY1MYOtxDqJLqQmiX8oTNcaN\nUe+JVIWME8LBhWEcAxhrZPPxMvVWvYbz5EVthio1ivnxCbybuUoNoFTzjNQOySY9hVYmJSATqW2M\n2nt1nmLIBG6/DhGRchUx/zlEahpVAjLxWSL/3jxpEZ0SywlfinkcWF4zTgkHCA4x8GUkO1olY4p4\nw3P+xJ2KTOzHyt5H9MA3oQjOBEMF2ZnIAZ7NY2yMjf328bGxAmw5cSECJDOxYfw4xZCJPicTAxN0\nKQAcwtzaLHAm1gZkwnyXotskTlSMkxbRtHCA0iHOxB5kIt3HmZg78ezoRUE7UEFBMgHB3jfoel5d\nEpbROr4iSEzQmKjCFMnTUSvmZgKRBiwxADjEBZQ9HZgJDHrztMwdcM6YOC4c4OgQz45ztIyWT1Ci\n6mjPn+EVk6tUnsoj3M0R/tGTvZEeKRIdlKiATCSrMPViz9aIuZlApKQVZwuH+CsfmgWNCV4Id9KD\ngAjcGDpKFvDWEQ6sVcTA3+Gs5QqMzppSzqG9z65So2XKrjbwwFeEd3aoKOx7MTGwsjAFZYKvrHcH\nxfKYGmp2IFK0vjQPC4fA56M1PDVaZgcWQLosgu2VxC1t0Gbjh5SuCAfWKmLgQ1SJFdAv/uZMGIGy\noB1+uxH9Frp7NrzLoxbs5qbm4KZmWSIDMsGrCCW+c/B+Y2jmi0MzN1/zKGQe65ddwj4c4vv/ecu6\nvnnH9c3bX32nzEpFPNK/8AIeicSZX3spTxcO3ykcFEa0DDQ+ue5eCGlLaOb2E1KlBnbMq5efAjIx\nUPBDhORE8hTCATDarZRtcAIo985LlUtLe62rkV+En3AVs7hswa/wVKl+gy3xUumnRMOWIyvWFiXb\nX7dV3saJ1CbEm8qpDU5UY/z+IuKxHEzYKDl0FlzY1+A1DlIwt/nFFinwcUvUeme/NTnR4geFitQm\nREzcHPBtcFiFsFhxILPGv0GUecTVim1x2NfgJS6WOXW5fq0rJbVR6Q6iH0mge1LjOsqTGpbbhOCT\nGkFtcMTADCOd3Ek/z++yp11Bbot5YGhIzjRllKNNSQFtbCpAKVU/RsFQFOxQIMaO1DbEMundFoe9\nhKTDgeFYWfFkIi2unNfu+SHn1La1u/GSyoRnTLbCFlGTzXMm7CI7V+PIe+YtPJ8v6GlviW1dn/Lr\nO5ybZSMhE6rIHq1wZ1fNrUevd9jwsmv1sC8B3q4WnK9wJh6zi+wRi+OFnEN3ezmrQxxSfcIYq4Zg\nFtHUW0QmnCKbT9nM9l+DO7++UUuPXSlNuboespiWDTgIcovsaYQ+TC/O3wWkCtq+MmVpOipWCy1s\n0Qx/GRx8525VZNtfhK+bv+hZu/AK6PEVrsRBzwTY/k+qR8R5GRfZ6ovwPF7Iw7q+jAJlOBOJYzrE\n4uunKWoX2akp/iLc/mtwF9eIL5PBdNHwlq+zhevk9wh9o3kPySJbfRHWV2DdX3K6oPEZb8rDTu1/\nrmr83rfqiPq3R80daP9zVYN/H6tfD7/299DkCL/QuVuxBoXRP8KMHzCb+Linu6dETWf5AOqjVshN\nqc+GIZAFMP0PEgj53IgqqawAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$\\frac{2^{- 2 L} \\sqrt{\\pi} \\left(2 L\\right)!}{4 \\alpha L!} \\left(2 \\alpha\\right)^{- L - \\frac{1}{2}} \\left(L + \\frac{1}{2}\\right) \\sin{\\left (\\theta \\right )} \\sin^{2 L}{\\left (\\theta \\right )} \\cos^{2 L}{\\left (\\phi \\right )}$$" + ], + "text/plain": [ + " -2⋅L -L - 1/2 2⋅L 2⋅L \n", + "2 ⋅√π⋅(2⋅α) ⋅(L + 1/2)⋅sin(θ)⋅sin (θ)⋅cos (φ)⋅(2⋅L)!\n", + "──────────────────────────────────────────────────────────────────\n", + " 4⋅α⋅L! " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAABJBAMAAAAnEvkVAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAALVUlEQVR4Ae1bbYycVRV+5uPd2Z2ZnZmC/AAi\nHajiF8pCrcTYwKRipQJ2gqYEqu4kQFUw7NBEVqClk0gESbRbktJKVQYI5Ycfu5UQxYbsgKbEUOmk\n/DCoZMePhBgIbFu2rRRcn3M/3ve+uzOzi3bTnW1PMveee85z7z3nvufcmZ05C8w3WjxpaXy+mXZ8\n7OkdvaT1Qt9orVogmjdb+xEZaa1bIJobHD8iJtCNKFp3dAuT/bHjVlfFGQDXh0adP4ict7Tge3Ha\nkk8C3lF/DJzt8GQ/DmQGjt0KJEthRYeO1iDhe+uVsKmCXvc2p7sOeVUOspL7i10Qxx1KNwE/sqZ3\nFdEzhIR4Z8irWk71iQa7zWU2/VU2nU8PAQMF40ZPFV1vIZVHvGIkU245lfYMDxMABtTB3Wgx8D41\nId5n6+iyDk255V4U+ag0PSVpFwINFwMvkhNYncMeK7ibzPrnKmy7panyhWPSJPPSLgDqPuI4sbmc\nGD7vZnsPelXgrIq3JVLHH3gbpGtEanhMkn8hULJPvDDP/7GQR+kGsA24OF2MfIUnkhSlfovoomZB\n0BLxwlPxjKg6CeXWfWzpbnQC+FiCgEPAq6JIjbMpJurCdz6lS8qHnar1Ex7pHRTs5fXG87jj2+Qv\nzOFxwWQpyOTSBeE7n/YjkrNe9OaxyPBn7KuoW66/DGyqU7g4nyHLt/saAwLxovAdT/E80jnrxQrg\nWcv390HczTbovci6JpLK400FQL3zibDT6azB9TdbHyJPDJ5ftYPoAYi70RIyF8qTjh/Zr1TDvPgm\ncN2IxXV0Pzo5+bZ1IMU/Z6t2EDmob7nBNxrRDSJ8+ik23i2T626ZLC2Uv3Kss9P7gaK65XzFWN5n\nTwJm9cawu6ncSeC072JqB7P9pKXegyet6+L4rpPa+1POnzqBE3YC9mejk7I/Yad+auNTJ3DqBObh\nCSTm8z04PtcHpr5gnOtN5u36+puGeWveHBvWITUU2dYJ+v8cUHXWk/nT5ImjuUnQ7tJsPeq+40R+\nKTE3CRqtzNZ7JE+k93OToO0j6i+HysHhnFDvq4Edx5Fb2XatzLijfs/er7+4iLZVa87iM7CzT9AZ\nFgqrN4aHU0ZdfY7gvXqfyHWNAG2q1py1Z2JbJehnQxMvbPHGEAK5g/B3za5G+J6aI3mv3idz8QPA\nDc4K/zvbIkFX3OwuGXGflatoxc8wob8STIzvCxceBZoWXE8lwh+r3aq1FkDMnCKtEvTP7pJTiuJc\nleEj+3R0mKH8rt6GxtroZqHqnQhXrbWYMosUaZWgIe/PbrG+ES/jT4iVECRV5DDecGWZQjC6K2Ah\np+YMffaL5L7WwD99QcAka0HVmmzivYFmuHCKaGtycO3QP3kHC/tcyPsPUyzWBOQaldk6rWhwrSBP\nB1R1oZl1I9C7lMSz+hVwheWDNRF3wllyp3vZfnQVLGBVzXLr+TslLz5F3ATf/WrewRkFb5dQiggQ\nKAO0w5JsoosbrcT0rvcetyJQV01qtDJqp8bewytuSu6+Lpqn4KnqQrNiyvTs3gVKdhQfvMhnnQId\niZ5VxWwR3FvT6rJhMnkWaNiqNf5yGd+KCQQ4g5LOTRECSVzDsUMlaNbfQSFU43ovRXG0xlRNajSn\nRHTEeo3H4VXVJL8Rh9Jl6OpCI+XYEO/sdMkOYrjfsm7/PAe/lOI0YcJ0L+7xq9Zk0WRDjnM6jpog\nRfTuXk3ZZddTCbrZt8vLWYXr/fUUcvGbdNWkRnPsva3AaWbPlKJBCRb05NCjqgvNkt6QYXDnwfOH\nK3aQQn/R8k6/hTE5JKUJ1zlCxcafeWWZX7XGTfAv4D+YjiPYSREBMpcKNNu3AypBN1VEI9RbUx0b\n13spiqM1D+mqSY0Wo44p8OcHtxXkfBzqLXPAR6qrC63iGctM6f2tIpf7mu4qn1weq+tI1Y3wfQXN\n9ExOjvtVaxI3jyLOZ+/jDJydmyI6wBKiDOxQCTrqTwi8f9mXwavy0PgyVZMaLZs9CE1bwfNJy+Wd\nR0aMVbcCD57E6kJTTihPohl5u0Qa/+PekdMP48oHrtpT4zCRZwFOBfsqisNVf30hclkp/vQH9n6L\nSoe4ZvdRpPkBQGbgmsGN8DY8X5EJgJsiavPLo+/PqYgwK0iCmmconO/9acPBNumGWRsYLlq0bPao\nTAGuPXhGFRi8+tnXF9HkEiVJvnCpNFIybMoJ8ZoSTGvWcE3gS1KC9CngyVrkAIexBgNt3bqH+QBH\nGKx9XDJWQmo5FueoDYibdB1e981xjYv8BKtyX8hhuZoQShGxJlo/R3759O1QCerUQvreB+trX8Qa\nHjMvZYMWoy6xKJ5PvJ7I/Z434WtDFO4XxS+U9jFbThhOJqVTjXfutdJ/v8j7mN5/WiUxknXgZ+qk\ne6vccysSiJUQUzcB0T5xk+SIqkkXXGoEY4UngNVXyIRQiog1l+JvuN2xQyVoUP8ePHuzvCf9Xr7E\nGrZ9jAJdJSqb3Sqy+0QuTEyuxXjiLbI3yni7NNE+v5wQq0UwjZK6LCd5iMVY9P4zwDvE9FSAIyqf\n4+McDi/PIUbv6WiNw4C4SbYs1wME159nzzsgW5YJISIwXeUBf5RnQ0WQoClZv6iw9tmrr7v4JHmO\nuvRXrAGW8GXQstlGDm3RILCSocvNxXZhtPd76EifKSds7r1ASd5dk3nxfqnvvfcuokPaK1w7ehAx\net/XxPv+Oj7CBcSgzTUmL59Otk8mhIjeJ/P0/kntfZCgWZqXyfGp7t792x27d+eDWT98mnymzEZ5\nny6RM2jZjJEOWzQIfEI+cwA7GKh5YSTWelld6JcT3i/SVrQImaOO94w1/gNKtsbFqnwKObxZiTXz\nnpv0V1S0CM5/9ufKBLhEYLYRaUQYGbTDSVA5sahG2mePxJIPimSswSMrkqE1ks40w6Bls89RBlM0\nCG8CPHXeVUVEK8LIPbOC1YXZhiknfE2krehPkkjBs5d7ZjvuJFouGCZoVyPWzHtu0lNJVA0uRWz9\nAcbAD2QChQERmKplcrG8uvWcBN1U8Isbfe9vwy0ys6cK0C6oOzieRzoHgxajFJmiQcmqbSIYq+tr\nQO4XVV3olxOepyc0b18u4iIJexP5iRJwQWQLsdE8vd+KWI5hPz3yuUm0dmbO4CI/RSK3qoBfqwkU\nBkRgdzWBiykRO4IEHWbQq6DlNjWDf0Q+TTGr+TjVB2SxRldNGjSN8hjk9NAUDfJpfV3G/SXeAEIM\nMF1daMsJf6PELZp7X3+1snLyppWTN77075KO92ueKxKcqnB094vfiQ8c/vLA4asvezC0ADfxLlBh\nKjisueA2RD50e0UmhHDq4/SKf7xaoFTZYRLUFjcqsO/9diyuUMKyb6YuSQJ9lFWTPlo22ykahvrj\nqjdNakjdhkwWOT6HIkPOYEb2EYs40zJN+2CT9jhtzTmyhrbDT1BnWd97RrAcPMu+U9oJ3xqDDjab\nUjTIt7wRhcnkDdR0iVp43H70d6veYJmmfbBJe5z61Mu3e5K2w09QZ9nuhj/4ueIWl7+nJb41BhBs\nNqVoMH6EfwwrushATbc2PJxhFNOHDvV5rA3WbjITjrcK+IZH0nb4Cdp07UReiWPju7TaWmPAzmZT\niwZ3pBsadLrBms7aGZa2GtkjnLLvNLjdZCacfNeCmkzXdvgJKqJpZJ519+GSVllrDNDdzJyPXWJg\nbdFMaViR9BnzMF1ZO17lKHSwtsFFGlqpgroNDhZo7PATtNmc7nJ3QcmHdQcYawy4zWZjLxjMvO6C\nBG1i5pWvvKSf4O+aKNuL+t9pr58nWj9Bm9jzcPMvWpsgp4lSB6aJ5qPAT9Djaxw/9HUCzVGCRgqd\n4Dw6JEHn6Cw7JEHnyPsOSdA58v44JOh/ASx7rK0dLSpXAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$\\frac{2 \\cdot 2^{L} \\sqrt{\\alpha} \\sqrt{L!}}{\\sqrt[4]{\\pi} \\sqrt{\\left(2 \\alpha\\right)^{- L - \\frac{1}{2}} \\sin{\\left (\\theta \\right )} \\sin^{2 L}{\\left (\\theta \\right )} \\cos^{2 L}{\\left (\\phi \\right )}} \\sqrt{L + \\frac{1}{2}} \\sqrt{\\left(2 L\\right)!}}$$" + ], + "text/plain": [ + " L ____ \n", + " 2⋅2 ⋅√α⋅╲╱ L! \n", + "──────────────────────────────────────────────────────────────────────────\n", + " __________________________________________ \n", + "4 ___ ╱ -L - 1/2 2⋅L 2⋅L _________ ________\n", + "╲╱ π ⋅╲╱ (2⋅α) ⋅sin(θ)⋅sin (θ)⋅cos (φ) ⋅╲╱ L + 1/2 ⋅╲╱ (2⋅L)! " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAALAAAAA7BAMAAAAp7h12AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAFkUlEQVRYCe1XTWhcVRT+5ufNTGYmyVTIQsV2\nYiW6UDL+oas6CCVYsI4uUhLBzKLtwkAz2kKQ2mQWIqKQtkgtNgufSlpdZRRUyiyS+lNBigluhBpw\nLKIIItPURDHa8Zx77n3vvsmkSa2u6oF33znfd+65571337xvANuiP+8Gem+zoJHGZSsKUDvP5i3q\nym6k/A2ww87p3pS1wgD1YHvFoq7sOhgFPrZz8vGCFQaoQrzkUaG5hrG6B9pOXw7YZgOYsKMAFd3v\nU7Gy77f2XgGOCpOQU8XOC1IJn7vZzmrhp9GZCdWRZiqieOeIlRakZnHR4+7yvNbOAB5Bew4HmJXC\nVlcIUs6RqLeo4/KMlha69d48zTy3G+c/Ouhyii6cZV9bkFrwt1t4HugYWdkHJIsmWZ/7Ef89CEnh\nIKajVdQgE518y7fUdY457QFeB+wto2Z3EuKlcqCiVdQYlzlcomHIZdeyE8BI3oqHq6er1ZoF+O5q\nynGZHS/ToNrmyNhMQRfu258RTF9vpGBSgOjo3SpoptI1hmd4aCvyGLRpVSKXcAXWs39iNH7nHQxG\n8LLibOolQpIKXeExmVW+PSR+U1HFWRJUZid2UOHQOxjI3AKkMKQWt6j0JGWf4xkyP1KS2daYzEnQ\n7spZZj+VpFJtJUSynzF8XnEW1TVXBlxGZVfFauwHbKuOnsyIE1anWS7cW0C4UqLYeU+BFoWhHDqY\nQqpOQyE+rzKsIV2UwKGN51t09Dl6DV+lfo4XCO3nQZtQCF8Er037IUcvSSad17R3+gqhDAdPdJU8\njJ22Ahy6+/G/yHe6d9kcUXT/F/GjAg/PAmFEGbMtmkVaFZ5uuDbujHcDl6lwndBk4w+LUxRt/8KU\nAsfzwJcWrd2bRp99ejWqkR7gQn37GvRjh7KKmabVlzBQaUqbaQS6CbLxsV9w+2wQ86LUJF+9M9zY\nO9worv4R8vKu2mlfvOopG5wgW3CDyf+nXd93gL9j/4Vd33f137363jUe0FqrhNbI1zD/OisL5Yy3\nwfP6glcKrZ+nWzStrCd4TX/r5N0PNC29nuA1hVnI7DSBOn9vRR3HgODSjosbtt5HKT80PqxZmcpN\n8SdDzKnQ5zQHEcEiZGN5ok4K/QJ9DoMthuedoqi8FUmxxtBDfmEWxrSOFsGiCHmthsp3alNwXOWa\nYRCxAtqOGOGk4MeFdO4Z8QsPEvY5sEeJYBGyHDvyNU9jCry0ZWNocxH7lT68vhxfMPycX5iF8VHg\nhIjg8TKnDNAh19k3ejzPS/vmuEgtqcJKOAnRojDlgdWpFsFKyCJFPb4mU3AMtDT/B9RfXBG8SdKd\nlpBuUZjz4llVZJouQ/pk4C2FYddil4v4gYXNZzKAL3hZ+/OhrUVhFsaRGvMsYrUQ7qgADzDGRks/\nimdY0fmCF28TIbeNUxAs7DDEwjg5zx6LYP08WAfvY8xrMUc7Ab7gDVOu0v+qRlPh0Cme6tLRVqYB\nLILV8yggWgcOUei1GC0qye8J3rM84RId39FxsFp9t1p9nyHwrpg4TU5HiQZVWIlgEbKq8KdEeC2G\na500wxO87VlsQuJPQrJ0kAVvxbc1un7OV7dCiWARsmqbbOcJpsWhQmeZQiN4HwbOIE5bWbZIc2Ha\n5/ia0tXDi2ZZBGshyw9PmdHkX/B/DTIRvKFToz0uknV6ftwXWbDjNL05HzAcLwIigkXIIpyFM8mM\naXEbYqqwCN4U7Wq3f2Z578wyJ5GZwj1vnKzRW3sJdLPI1AtCItgIWaTK5kfIaHJOIxPBK749msIa\n682kMsp9004i/0YvNppcA2sJ3glvgnK2lF4U4EIQx/Ne3NzixgRvpK7zItK4Kcc/0drWatHwrc+J\n5aIQTd9Te52Ntdhcfzqvkc0Bhl+na7NPrm36P539NwwXGDWfoeilAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$\\frac{2^{\\frac{3 L}{2} + \\frac{7}{4}} \\alpha^{\\frac{L}{2} + \\frac{3}{4}} \\sqrt{L!}}{\\sqrt[4]{\\pi} \\sqrt{2 L + 1} \\sqrt{\\left(2 L\\right)!}}$$" + ], + "text/plain": [ + " 3⋅L 7 L 3 \n", + " ─── + ─ ─ + ─ \n", + " 2 4 2 4 ____ \n", + " 2 ⋅α ⋅╲╱ L! \n", + "────────────────────────────\n", + "4 ___ _________ ________\n", + "╲╱ π ⋅╲╱ 2⋅L + 1 ⋅╲╱ (2⋅L)! " + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# use gamma(x+1) = x * gamma(x) and the previous formula\n", + "tmp00 = val_L00.subs(gamma(L+3*S.Half), (L+S.Half)*gamma_half)\n", + "display(tmp00)\n", + "norm_radial_L00 = 1/sqrt(tmp00.rewrite(factorial))\n", + "display(norm_radial_L00)\n", + "\n", + "# Remove sin and cos terms\n", + "tmp = simplify(norm_radial_L00.subs({sin(theta):1,cos(phi):1}))\n", + "tmp" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAADYAAAA1BAMAAAD4/wVAAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAInarRM2ZVBDdiWbv\nuzJCz3LGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACMklEQVQ4EdWTP2gUQRTGv937N7fn7W2nIYSc\nCKJEMIWFouB1xu4CFhZCLp5/QDQJKQJWGQwKimIKC7usKVJEEq5KYYpLI0FROMFOA9uIGpsTjF5O\nZH1vs7M3kNsiZV4x73vfb+bNzO0c0DV+HPG6+mx+zJQ11u9zKENarpKUT2kasJe00tB7kJ/U6vSE\nNhEo4menftiRpMSavdYxbsK4seFF9aR2BzGLGVi/IqYLq4bTwDtlFeg6TVXQdm+A856qN5ccJTEA\njEiNNZKzzJIlGgKFOYmT1O03UBYtZifKgFlklWzDWp7sPcztMjzbuEenyzJCtoGveGF6rO/zBIh/\nwCsWuMNDIyFpFHRqjnEHzzmbgzTYg0GPB4dcttBfyQViCoaDdK0gyZvzeT8g0cpKSnYFpoOnslCi\nQoXdnmLZs7hwBjiOA64CnK9c5nHE97eBT0i4XKmYrii1O6eCu+z2941D3zIu9nyHzMSF2DU5nz9l\n98jJ7j67udWxeIjvsSyFQizL40vE+Ltb88tRbc+PKt3Xogd5CddVreUnVWIpFxuRxz9iWJjE6jHv\ngtnnaBEJe/GDKpltbR4t0YGuXmMzj4saE1s1rMN4hsfOQdqcnncIaZ3wJW5LetL5yjF274YI3JOe\ndL00LpEuu2SLYZ29JlZ7Tzu+lWTP8BBEeM76UJvYH7LE6KMQBT2naT/8JdYkN8v/mp3gdXlXrKMK\nrDS/KZtz+tz2WWCh6sF6uYpbRZ3F6v8HgK9qIjhYdQAAAABJRU5ErkJggg==\n", + "text/latex": [ + "$$\\frac{\\sqrt{2} \\pi^{\\frac{3}{2}}}{16 \\alpha^{\\frac{5}{2}}}$$" + ], + "text/plain": [ + " 3/2\n", + "√2⋅π \n", + "───────\n", + " 5/2\n", + "16⋅α " + ] + }, + "execution_count": 16, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Integral for general L doesn't seem to work in Sympy, but any particular value of L works\n", + "integrate(val_L00.subs(L,1), (theta, 0, pi),(phi,0,2*pi))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Sum of primitive basis functions" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIAAAAA+BAMAAAAR0A79AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAiXarEJlmu83vMlQi\n3USF4r00AAAACXBIWXMAAA7EAAAOxAGVKw4bAAADCUlEQVRIDe1VS2gTURQ9SV5r2kmnofjZCKYt\nLerGWNypOIgFXZkiWOkmwR8UQYsU/GwaXEhX9YMfsIumC3HZbF2UzqJdCaYLN6KQbNSFKFVrtZp0\nfJ95bybp9JMZV+Jd3HvuuZ/c94bcBzAJtR/g1rc6jX7ftbywD1eCNQDOBGxAngRrQHYngjVA0EsE\npvJBRqA3UMgEaYC3GDQCNZhsPxSo/h8vntxhBjvhkD4frEGWLK7f4I7lyFevVG2DCWasHC879vLc\nxLJXg+OmF+twEUv9QmPR4w9BNlx2xW+qW0NWQQV6Z+MKe4MpK6MCAwopUHQmVFw1aLTKiggrxMGF\ni50mB23bYHDgqUYqnjQlL6+0lFispRRJxBGlKDyyAkwvDzNWyQkrpXAViB0dEn4/GsdKIMyJ/qDq\nlaCl1q3fEkqrXaMyjEhCEOPQH9qh9yMmUHutu5YNO1xjmnKCKKN5EVoHc3bOpNCcF7TSUSurcBUo\n2N4YQguG9o55yZYFhA0K3EImvrhc0tZl2u4H2851Z/ryuMm8VHgJDYIOte2145ixShICn8xY7Rl5\nkMzTJD1PKgjRr5FC7A1GmWXSYD3glqnYOHpzDKySHoM+4wbSpc9As4mmOA4yy6XonCHyqOO8IL00\nnb4wtJVH0oYroeB82daki18FaYMtZTH2Y3dwj6m81riCHuAjvYbvceidIGxoarnov5xcejQkHLcG\nsemLJYROGlgCNGaZuB/GcBa9Jme91Dwl9wGz08ApkAFmqZD73Niqqzvjdt2Y3KjQWIRSV+kH2X6W\nW+pFk+6szWCymONptr2+mZqqHC0rrlrYkPMNtaq0tR2y3+RBYQeFw5hLnK5Tub6hfrfOWp5eyMsq\nkk5JWIcl92TyiyP2GyOJzdnoz1Eur+kLt+Z2Xa9V2nkaLc+3cb3i/7G/dQP2YqXtnvlrKRcrnh72\n18CpEuvF8etG/hrIxUp/zl8Dulgj7M9w228De7Gy8/qbgC1YW/w1oIuVH+GW3wnkgqVTPJej1GXl\nggXmij11VYrkP3X0x7vQxbBeAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$\\sum_{i=1}^{3} e^{- r^{2} \\alpha_{i}} N_{i} c_{i}$$" + ], + "text/plain": [ + " 3 \n", + " ____ \n", + " ╲ \n", + " ╲ 2 \n", + " ╲ -r ⋅alpha[i] \n", + " ╱ ℯ ⋅N[i]⋅c[i]\n", + " ╱ \n", + " ╱ \n", + " ‾‾‾‾ \n", + "i = 1 " + ] + }, + "execution_count": 17, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "c = IndexedBase('c')\n", + "alpha2 = IndexedBase('alpha')\n", + "norm2 = IndexedBase('N')\n", + "i = Symbol('i',integer=True)\n", + "# Fix to size 3, could be general\n", + "cg_sym = Sum(norm2[i]*c[i]*gto000.subs(alpha,alpha2[i]),(i,1,3))\n", + "cg_sym" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0.15432897 3.42525091 6.36113067148303\n", + "0.53532814 0.62391373 1.77361123608569\n", + "0.44463454 0.1688554 0.665504872778405\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA38AAAAYBAMAAAC4gFnDAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJdjLNVN0iZu+7\nq0QgoRR7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAALwElEQVRoBc1ZfYxcVRX/vfncmdmZebRNhVSz\nUxYIH/0YdlcMpE3H3QIBbBxbgxglHaQuQvgYFYQoyEZBU/noUBMTKaZDqSS20V2iQIoIrxgTiWgH\nUkDA2omfqIR2oPRjKTv+zjlv5o2EP+xuGPYk797z7jn3nd/9nXfvfXcGmE0SHXx6NsHpwLJmWaXj\nbpaos4+uXrw6S7h5NwwvM/bupg/+fvbRFUOf+8Hz8l4ISk7jvZo/2Lbu0BUeeVaH2QuEFy3EpZ/d\n1q4xMgJ8ZRCIei0mbmopXagjQ3mJEh6qYuApF0PLXQWDz7ONoD923HHQRoUSH9OqO4XwQiEQDF6I\n5NBC4IFhN/r8vVW5AHJp0g26TsciDXari7PxAvY1a+06mY/XIrVILrp4zAfkXOkr3agewC0S5trE\nRLiSLoWqIU/AYPUfAQG9t9n0pNGgzO/iHii8MKoAidU34cf4MJJFx8s034FcALlU6Qpd/RiXYM4K\nF+uxtvLPT6Fd9xTRyNYTDbRf74/4yAzf+1zejizDRWpOoacSPRLKJ/crGPC1FtB/fv01VxoVhWPL\nyPuMyH+88CIqgazCPDyNGHorKGVucCGXcqmuXaLrFAkW4e72BDbXX6feqvsqeHRzEVPtBDqnzVVg\nXSmcI0iVYVc9MQWESgpGeANOQRnfZB0qyR3m3JPXuiuF8CKBCETqU4s/Quj2RC1TB+RSLlmhS3Q5\nD0uwzzCBMvf/spILul+vdfG9viIOaQK3bFtV6WkeFt/uSLSBniqQvX6YyOJ85RcUFYwmUEA7krMF\nxeTpW1/ChuZYd1BJFOFFaibwEOk6/o08nA39yJw07MqlXAZ0zRk6WZwxZ+g0YN7IEqoxF9G7+gcH\n2Sa28NBHcclzg4McTzTHe6wG1rmr62aJfTy6BInBZeVWbzF3iHNeQe7ymsDbcAJu5p3VBHpbqhyd\n1AT+bKqnLo5tcfqHi3LDusKKoQONm/uIGBlKMfq+5kODDkZ9nJUDdXRglCeaSAILwOacnBCozVvH\n2UgwwpuCTrsc0jok99Xv9rv4VYuyjYt+44dSOjSoIRHKKC2i1FHRG4/qroUR7vPaiiK8iH4TEm/w\nYJW8bj3w4OGKU4l7cimXAV1X4Hjhx3kMmyrIIV1GYpeLZLPJ7VJtzrVcT8Z5X+M4PXnwncCTzQPs\nI5ZUs1nARsQ5cu2tZvFqi3zEJCuSwHCOKhcIv15bwS/iXkSX0MQ8j7ZOiZUdAuccqCZ2sGLoQMP8\neqjERiJRjL4vfVSzwahPCHixEyM7+RLlElognDLXcAywMeQpGJ2B8uUlH4FsvOcTUndKi7IaVtQt\nlNBhmsUXYiiEpxZxNPTGo7KnhU+Z8dqKIbyIfhOcQ8jWz4ne4YarlzxEwo/IJVwGdCUnEC7ROT2B\ntBevIlFyntrrIsSk1swWy+MZcHam6XWRxyLyOHD5bjqoJbb8AuCTZRyG9jYzvRB/hrKTytoy32Q9\n4W0Bd+S+OvyaK9YUPn3uEZ2B6YJ0MjlBqu3AWVK/AvyJFUMHGi5Hr6ehDKPvSx/VdDDm83XgBwgw\nygN9ZM5bJIdLaB2TSObYHD2kYCyBBM3zjTbultrE4ajRpsxFasIPRTpM0/hKTEAUp3JqwtAbj8qe\nFD5lxisffbNw5iovvONSMIXs/SW+lfP5pvZWnCm5lMs2Xb01Xf+RLSHcCHtwxpgOF2FOl6LZxgV0\nHuCikHjQo37v1Uwua6hF8or/1J2D8HuLORBuJpurPO2duPR0ZPKZ4npsdls118yGg4ynCUwVg07X\ni/p7/3P5ILDC1dCBluSrSGEow2i+Ak+1rAzGfJYBx+u7Zxi1m1/crV+hsSpZ+iv+vqCYmFQw8uIr\n6HWyAyYm8XDQKV6g3qIsNYb0WxpK6TBN4ysxAVHqCEVvPCp7UviUCfAOEV7klgm8Ctl/eYiU7wPO\n76kkGnIpl226sjVkDtG5j25TmaOcq5pAtvwNZvsddUqmwHgxj1r+aj+BakmzhcIl1O8t5g5Zj+3u\nfS73GhcX33izeyZXq1YdycVrvYXeIuITtHf00QTywLOvzCXjKAHVJXSgoZcdKAxlGM1X4KmmgzGf\nfTtxjiVQMWo3vxjmOfA+N+PFS9ErbtzZU0k+pGDwDTnobHfxGNdtNib2B500gRaSC1CDCdRQmkDT\nNL5PWYsodfTRK481eaIUPmXGazuM8JLg2kUgr+FruBJpN1bBs2E3nJNLuWzT1ZdD5gC7pvjSTmLv\n5DnUOQMpHsx2cM1wnbe9vH4S82Q7Y4Ze3noyoJb0Nv3WiOVgvdXsDA1X6C+yZuQlfK6C5K4XuR80\nER7It2sMDLiZwdOQOPWdAqEGIgl03mACC1Qe4QysauhAy35xaIkhUYy+L+GZpoMxn3TzgSo6MHYg\nC/cvFGRbBlzu4/ujQ8srAgYbN+xW0LibC6g0egThiyawRRnbehrGm9ChmsU3ygKixNEfh/Ko7Glh\nhCuvc1cyuAl5ie5QIElinDN0Npz+kQoG+/nxxku4bNO1OY/om+wW3o/YAYSal1HXBHIqqs05WNXT\nCDPrFGKebGdM4E4e6czS68oZas232E97q/nVSiLPhumKJDBBVGvlISuA6woaOtD6rpLzm4RSjOYr\nPqbpYMwHSw+66MA4M2SawBZlxLYvb7wJHapp/DONsjZR6mjjANcDY08pNMp0cJn1+D5NxyxtNOvw\n1Ulsufjtup/APs5USa7TdHGGC2c/k4yYx3VLEsgZOeZbgN/yNsydQnuLOfEo5hTFaZryPwmMFeIr\nqho60PoaCD2kSDoSKD5+2mUw5pM8+bonBIWPcYbI3p3AO3mgEt6EDtUsgUZZmyia6ajolUdz1z5K\nmQ5uU9nhonbs0pe3JRTJlWsOJWsY59eqzsDzuC+q7TC/VuqIlIAbEOMHR9ESGDrCD0+xAEsrLH7p\nam81px9beT+bpi26hHIG6hKKreeuKGvoQMuOITmloRSjo77iY5oORn3wUyTfdgnExzhDZLaEtijj\nGYPrGXnj80mHaBZfiekgSh0VvfJo7tLHCNfBvbL1QsF5zJLNIS4AKOlGiid43mgC+c6YjTvPeFWO\nYfxlIuZhLpjAdAnhKVnVx6svcJGr/xDYVdbeauY2P305a3T0jtHRL4B74K1le8ze4zQ0b3wtVUJy\n0kIpfvH9svq0enEw4gO+j2vLAcaZIAuPjn7pkdHRXEDZKoWXbhgdwqChVmI6iOJPmip7XeHR2NM+\nSpnxyq+faQkfmNRvVs7xGpdNvGwJdPhAs/GTivOsj8eNPXv2/qpwy549b+8mOaEG1EKwS12uGbuq\n2lvN2fy0sLQ7yQyUVWfcfyevsdBs9DV+YSanNJRhFN8PKbxWr1hNfWTlT1cDjDNEpjOwTVk8B765\n5M3oEM1QKzEdRPmOuIarWk0/GEPyfoUaSpkN7k0Z8zQkMmG/1yZ38K+CVEVO2joDo3yVzXYr90C+\nxjl5eMqT8nGEmN2S/KtxBnLAd9iCDa3evEnlgYJ4TlM0gRfJL8qU7fyrgBVDBxp/ugxNsPFxH6Pv\nSx/VdDDms4NbTSUXYJwZMk1gizL5pewCDaV0qGbxjbKAKHH00QuP6q6FT7gM7ufkmyM6dvkuT6s8\nhUcuc65C5jIkC5bAOBMIsaE37zzKFSonj856Uh5FIkc8ZtkIIt+G+IFWb5q5Yc7huzBt0QT2lHme\n7ivhPHd+lU9i6EDDr7GgzEaeERWj+qqPajoY8/mHi638pa+FcYbINIEakpQlvj24eExDKR0WVOMb\nZW2i1NFHLzyquxY+4To4OO1jBEf2/8vclc8Dd/HXwEV14NL+JcDiJ++qyockoDYMLC4Cm4Sv8K53\nqsCy5svYuHwh78WS6H/ORebEdm81D19Iv+mLJtAZWlVETw3hE//AJ0noQGMjUWsoxai+6mOaDkZ9\nos8N1DsxzgyZJdCnLMbfh8eMN6VDg/rxhbI2UeZo6JVHdddCCVfgkZO2TZ+w2dZTEzjbQAkeS+Bs\nRDa7MP17dsEJ0GSqgT47tP8CCynLwhqWBXMAAAAASUVORK5CYII=\n", + "text/latex": [ + "$$0.981706744565384 e^{- 3.42525091 r^{2}} + 0.949464004096853 e^{- 0.62391373 r^{2}} + 0.295906452975584 e^{- 0.1688554 r^{2}}$$" + ], + "text/plain": [ + " 2 2 \n", + " -3.42525091⋅r -0.62391373⋅r \n", + "0.981706744565384⋅ℯ + 0.949464004096853⋅ℯ + 0.2959\n", + "\n", + " 2\n", + " -0.1688554⋅r \n", + "06452975584⋅ℯ " + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# STO-3G for Hydrogen\n", + "h_alpha = [3.42525091, 0.62391373, 0.16885540]\n", + "h_coeff = [0.15432897, 0.53532814, 0.44463454]\n", + "nbasis = len(h_coeff)\n", + "\n", + "cg_unroll = cg_sym.doit()\n", + "for idx in range(nbasis):\n", + " cg_unroll = cg_unroll.subs(c[idx+1], h_coeff[idx])\n", + " cg_unroll = cg_unroll.subs(alpha2[idx+1],h_alpha[idx])\n", + " #cg_unroll = cg_unroll.subs(norm2[idx+1],norm_sym.subs({i:0,j:0,k:0,alpha:h_alpha[idx]}))\n", + " print(h_coeff[idx],h_alpha[idx],radial_norm.subs(alpha,h_alpha[idx]).evalf())\n", + " cg_unroll = cg_unroll.subs(norm2[idx+1],radial_norm.subs({alpha:h_alpha[idx]}))\n", + "\n", + "cg_unroll.evalf()" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKAAAAAPBAMAAACRq9klAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEJmJdjLNVN0iZu+7\nq0QgoRR7AAAACXBIWXMAAA7EAAAOxAGVKw4bAAACo0lEQVQ4Ea2TTWsTURiFz8w0SZPJx2BLFSok\ntujKlkCLC1EyWBWkFEMX7sQUlYIoxoVLacCFCBVCd4rQIEWQLpouFATFEXHhQs2i3bSG5h+U+NHa\nz3je90Z/gYE8HM7JPdx73xvgwPBRyOeaN94wwNBIHjPHP4obysDqP5M3ihwXUzEhysTlttk9Mkg1\niUNqvGv9BBQ9Dbtg1ZBrMO0MEC5aD40ip/lVPG2KYtxbib5umxkkiohV4RTEuLrEXgPEg4QHt0p3\nNMA8cMIooOMNleDwgBYyXgW+GTNSQbSAeA0RzWp0IYhtEW4JiV9A9HmAz8CUpwp4cpOZokMWSbwB\n5Dw1nQBWCakakpvMtEsRl52Fm1rohAPsAetFiAKyUqjQQprWLnCxoWZy13MqSGeQ5OUBK3McjiB1\nfVjuFp3cwwuu+M5CXxQ3X2ahAloo5ivusGLMte2TwGwWoR9cjveYbSjSN+AWaaxnYfnhIMr0SlYU\n0A0WKrRQzRxw2zem3eLo/xUC8RIE6SbsRappwMHfQlE8lxQqtFDNsB/hDtV8dmmngXS2fWTAlmnY\nW6kSYvsULLgLHpk7XPdFIZJnoYKj5oWoiblzuaKasRoWFpHKICJDSRTg7CvcAmLbwBhgZaWGdzh1\nR1UXWKjQQhNz7ZqnpushtClPMyYvgDV2U8Epc4eRDB5H6vW1tz6PvnBQ1b16fWdJoYUm5tpbUDNN\nuYKOKuwClc2xFhR8lnZV/l4XaLsBRoFjRpF80wY6ZYnnyyG5KiZuGfgEPEBv3tpANMNIgQ/oLUbv\nDw2U+LNUgM6i9dIokq/OwBQyPu/1VIyZnEDMB7rOLgOPgJnTR9pw+pYRbrVaJY75617FGh7LGwWc\naq0Y2Ku/qSR2+r5InySX+wdF/t/PH/PtA7AHgkKuAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$0.55624044414948$$" + ], + "text/plain": [ + "0.556240444149480" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Evaluate at a concrete value for r\n", + "cg_unroll.subs(r,1.3).evalf()" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAALkAAAAPBAMAAABKEHMHAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMpmJdlQiZu+7\nq0TEZSulAAAACXBIWXMAAA7EAAAOxAGVKw4bAAACw0lEQVQ4EbWUy2sTURTGvzt5NI8mjC7sRmgL\nLosERUFQDIrdiUUwIH0FqSlIxSpUECnOPyDNpouqi2zc1EVUuqjYRUA3otIgFKFSnJ0LFzWtL/pw\n/M65Yxz/AA/Jd+85353fzJy5M0AYpjRV1GlyfAxQQS9gy7NDL4GumyOy4AyQrkx4cK7HR3B6pVIp\nx0+UKhVra26PtghFAo5njsnUHMR5K3jQCsumiWs+BpDzuOA4MIvUFrJBMIl6EATNDGXX2ppHEEKU\neARcktEp442V/cOkaznnIvs4VUO6B+hYAu542IRz9RZwEcghUQWa1tb8L0KAGm+BGZezOpdaQQfp\nWs4WkPuabMAUgHt9wGffbJDKKANjSLKTxaTamkcQCGMXWPc4fy25itK17LRIz++4yRp5pEM6o3Qg\nPyn5J4S25BGEeAzzjXRZuNE/5VsRersca+Hj1hUgU1W6w4cwr9ugUw9vwNoA8whCTW6D78C5Ms+y\nUcOiClvcapexXkYi4B7qQh+P6D9CjJv+xSnPyMbUYG3JowgxGW164GJwj4j7D5075eGpbZ+NFjqS\nC6IveDVfZNLNv9qSmwhi72GJA4bXrp3Z5DP1VbQzYTnRQKaJ+pNU0dLx3CVwtIqOHo64wZaJzftl\nHkGIKcG+z3gcn5JeUxH6n/I0kHUR/7kPQp8Dlr13kJcg1uAx8g6oDc0jCDEluKAu18MHXvdVlG7L\nqQHMyd2v3l1b234P3vqyXMGoi+4my4Y7S20OzCMImhq3wXeUwU0/aEXptsy3/1a2CrzigiX9nawO\nAIe4EzggzhsPbckjCCFKxDyzgO4edJbNohWlazl9tDJcyPciM8mFO8A8Uj/4Ocg8I0noKdJDW/II\nQsgSZny6iFgTuDxctJL4sLlqyw4/JAWcLclXbCJYRf7CkI90acUF7nuspRcp1ta8jWD5P8ZvGoQg\nxMf+bV4AAAAASUVORK5CYII=\n", + "text/latex": [ + "$$-0.661028435778766$$" + ], + "text/plain": [ + "-0.661028435778766" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Can check derivatives as well\n", + "d_cg = diff(cg_unroll, r)\n", + "d_cg.subs(r, 1.3).evalf()" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 23, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgUAAAFkCAYAAACw3EhvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAAPYQAAD2EBqD+naQAAIABJREFUeJzs3Xl8TPf6wPHPySZiSWissSSxxb7EVkSEqK2l6C1KbbXe\nLqTbrdLWbRU/1VJb61ZvSbloqV0tIUi0tcROEGusIdZIIpKZ8/vjSCr7zGRmMkme9+s1L3LOd3mm\nyjxzvpuiqipCCCGEEHb5HYAQQgghbIMkBUIIIYQAJCkQQgghxFOSFAghhBACkKRACCGEEE9JUiCE\nEEIIQJICIYQQQjwlSYEQQgghAEkKhBBCCPGUJAVCCCGEAIxMChRFmaAoyn5FUR4qihKjKMoaRVFq\n51LHX1EUfYaXTlGU8nkLXQghhBDmZOyTAj9gLtAKCAQcgW2KohTPpZ4K1AIqPn1VUlX1lpF9CyGE\nEMKClLwciKQoijtwC2ivqmp4NmX8gZ1AGVVVH5rcmRBCCCEsKq9zCtzQngLczaWcAhxRFOW6oijb\nFEVpk8d+hRBCCGFmJj8pUBRFATYApVRV9c+hXG3AHzgIFANGAq8DLVVVPZJNneeALsAl4LFJAQoh\nhBBFkzPgCWxVVfWOMRXzkhR8h/bB3VZV1RtG1t0FXFZVdUg2918DlpkUmBBCCCEABqqq+j9jKjiY\n0ouiKPOA7oCfsQnBU/uBtjncvwSwdOlS6tata0LzBUdQUBCzZs3K7zAsTt5n4SLvs3ApKu8TisZ7\njYyMZNCgQfD0s9QYRicFTxOCXoC/qqrRxtZ/qgmQUzLxGKBu3bo0a9bMxC4KBldX10L/HkHeZ2Ej\n77NwKSrvE4rWe8WE4XejkgJFURYAA4CeQLyiKBWe3nqgqurjp2WmAh6pQwOKoowDLgIn0cY5RgIB\nQGdjgxVCCCGE5Rj7pGAM2mqDXRmuDwOCn/6+ElD1mXtOwNdAZSABOAZ0UlV1j7HBCiGEEMJyjEoK\nVFXNdQmjqqrDMvz8FfCVkXEJIYQQwsrk7IN8NmDAgPwOwSrkfRYu8j4Ll6LyPqFovVdT5GlHQ0tR\nFKUZEBEREVGUJoQIIYQQeXbo0CF8fX0BfFVVPWRMXZOWJAohhLAt0dHRxMbG5ncYwgrc3d2pVq2a\nRdqWpEAIIQq46Oho6tatS0JCQn6HIqzAxcWFyMhIiyQGkhQIIUQBFxsbS0JCQpHY8K2oS92YKDY2\nVpICIYQQ2SsKG74Jy5LVB0IIIYQAJCkQQgghxFOSFAghhBACkKRACCGEEE9JUiCEEEIIQJICIYQQ\nNmzJkiXY2dlx6FDWG/N16NCBRo0aGdxeWFgYr776KlWqVKFYsWK4ubnRunVrvvjiC27dupWuz9xe\n3t7e6dreu3cvvXv3pmLFijg7O+Pl5cWYMWO4cuWK6f8BrEyWJAohhLBpiqKYdC+jTz/9lClTplCj\nRg2GDRuGt7c3jx8/JiIigm+++Ybg4GCioqJo3749S5cuTVf3jTfeoFWrVowaNSrtWsmSJdN+P3fu\nXMaPH0+NGjV45513qFSpEpGRkSxatIiVK1fy+++/07p1ayPedf6QpEAIIUSht3LlSqZMmUL//v0J\nDg7GwSH9x9+sWbOYNWsWAF5eXnh5eaW7P3r0aLy9vXnttdcytb13716CgoJo3749v//+O87Ozmn3\nxo4dS5s2bXjllVc4efIkrq6uFnh35iPDB0IIIQq9Tz/9lHLlyrFo0aJMCQFAqVKl+PTTT01q+4sv\nvsDOzo4lS5akSwhASzBmzJjB9evXWbhwoUntW5MkBUIIIWzegwcPuHPnTrpXbGwsycnJudaNiooi\nKiqK3r174+LiYta4EhMT2blzJ35+ftluO9yvXz+KFSvGxo0bzdq3JcjwgRBCCJumqiqdOnXK9n6D\nBg1yrH/69GkA6tevn+nenTt30v3s5uaGvb29wbFFRUWRkpJC48aNsy3j5OREnTp1iIyMNLjd/CJJ\ngRBCFCEJCfD0M9KifHzAXF/KFUVhwYIF1KpVK9O9d999F71en2P9hw8fAuknBoL29KFcuXIoioKq\nqgAcPHjQqPMj4uLiAG34ISelSpVKi8OWSVIghBBFyOnT4Otr+X4iIsCcZzO1aNEiyw/rMmXKZPq2\nn1HqB/ajR4/SXS9ZsiQhISEAbN26lZkzZxodV2rbqclBduLi4nJNHGyBJAVCCFGE+PhoH9jW6MdW\n+DwN5sSJE+mu29vb07FjRwCT9xKoWbMmDg4OHDt2LNsyT5484cyZM7Ro0cKkPqxJkgIhhChCXFzM\n+w2+IKhduza1atVi7dq1zJ49m+LFi5utbRcXFwICAggNDeXKlStUrVo1U5mVK1eSlJTESy+9ZLZ+\nLUVWHwghhCj0Jk+ezO3btxkxYgQpKSmZ7uc2LyEnkyZNQq/XM3ToUB4/fpzu3sWLF/nwww/x8PBI\nt/GRrZInBUIIIWxa6iTAvBgwYAAnTpxg+vTp7N+/n/79++Pl5UV8fDwnTpxg+fLllC5dmjJlyhjd\ntp+fHzNnzuS9996jUaNGDB06NN2OhgBr1qyx+Y2LQJICIYQQNi63rYwN3er4yy+/pGvXrsybN4+f\nfvqJ2NhYihcvTu3atfnggw8YPXo05cuXz7aPnPoZP348LVq04Ouvv+bbb7/lwYMHVKpUiX79+vHx\nxx9nOaxgiyQpEEIIYbOGDBnCkCFDsr0fGhpqVHt+fn74+fkZHYchywnbtm1L27ZtjW7blsicAiGE\nEEIAkhQIIYQQ4ilJCoQQRjl2DDZvBgO2nBdCFDCSFAghcqXTwapV4O8PjRtDjx7g5QVTpkBsbH5H\nJ4QwF0kKhBC5+uc/4R//0H7/yy/ajnjdusHUqdC8Ody+nb/xCSHMQ5ICIUSOfvgB/vMf7dfdu7Xk\noFkz7efISO2Anf79IYv9YIQQBYwkBUKIbP31F7z1FowZAyNGZL5fvbr25GD3bpgwwfrxCSHMS5IC\nIUSWbt2Cvn214YFvv82+XIcOMHOm9vrlF6uFJ4SwAEkKhBBZmjpVGxr49Vdwcsq57LhxWgIRFAQZ\ntn4XQhQgkhQIITK5eRMWLtQ+5CtXzr28omhJxM2b8N//Wj4+IYRlSFIghMhk5kzt6cA77xhep3Zt\nbcLhtGmQlGS52IQQliNJgRAindu34bvvtCEBNzfj6k6cCNeuwZIllolNCGFZkhQIIdL55huws4Px\n442vW6+etmRx2jTZ8VCYx5IlS7Czs+PQoUNZ3u/QoQONGjUyqK0NGzbQoUMHKlSoQIkSJahRowb9\n+vVj27ZtAAQEBGBnZ5fr6/PPP09rMyUlhTlz5tCyZUtKly5NqVKlaNmyJXPnziWlAK7TlVMShRBp\n7t6FefO0ZYhly5rWxqRJ0KgRLF0Kw4aZNz5RNOV0ZLGhxybPnDmTDz/8kA4dOvDxxx/j4uLCuXPn\nCAkJYcWKFbzwwgtMmjSJkSNHptU5cOAAc+bMYeLEifj4+KRdT01CEhIS6N69O2FhYbz44osMGzYM\nOzs7tmzZwrhx41izZg2bNm2iePHiJr5z65OkQAiRZulSbT5AUJDpbTRsqG2D/P33khQI26DT6Zgy\nZQpdunTh999/z3Q/9ule3Z06dUp3vVixYsyZM4fAwEDat2+fqV5QUBBhYWHMmzePsWPHpl0fPXo0\n3333HW+++Sbvv/8+8+fPN/M7shwZPhBCpAkOhu7doXz5vLUzdCjs3w+nT5slLCHyJDY2locPH9Km\nTZss77u7uxvd5rVr1/jvf/9Lp06d0iUEqcaOHUtAQACLFi3i+vXrRrefXyQpEEIAcOqUdqbB4MF5\nb+vFF7VJij//nPe2hAB48OABd+7cSfeKjY0l2YDJK+XLl6d48eJs2LCBe/fumSWe33//Hb1ez+uv\nv55tmcGDB5OSksKWLVvM0qc1SFIghAC0D/AyZbRH/3nl7Az9+mlt6vV5b08Ubaqq0qlTJ8qVK5fu\nVb58ef74449c6yuKwgcffEBERATVqlWjR48eTJs2jcOHD5sc06lTpwBo3LhxtmUaN26MqqpERkaa\n3I+1yZwCIQQ6nTafoH9/KFbMPG0OHqxtgLR7NwQEmKdNkXcJyQmcjrX8uI6Puw8uji5maUtRFBYs\nWECtWrUy3Xv33XfRG5B5Tp48mbp167JgwQK2bdvGli1bmDhxIk2bNmXZsmXpJhIaIi4uDoBSpUpl\nWyb13sOHD41qOz9JUiCEYNcuuHrVPEMHqZ5/HmrU0OYpSFJgO07Hnsb3P74W7ydiVATNKjUzW3st\nWrSgWbPM7ZUpU4Y7d+4Y1Ea/fv3o168fjx49Yt++fSxevJhly5bRs2dPTpw4gVNu+3k/I/UDPzU5\nyIohiYOtkaRACEFwMNSqBa1ama9NRdGSjK++gvnzwcU8XxpFHvm4+xAxKsIq/diqkiVL0qlTJzp1\n6oSDgwPBwcHs27cPPz8/g9uoW7cuqqpy7NixbPdJOHr0KAD16tUzS9zWIEmBEEVcfDysXg0ffaR9\nkJvToEHw2WewZg0MHGjetoVpXBxdzPoNvqBr3rw5wcHB3Lhxw6h63bp1w97enp9//plBgwZlWSY4\nOBhHR0e6du1qjlCtQiYaClHEbdumJQb9+5u/bW9vaN0afvvN/G0LYajExET++uuvLO9t3rwZgDp1\n6hjVZpUqVRg2bBghISF8//33me5///33hIaGMmLECCobcqqYjZAnBUIUcRs2aNsT16xpmfZfekk7\nQfHxY21VghDGUlU1T/UTEhJo06YNrVu3pmvXrlStWpX79++zdu1awsPD6d27d7arCHLqe9asWZw5\nc4Y333yTLVu2pD0R2LJlC+vXrycgIICZM2fmKXZrkycFQhRhOh1s3Kh9cFtKz57ak4hduyzXhyjc\nctvKOLf7bm5uLFq0iEqVKrF48WLefPNNPv30U+Lj45k5cyYrVqwwqe0SJUqwY8cOZs2axfXr1/nw\nww/54IMPuHbtGnPmzGHbtm0FaotjACWvGZglKIrSDIiIiIjIcrapEMI8/vwT2rSB8HBo29Yyfaiq\nNozQvbs24VCY36FDh/D19UX+zSz8DPmzTi0D+KqqmvVJUtkw6kmBoigTFEXZryjKQ0VRYhRFWaMo\nSm0D6nVQFCVCUZTHiqKcVRRliDH9CiEsY/16cHfXxv0tRVG0JxEbNmgJghDCdhk7fOAHzAVaAYGA\nI7BNUZRsn48oiuIJbAR2AI2Bb4FFiqJ0NiFeIYQZbdig7WBob2/Zfnr2hCtX4OkKLSGEjTJqoqGq\nqt2f/VlRlKHALcAXCM+m2ljggqqqHz79+YyiKO2AIGC7UdEKIczmwgU4eRKeORreYtq3h9KltSSk\nSRPL9yeEME1eJxq6ASpwN4cyrYGQDNe2As/nsW8hRB5s2ABOTvDCC5bvy8kJunTR+hRC2C6TkwJF\nm5I5GwhXVfVUDkUrAjEZrsUApRVFMdMu60IIY61fDx07QsmSxtVLSE7gepzxR8H27AkHDkABOkVW\niCInL/sULADqARaaswxBQUG4urqmuzZgwAAGDBhgqS6FKBIePoQ9e+Dbb42v+/nuz1l2fBnR46Nz\nXQr2rG7dwM4ONm+GESOM71cIkdny5ctZvnx5umsPHjwwuT2TkgJFUeYB3QE/VVVz2xvyJlAhw7UK\nwENVVZNyqjhr1ixZXiOEBezeDSkp2iN9Y+24uIOrD68SdTeK2s/luvgozXPPQYsWEBIiSYEQ5pLV\nF+VnliQazejhg6cJQS8gQFXVaAOq/Al0ynDthafXhRD5ICQEvLy0UwyNEZcUx+Eb2hn0YZfDjO63\nc2fYsQMMOOlWCJEPjN2nYAEwEHgNiFcUpcLTl/MzZaYqirLkmWrfA96Kovyfoih1FEX5J/AK8I0Z\n4hdCmCAkBAIDja/359U/0ak6niv+HHui9xhdPzAQYmPh2DHj+xZCWJ6xTwrGAKWBXcD1Z16vPlOm\nElA19QdVVS8BPdD2NTiCthTxDVVVM65IEEJYwbVrcOqUaUlB2OUw3F3cGdhwoElPClq31o5Q3i6L\nkYWwSUYlBaqq2qmqap/FK/iZMsNUVe2Yod4eVVV9VVUtrqpqLVVVfzbXGxBCGGfHDu3Xjh1zLpeV\nsOgw/Kr50b56ey7ev8jVh1eNql+smLZnQYh8JRDCJsmBSEIUMSEh0LSptr2xMZJSkvjr6l+0r94e\nv+p+gGnzCgIDISxMOzVRCGFbJCkQoghRVdPnExy8fpAkXRJ+1fwoX6I8Pu4+7Lls2ryCxETtMCYh\nhG2RpECIIiQyEm7cMC0p2HN5D6WcStG4onbuvF81P8KijX9S0LAhlCsnQwhC2KK8bF4khChgQkK0\nLYfbtTO+blh0GG2qtsHBTvtno3319vxw6AfuJNzhOZfnDG7Hzg46ddJi+fJL4+MQwhri4+OZP38+\ncXFxJCYmcu7cOaZMmUKDBg1Mam/kyJEMHz6c559/3iLtm4s8KRCiCAkJgbZttRUAxtDpdey9spf2\n1dunXfOrps0rCI/O7iy07AUGwsGDcO+e0VWFsIoJEyawe/duvvjiC2bOnImnpycBAQE8fPjQ6LZ2\n7drFjz/+SHJyskXaNydJCoQoIlJSYNcu7Vu6sY7FHONh0sO0RACgult1qrlWM3legV6vxSOELVJV\nlZiYv4/t8fHx4e7du5w5c8aodpKSkti5c2emLcHN1b65SVIgRBFx+DDExUFAgPF191zeg5O9Ey08\nWqS73r56e5M2MapeHTw9te2WhbBFc+fO5eDBg2k/R0VFUbJkSerWrWtUO/PmzeOtt95CVVWLtG9u\nMqdAiCJizx4oXhyaNze+builUJ6v8jzODs7prvtX9+d/x//Hg8cPcHV2zaZ21vz9tZiElSUkwOnT\nlu/Hx8f4caocxMfHM336dJycnEhMTOTOnTvMnTsXJycnEhIScDFjXxk9ePCAlStXEhwcTEkjjhU9\nceIEFSpUoHz58hZp3xIkKRCiiNi9G9q00SYaGkOv6tlzeQ/jWo3LdC/AMwC9qicsOowXa79oVLv+\n/hAcDPfvg5ubcTGJPDh9Gkw8LMcoERFgpgPt4uLiCAgI4JNPPqFXr14AjB49mkmTJjFjxgxmz57N\nxx9/bJa+npWUlMT333/Prl27ePfdd9P6NoSqqixbtoxp06ZZpH1LkaRAiCJAp9M2DAoKMr7u0ZtH\nuff4HgFemccdvMt4U7V0VUIvhhqdFLRvr+2bEB4OLxpXVeSFj4/2gW2Nfsxk/PjxeHp6pvvQ7Nix\nI+PHj2f69OnY29unK6/X6+nbty9JSdpBvBkf3aeO76uqipubW6ajh1MVK1aMcePGMW7cOPr378/L\nL7/M2rVrDYr5xx9/5I033sixTF7atxRJCoQoAo4f176R+/sbXzf0UijODs608miV6Z6iKHTw7MCu\ny7uMbtfbGzw8tCcYkhRYkYuL2b7BW8PNmzf5+eef2bhxY7rr7u7uxMTEsHDhQvr165funp2dHWvW\nrDFrHCNHjqRz584sWbKEIUOG5Fj26tWrJCYmUrNmTYu0b0ky0VCIImD3bu3cgVaZP9dzFXoplDZV\n21DMoViW9wM8Azh84zD3Eo1bX6goWpIikw1FTvbv349er8fPzy/ddQcH7TttbGwsnp6eZu0zJiYG\nDw8Pvvjii7Rr1apVS4snN1u2bGH//v0MHz6c4cOHM3DgQACmT5/O+++/T0xMDJUrVza5fUuSJwVC\nFAF79mgJgbNz7mWfpdPr2HN5Dx+0+SDbMgFeAaio7Lm8h14+xo2J+vvDypXaqohSpYyLTRQNOp0O\nNzc3ihcvnu66nZ0diqLw9ttvZ6qTcfggO9kNH9y8eZMbN25w//79tGuxsbEAeHt75xrziBEjGDFi\nRNrPly9fZvny5UyYMAE/Pz+OHj3KzZs3TW7fkiQpEKKQU1UtKRgzxvi6h28e5mHSQwI8s1/H6Onm\nSXXX6uy6tMvopKB9e22+wx9/QJcuxscnCj9/f39UVSU2Nhb3p6d4XblyJW3sPSUlhePHj9OwYcO0\nOnkdPmjUqBGdO3fmnXfeSbu2atUqqlSpwrBhwwBYuHAh06ZNY9++fVSoUCHH9lJSUgAtwTG0/fwi\nSYEQhdypUxAba+J8gouhuDi6ZNqfIKMArwBCL4Ua3X6dOlChgjaEIEmByErZsmX59ddfGTduHPXq\n1UNVVSpWrMiMGTPQ6XRMmDCBwMDAdElBXimKwrJly5gyZQp6vZ7ExETi4uIIDw+nbNmygPaUISkp\nCb1en2NbX375JStXrkRRFEaNGkWXLl2YO3cuS5cu5csvv8y2/fyiZJyVaQsURWkGRERERNCsAE2I\nEcIWLVgA48ZpEw1LlDCubo//9SBZl8y217flWG7JkSUMXTeU2A9ijToHAeDVV+HaNdi717jYxN8O\nHTqEr68v8m9m4WfIn3VqGcBXVdVDxrQvEw2FKOR279Y2LDI2IUjRpxB2OSzHoYNUqcsVTdnyuH17\nOHBA21NHCJG/JCkQohBTVW1/AlOGDiKuRxD3JC7L/QkyquZaDe8y3iYNIfj7Q3Iy/PWX8TEKIcxL\nkgIhCrGLF+HGDdOOSt55cSelnErhW8mw3e86enZk58WdRvdTv762o2G48YctCiHMTJICIQqx1A/a\nNm2MrxtyMQR/T38c7R0NKh/oHcjJ2ye5EXfDqH7s7LTjnCUpECL/SVIgRCG2dy/UqwfGTmhOSE4g\nPDqcQK9Ag+t09OoIwI6LO4zrDO1Jxp9/asc7CyHyjyQFQhRi4eGmDR3sjd7LE90TAr0NTwrKlShH\nk4pNCLkQYnR/7drBo0dw7JjRVYUQZiRJgRCF1N272h4FpiQFIRdCqFiyIvXK1TOqXqBXINsvbM90\nAE1umjfXTm+UIQQh8pckBUIUUn/8of3atq3xdUMuhhDoHZh2mpyhAr0DuR53ndOxp42q5+wMLVtK\nUiBEfpOkQIhCau9eqFQJvLyMqxebEMvhG4eNmk+Qql21djjZO5k8hBAeri2jFELkD0kKhCikwsO1\npwRGftkn9GIoKiqdvDsZ3WcJpxK0qdqGkIumJQU3bmjLKIUQ+UOSAiEKoaQkbZdAU+cT+Lj7UKV0\nFZP6DvQKJPRiKCl645YSpC6bDAszqVshhBlIUiBEIRQRoSUGJiUFF0NMGjpIFegdSNyTOA5cO2BU\nvTJloEEDmVcgRH6SpECIQig8XDvroHFj4+pdvHeRC/cuGLUUMSPfyr64FnNl+4XtRtdNnVcghMgf\ncnSyEIXQ3r3QujU4GPk3fNv5bdgr9nTw7GBy3w52DnT06sj2C9v51P9To+q2awfffw+3b0O5ciaH\nIIRZTJs2jTt37lCyZEkuXrzIvHnzKFWqlMH1ExMTmTx5MklJSZQtW5Zy5coxduxYABISEpgxYwZ3\n797lyJEjeHl5MWPGDCpUqGCpt2MQeVIgRCGjqlpSYMpSxC3nt/B81edxdXbNUwxdanThzyt/8uDx\nA6PqpQ53pC6nFCK/zJ8/nz179jBz5kwmT55M3bp1ef311w2ur9fr6dOnD/Xr12f27Nl07tyZoKAg\nTp48CcAXX3zB6NGjmTNnDnv27CEmJoaOHTuSnJxsqbdkEEkKhChkzpyBO3eMn0+QrEtmx4UddK3R\nNc8xdKnZBZ2qM3rL4+rVoWpVGUIQ+W/GjBkMGTIk7efBgwezfv16zp07Z1D9n376iYSEBAYPHgxA\n1apVGTBgANWrVycpKYl58+bx448/ppV/7733iIyMZP369eZ9I0aSpECIQiY8XDtkqHVr4+r9efVP\n4p7E0aVmlzzH4OnmSZ3n6rDl3Baj68q8ApHfoqKiuHLlCvXq/b2jZ+XKlXF1dSU01LDjwefNm0e3\nbt3Sfq5SpQo//fQTJUuWRKfT4e7uTnx8fNr96tWrA3D+/HkzvQvTyJwCIQqZ8HBtgqERQ58AbDm3\nBXcXd5pVamaWOLrU6MLaM2tRVdWonRHbtYNff4WEBHBxMUsoohCIj49n+vTpODk5kZiYyJ07d5g7\ndy5OTk4kJCTgYsb/Wc6fP4+iKJQuXTrd9VKlShEdHZ1r/djYWI4ePcrIkSOZOXMm8fHxnDlzhi++\n+IIaNWrg4uLCxQwbcly6dAkAL2N3GzMzSQqEKGT27oVnvqAYbOv5rbxQ4wXsFPM8QOxasytz9s/h\ndOxp6para3C9du200xL374cOHcwSinhGgk7H6YQEi/fj4+KCi729WdqKi4sjICCATz75hF69egEw\nevRoJk2axIwZM5g9ezYff/yxWfoCuHfvHgAlSpRId71kyZJp93Jy+fJlANauXcvmzZtxcHDg4MGD\ntG3blqioqCwnKy5fvpw6derw8ssvm+EdmE6SAiEKkZgYOHfO+EmGMY9iOHTjEONbjTdbLP6e/hSz\nL8bW81uNSgrq1wdXV+2JhyQF5nc6IQHfiAiL9xPh60szYx9XZWP8+PF4enqmJQQAHTt2ZPz48Uyf\nPh37DMmHXq+nb9++JCUlAWQ6oCv1yZWqqri5ubF8+fJ091Pby9hucnIyKQac763T6QDw9fXF4ekS\noObNm5OQkMDChQt5//3305U/evQoa9euZfv27Tg6OubaviVJUiBEIbJ3r/arsUlB6p4CL9R4wWyx\nuDi60L56e7ae38r41oYnG/b22u6GMq/AMnxcXIjw9bVKP+Zw8+ZNfv75ZzZu3Jjuuru7OzExMSxc\nuJB+/fqlu2dnZ8eaNWtM7rPc0/Wwer0+3fX4+HhcXXNfmVO2bFkAvL290113c3PjwIH0m3o9evSI\nUaNG8dtvv9G8eXOTYzYXSQqEKETCw8HTE6oYuUPxlnNbaFqxKRVKmneNdJcaXZgUOonE5ESKOxY3\nuF67djB9Ouh0WpIgzMfF3t5s3+CtYf/+/ej1evz8/NJdT/0GHhsbi6enp1n7TB3Xj4mJwd3dHdCe\nKty/fz/TB31WqlevTokSJdKeGDzryZMn6X4eO3YsM2fOTHt/Fy9ezNd5BZIUCFGImLI/gV7Vs+38\nNt5o+obZ4+lSswvvb3+fsOgwo55C+PnBxIlw/Dg0aWL2sEQBotPpcHNzo3jx9EmlnZ0diqLw9ttv\nZ6qTcfgrxHqKAAAgAElEQVQgO9kNH3h6elKzZk3OnDlD/fr1ATh9+jRJSUl07Ngx15gdHR0JCAjg\n6tWr6d7HnTt3aJN6yAfw5ZdfMmTIkLSEIDo6mt27d0tSIITIu/h4OHQIhg0zrt6hG4e4nXDbLEsR\nM6pfrj4epTz4Pep3o5KCFi3AyUl78iFJQdHm7++PqqrExsamfWu/cuUKa9euBSAlJYXjx4/TsGHD\ntDp5HT4AbV+C4OBg+vTpA8DixYvp2bMntWrVAmDhwoVMmzaNffv2ZbkLYVBQEO+88w6ff/45dnZ2\nbNiwgTJlyjBq1CgAVq5cSWhoKA4ODkQ8neNx8uRJxowZk6e480qSAiEKif37tVn7xj4p2HR2E67F\nXGlb1YQtEHOhKArda3VnU9QmZnWdZXA9Z2do3lw7MfGtt8welihAypYty6+//sq4ceOoV68eqqpS\nsWJFZsyYgU6nY8KECQQGBqZLCszhX//6FxMmTGDcuHG4urpy8+ZNFi9enHZfVVWSkpIyzTtIFRAQ\nwEcffcSAAQOoVKkSt27dYs+ePbi6unL37l2GDx/O48eP0+17oCgK33zzjVnfh7EkKRCikNi7V5u1\n//Rpp8E2RW2iS80uONpbZtbzi7Vf5IdDPxB1J4paz9UyuF67drB0qbZtsxHbHIhCqGPHjlk+tp89\ne7bF+nRwcOCrr77K9v6YMWNy/VY/cOBABg4cmOl62bJl021cZEtkR0MhConwcG3Wvp0Rf6tjHsVw\n4PoBetTqYbG4Onl1oph9MTZFbTKqXrt2cP06PN3TRQhhBZIUCFEI6HTaIUIZJmjn6vdzv6Og0K2m\nCbsdGaiEUwk6eHYwOilInY8lSxOFsB5JCoQoBI4ehbg445OCTVGbaOnRknIlLHtOcY9aPdh9aTdx\nSXEG13nuOahXT5ICIaxJkgIhCoHwcG22vjF7nzzRPWHrua0WHTpI1aN2D5L1yWmbJBnKz0+SAiGs\nSZICIQqBsDBo2VKbtW+o8Ohw4p7E0aO25ZMC7zLe+Lj7sOms8fMKTp3SjoIWQlieJAVCFHCqqiUF\nRg8dnN1EpZKVaFqxqWUCy6BHrR5sPrcZvZr1Eq6stGun/Zq6fbMQwrIkKRCigDt/XjsIKfUD1FCb\nojbRo1YPo441zosXa7/IzUc3OXzjsMF1qlcHDw8ZQhDCWiQpEKKACwvT1vE/s3tqrs7dPceZO2es\nMnSQqm3VtrgWc2XD2Q0G11EULdmRpEAI65CkQIgCLiwMGjUCNzfD66w7vQ5nB2c6e3e2XGAZONo7\n0r1Wd9adWWdUvXbt4OBBSEy0UGBCiDRGJwWKovgpirJeUZRriqLoFUXpmUt5/6flnn3pFEUpb3rY\nQohU4eHGDx2sO7OOzt6dKeFUwjJBZaNXnV4cuXmEy/cvG1ynXTtIToYMJ84KISzAlCcFJYAjwD8B\n1cA6KlALqPj0VUlV1Vsm9C2EeMbNmxAVZdwkw9iEWPZe2UuvOr0sF1g2utXqhqOdI+vPrDe4TsOG\nULq0DCEIYQ1GJwWqqm5RVfVTVVXXAcbMULqtquqt1Jex/QohMkv9oDQmKdh4diOqqvJi7RctE1QO\nShcrTYBXgFFDCPb22nwJSQqEsDxrzSlQgCOKolxXFGWboihGTIkSQmQnPBy8vaFyZcPrrDuzjtZV\nWlOhZObjXq2hV51e7L68m/uP7xtcp107bVmiTmfBwISwsrt373Ljxg2uXLnC5cuX0175yRpJwQ1g\nNNAX6ANcAXYpiiKnpAuRR2Fhxs0nSExOZNv5bfkydJCqZ52epOhT2By12eA67drBw4dw4oQFAxMi\ng2nTpvH+++8zefJkhgwZQlyc4dt0p9qxYwcdOnTIdH3w4MG4u7vj4eFB9erV8fLywsvLi5YtW/Lk\nyRMzRG8aix+drKrqWeDsM5f+UhSlBhAEDMmpblBQEK6urumuDRgwgAEDBpg9TiEKmocP4cgRGDvW\n8DohF0JISE6gl0/+JQVVSlfBt5Iv686s47WGrxlUp0ULcHTUnow0bmzhAIUA5s+fz549e/j9998B\nmD59Oq+//jpr1641qP6vv/7K5s2befToEdHR0ZnuFytWjJUrV+Lo6Ijd06NN169fT5cuXXBycjI4\nzuXLl7N8+fJ01x48eGBw/UxUVTX5BeiBnibUmwHszeF+M0CNiIhQhRBZ27pVVUFVT582vM4b695Q\na8+tbbmgDPT5rs/VUlNLqY+THxtcp3VrVe3f34JBFWARERGq/JtpXtWqVVOXL1+e9vO1a9dURVHU\nqKgoo9pZvHix6uXlle6aXq9Xv/7663TXbt++rY4cOTLX9gz5s04tAzRTjfx8zq99CpqgDSsIIUwU\nFgblykHt2oaV1+l1bDi7IV+HDlL18ulF3JM4Qi+FGlynXTvtPauGrnkSwkRRUVFcuXKFevXqpV2r\nXLkyrq6uhIYa/v9sdhRF4c0330x37bPPPuPzzz/Pc9t5ZfTwgaIoJYCa/L3ywFtRlMbAXVVVryiK\nMg2orKrqkKflxwEXgZOAMzASCACst2uKEIVQ6nwCQ3cp3ntlL7fib9Hbp7dlAzNAw/IN8S7jzW+R\nv9G1ZleD6rRrBzNnQnS0tv2xKFri4+OZPn06Tk5OJCYmcufOHebOnYuTkxMJCQm4uLiYra/z58+j\nKAqlS5dOd71UqVJZDgWYolixYmm/DwsLo3LlylSsWNEsbeeFKXMKmgOhaI8mVODrp9eXAMPR9iGo\n+kx5p6dlKgMJwDGgk6qqe0yMWYgi78kT2LcPpk41vM7qU6upXKoyraq0slxgBlIUhb51+7L4yGK+\n6/Ed9nb2udZp21b7NTxckoK80CXoSDidYPF+XHxcsHfJ/c/VEHFxcQQEBPDJJ5/Qq5f2pGv06NFM\nmjSJGTNmMHv2bD7++GOz9AVw7949AEqUSL+5V8mSJdPumdNHH33Exo0bzd6uKYxOClRV3U0OqxZU\nVR2W4eevgK+MD00IkZ2ICHj82PCVB3pVz2+nf6OPTx/sFNvY3bxv3b589cdXhEWH0cGzQ67l3d2h\nbl3tCcnAgZaPr7BKOJ1AhG+ExfvxjfClVLNSZmlr/PjxeHp6piUEAB07dmT8+PFMnz4de/v0yYde\nr6dv374kJSUBpM5VS5N6CJiqqri5uWWaqJfaXsZ2k5OTSUlJMct7ShUREUFiYiJlypQxa7umsvjq\nAyGE+YWFQYkS0NTAU48PXDvA1YdX6Vuvr2UDM0ILjxZUKV2F1adWG5QUgByOZA4uPi74RvhapR9z\nuHnzJj///HOmb9Lu7u7ExMSwcOFC+vXrl+6enZ0da9asMbnPcuXKAVpy8az4+PhMK+Ly6n//+x8+\nPj5mbTMvJCkQogAKD4fnnwcHA/8Gr45cTTmXcvhVM2LrQwuzU+zo49OHVZGr+LbbtwY9wWjXDn74\nAe7ehbJlrRBkIWTvYm+2b/DWsH//fvR6PX4Ztu10ePo/f2xsLJ6enmbt08vLC4CYmBjc3d0B7anC\n/fv38fb2NmtfoaGhtDP28BILkqRAiAJGr9eSgnHjDCuvqiqrI1fT26e3QWP31vRKvVeYs38O+67u\n4/mqz+daPvXfzj/+gBetv0uzyAc6nQ43NzeKFy+e7rqdnR2KovD2229nqpNx+CA72Q0feHp6UrNm\nTc6cOUP9+vUBOH36NElJSXTs2DGP7yh9/8ePH6dLly5mazOvJCkQooA5dQru3TP8vIOjMUe5cO8C\nfXvYztBBqjZV21ChRAVWR642KCnw8oJKlbSkSJKCosHf3x9VVYmNjU371n7lypW0TYRSUlI4fvw4\nDRs2TKuT1+ED0HYcDA4Opk+fPgAsXryYnj17UqtWLQAWLlzItGnT2LdvHxUqZL9luE6nyzQMkeru\n3bvodDqjNiuyNEkKhChgwsO1YYNWBi4iWH1qNWWcyxDgGWDZwExgb2dPb5/erI5czVedv0qbAJYd\nRfl7vwJRNJQtW5Zff/2VcePGUa9ePVRVpWLFisyYMQOdTseECRMIDAxMlxSYw7/+9S8mTJjAuHHj\ncHV15ebNmyxevDjtvqqqJCUlZfuBv2XLFhYtWkR4eDi3b9/Gz88PHx8ffvjhh7QyJUqUwMPDgwYN\nGpg19ryQpECIAmbXLmjeXJtomBtVVVkVuYqedXriaO9o8dhM0bdeX76P+J5DNw7hWzn3CXDt28O7\n70JCAphxabqwYR07dszysf3s2bMt1qeDgwNffZX9wrkxY8YwZsyYbO937dqVrl1z3oPD2dmZK1eu\nmByjJdjG2iQhhEFUVUsKAgz80n/81nFOx57mH/X+YdG48qKDZwfcXdxZeXKlQeUDAiA5WTs1UQhh\nXpIUCFGAREZCTIzhScHKEysp41yGzjVsdwNRBzsH+tbtyy8nf8m0njwr9epp2zubYbdZIUQGkhQI\nUYCEhmqnBbZpk3tZVVVZeXIlfer2wcnediYyZaV/g/5cfnCZfdf25VpWUaBDB+2JiRDCvCQpEKIA\n2bULWrY0bD5BxI0Izt87T7/6/XIvnM/8qvlRsWRFVp4wfAjhwAF49MjCgQlRxEhSIEQBodcbN59g\n5YmVlHMpR4CX7a06yMjezp5/1PsHv5z6Bb2a9WzuZwUEQEqK7G4ohLlJUiBEAXHyJMTGGpYU6FU9\nv5z6hVfqvYKDXcFYZNSvfj+ux10nPDr3T/o6daBiRRlCEMLcJCkQooAIDQUnJ21749z8dfUvoh9E\nF4ihg1TPV32eKqWrGDSEoChaciSTDYUwL0kKhCggQkOhdWvIsNtrllaeWEnlUpVpV8129lTPjZ1i\nR7/6/VgVuYoUfe4n0XXooJ0W+fCh5WMToqgoGM8VhSji9HrYvRveeSf3sin6FFaeXEn/Bv1t7qyD\n3PRv0J+v//yanRd38kKNF3IsGxAAOp22u2GPHlYK0MZFRkbmdwjCwiz9ZyxJgRAFwLFj2nkHhswn\n2HlxJzHxMQxqNMjygZmZbyVfaj9Xm2XHl+WaFNSsCR4e2ryCop4UuLu74+LiwqBBBe/PXBjPxcUl\n7RwIc5OkQIgCYOdOcHY27LyDpceWUue5OvhWyn3LYFujKAqDGg5ixh8zWNB9ASWcsl97mTqvYMcO\nKwZoo6pVq0ZkZCSxsbH5HYqwAnd3d6pVq2aRtiUpEKIA2L5dOxXR2TnncvFP4vkt8jc+avdRrocL\n2aqBjQby6a5PWX9mPQMaDsixbGAgLFumrcqw0BenAqNatWoW+6AQRYdMNBTCxiUlafMJOhuwU/G6\nM+uIT47ntYavWT4wC/Eu402bqm1YenxprmUDA7XzIORpgRDmIUmBEDbujz8gMdGwpGDZ8WW0qdoG\n7zLelg/MggY2HMjWc1u5HX87x3IeHtpZCNu3WykwIQo5SQqEsHHbt2sHADVqlHO5W/G32HpuK4Ma\nFvzJZq/WfxVFUQw6ObFzZ+2/kQFnKQkhciFJgRA2bvt27TG5XS5/W1eeWImiKLxa/1XrBGZB7i7u\ndKvZjaXHch9C6NwZoqMhKsoKgQlRyElSIIQNu3NH26DHkKGD4GPBdK/VnedcnrN8YFYwqNEg9l3b\nx5nYMzmW8/fXTo6UIQQh8k6SAiFs2M6d2mPx3JKCE7dOcPD6QYY2HmqVuKyhZ52euDm7seTokhzL\nlSypbf0sSYEQeSdJgRA2LCQEfHygSpWcyy0+shh3F3d61C48u/g4OzgzoMEAgo8Go9PrcizbubO2\nDXRK7rsjCyFyIEmBEDZs+/bcnxIk65JZemwpAxsOxMneyTqBWcmwJsO4FneNkAshOZbr3Fk7A+HA\nASsFJkQhJUmBEDbq/Hm4eFGbZJiTLee2EBMfw9AmQ60SlzU1r9yceuXq8dORn3Is5+sLrq4yhCBE\nXklSIISN2roVHBy00wBzsvjoYppUbEKTik2sEpc1KYrCsCbDWHt6LfcS72VbzsFBS562bLFicEIU\nQpIUCGGjNm2C9u2hdOnsy8QmxLLhzIZCNcEwo0GNBpGiT2HFiRU5luveHfbt01ZsCCFMI0mBEDYo\nMVFbedC9e87l/nf8f4B2XkBhVbFkRbrV6pbrEELXrtoR01u3WikwIQohSQqEsEG7dsHjxzknBaqq\nsujQIl6q8xLuLoX7NKBhTYZx4PoBjsUcy7ZM5crQtCls3mzFwIQoZCQpEMIGbd4Mnp7acsTs7Lu2\nj+O3jjOq2SirxZVfXqr9EhVKVOCHiB9yLNe9uzavQJfzCkYhRDYkKRDCxqiqlhR07w45nX78Q8QP\nVHetTucaBmx3WMA52jsyrMkwlh5fSkJyQrblunfX5hTI0kQhTCNJgRA25uxZuHAh56GDh0kPWXFy\nBW80fQM7pWj8NR7RbAT3H99n1alV2ZZp1QrKlpUhBCFMVTT+NRGiANm0CZydISAg+zL/O/4/Hqc8\nZnjT4dYLLJ/VKFuDTl6d+OFQ9kMI9vbahMNNm6wYmBCFiCQFQtiYzZu1hMDFJfsyPxz6gR61euBR\n2sN6gdmAUb6jCI8O59TtU9mW6d4dDh2CGzesGJgQhYQkBULYkLg42LMn56GDiOsRHLpxiFG+hX+C\nYUa96vTC3cWdRYcWZVumSxdtLoZsZCSE8SQpEMKGbNsGyck5JwULIxbiUcqDrjW7Wi8wG1HMoRhD\nGw9l8ZHFJCYnZlnG3R1at4b1660cnBCFgCQFQtiQNWugUSPw9s76/v3H91l2fBmjfUfjYOdg3eBs\nxJjmY7j3+B4rT67MtszLL2ubGCVkv1BBCJEFSQqEsBFPnsDGjdoHWnaWHFnCE90TRjQbYb3AbEyN\nsjXoWrMr8w/Mz7bMyy9ru0Ju22bFwIQoBCQpEMJG7N4NDx5A795Z39erehYcXEDfun2pVKqSdYOz\nMW+2eJOD1w+y/9r+LO/Xrg316sHatVYOTIgCTpICIWzEmjVQvTo0bpz1/Z0Xd3L2zlnebPGmdQOz\nQd1qdsPTzZMFBxZkW+bll2HDBkhJsWJgQhRwkhQIYQP0eli3TntKkN0uhvMPzKdB+Qa0q9bOusHZ\nIHs7e8b4jmHFiRXEJsRmWaZ3b7h7F8LCrBycEAWYJAVC2IADB+D69eznE1x5cIX1Z9bzZos3UXLa\n+7gIeaPZGwD8dDjr0xN9faFKFe0JjBDCMJIUCGED1q7VltK1bZv1/e8OfkcJxxIMbFh4j0g2lruL\nO/0a9GP+gfmk6DOPESiKlmStXaudJyGEyJ0kBULYgDVr4KWXwCGLVYYJyQksjFjIiGYjKFWslPWD\ns2HvtHyHyw8us/5M1psSvPwyXLmi7XAohMidJAVC5LNTp+DMmexXHfx89GfuP77P2y3ftm5gBYBv\nZV/8qvkx+6/ZWd5v3x7KlIHVq60cmBAFlCQFQuSzFSvA1RVeeCHzPb2qZ/a+2fSq0wuvMl7WD64A\nGN96PGHRYURcj8h0z9ER+vSBlStlCEEIQ0hSIEQ+UlXtA6t3byhWLPP97ee3czr2NEGtg6wfXAHR\nq04vPN08+Xbft1ne799fO4r64EErByZEASRJgRD56MgROHtW++DKyux9s2lWqZksQ8yBvZ09b7d8\nmxUnVnAjLvPRiB06QPnyWvIlhMiZJAVC5KMVK+C556Bjx8z3Im9HsuXcFsa3Gi/LEHPxRtM3KOZQ\nLMvNjBwc4B//0JICvT4fghOiADE6KVAUxU9RlPWKolxTFEWvKEpPA+p0UBQlQlGUx4qinFUUZYhp\n4QpReKQOHbzyijb2ndHXf35NpZKVeLX+q9YProBxdXZleJPhLDi4gPgn8Znu9+8PV6/CH3/kQ3BC\nFCCmPCkoARwB/gnkOnVHURRPYCOwA2gMfAssUhSlswl9C1Fo7NsHly9nPXRwPe46Px/7mfGtx1PM\nIYvJBiKToOeDePD4Af89/N9M99q0AQ8P7cmMECJ7RicFqqpuUVX1U1VV1wGGPNMcC1xQVfVDVVXP\nqKo6H1gFyMwpUaStWAEVK4KfX+Z7c/bNoZh9MUb7jrZ+YAWUp5sn/Rr045u/vsm0mZGdHfTrB7/+\nKmchCJETa8wpaA2EZLi2FXjeCn0LYZN0Ou0D6tVXwd4+/b2HSQ/5/uD3jGk+Bldn1/wJsID6oM0H\nXLp/iVWnVmW6178/3LoFu3ZZPy4hCgprJAUVgZgM12KA0oqiyHNRUSTt2KGddfDaa5nv/RDxAwnJ\nCYxrNc76gRVwTSo24YUaLzBj7wzUDBsTNG8OtWpBcHA+BSdEAZDFpqq2IygoCFfX9N+UBgwYwIAB\nA/IpIiHMY8kS8PGBli3TX3+ie8Ksv2YxsNFAPEp75E9wBdwHbT6g88+dCbkQQucaf09dUhQYMgSm\nToX586GU7BgtCoHly5ezfPnydNcePHhgcntKxmzaqMqKogdeVlU1643HtTK7gQhVVd995tpQYJaq\nqmWyqdMMiIiIiKBZs2YmxyeELXrwQJtLMHky/Otf6e/99/B/eWP9G5wYe4L65evnS3wFnaqqNP+h\nOW7ObuwYvCPdveho8PSEH3+EYcPyJz4hLO3QoUP4+voC+KqqatTJH9YYPvgT6JTh2gtPrwtR5Pzy\nCzx5AoMGpb+eok9hWvg0+tTtY5sJQUwMzJkDrVpBtWrw7rsQEWFz+wcrisJEv4nsvLiTP66kX4NY\nrRp06gSLF+dPbELYOlP2KSihKEpjRVGaPL3k/fTnqk/vT1MUZckzVb5/Wub/FEWpoyjKP4FXgG/y\nHL0QBdDixdo5Bx4ZRgd+OfkL5+6eY6LfxHyJK0ezZ2sBv/++9pijZ09YtkwbqA8MhPjMewPkp5d9\nXqZeuXp8GfZlpntDh8KePXD+vPXjEsLWmfKkoDlwGIhA26fga+AQ8O+n9ysCVVMLq6p6CegBBKLt\nbxAEvKGqasYVCUIUemfPahvoDB2a/rpe1TM1bCrda3WnWSUbGzJbuBCCguDtt+HGDVi3DubNg2vX\ntDOf9+/XDm94/Di/I01jp9gx0W8im6M2c/jG4XT3evfW5hPIhEMhMjNln4Ldqqraqapqn+E1/On9\nYaqqdsxQZ4+qqr6qqhZXVbWWqqo/m+sNCFGQBAdrJyL26pX++rrT6zh5+yST/CblT2DZWboUxo6F\nd96Bb77R9mRO5eAAL78MGzZAWJi25i85Of9izeDV+q9Ss2zNTE8LXFy0paDBwbLtsRAZydkHQliJ\nTqetOujXD5yd/76uqipTwqYQ4BnA81VtZ/uOR7t2sXnBAr78+mv+MXIkjQ8epPqff+IaFkbxPXuo\n9Mcf1N2/n0A3N95fv56lSUlcmjAhv8NO42DnwEdtP2J15GpO3jqZ7t7QoXDpkuxZIERGeVp9YCmy\n+kAURhs2aEPxBw+CNjFYs/HsRl5a/hI7Bu+go1cWJyNZUVxKCitv3WLt7duE3LpFkqMjbg4ONClZ\nkrouLjz39GdHReFBSgr3UlK49PgxRx494uLT4YNGisLL1arxWoUK1HFxydf380T3hFpza9GmahuW\n9/172ZaqQr160KiRnJ4oCp+8rD6w6X0KhChMvvtOm5f3bEKgqiqfhn6Kf3V/AjwD8i22C4mJzL12\njf/euMEjnY52jx4xLTiYlz78kBq+vgad0ng/KYmQUaNY27gx39rZ8fnly3QtW5ZxHh68ULYsdvlw\n0qOTvROT/CYxeuNoJvpNpEH5BoC2Z8GYMdq8yZs3tbmTQggZPhDCKi5ehC1btOH5Z607s47DNw/z\n7w7/zpfjkW8kJTHyzBlq7dtH8M2b/NPDg0s1a7J74ECCnnuOms2bGxyXW7FivDJiBEvfe4+YqCgW\n+/gQ8+QJ3Y4fp2VEBLvu3bPwu8na0CZD8XTzZPKuyemuDx6snU7538znJwlRZElSIIQV/Oc/ULq0\nNp8glV7V89muz+jk1Ql/T3+rxvNYp+Pfly5Rc98+frt9m5k1anDl+eeZ5u1N1U8/1SYRTplifMN+\nfjBwIMU++oghzs5E+Pqys3Fj7BWFgKNHeen4cc4lJJj/DeXA0d6RT9p/wurI1Ry5eSTtepky2tzI\nhQu1+R5CCEkKhLC4J0+0HfSGDIESJf6+vvrUao7FHOPfHf6dfWUL2PfwIc0iIvjy8mXe9PDgfKtW\nBFWtiou9PRw5Aj/9pCUEz640MMaMGZCYCFOnoigKAWXK8FezZqyoV48T8fE0OniQ2VeuoLPifKbX\nG79OzbI1Mz0tGDtW2+VwyxarhSKETZOkQAgL++03uH1bG8NOpdPrmLx7Ml1qdKFttbZWiSNZr+ej\n8+dpc+gQLnZ2HPL1ZUaNGrg5Ov5d6OuvoXp1GDnS9I4qV4a33oLvv4eHDwFtl8F+5ctzvHlzRlSq\nRND58/gfPszFxMQ8vivDONg58Jn/Z6w7s46D1w+mXW/RQpvj8d13VglDCJsnSYEQFvbdd+DvD3Xr\n/n1t6bGlnLp9is8DPrdKDNeSkgg4coSvr15lipcXfzVrRoOSJTMUugYrVsC4cdrwQV689Zb2tCDD\ngH1JBwfm1KrF7iZNuP7kCc0iItgYG5u3vgw0oMEA6rrXZcKO9Msmx4yBzZu1eR9CFHWSFAhhQYcO\naVvqvvXW39cepzzm012f0rduX1p6tMy+spnsvHePZgcPcunxY3Y3acKE6tVxsMvir/78+VC8OLzx\nRt479fDQJlB8+22WA/bt3dyI8PWlvasrL504wccXLlh8OMHezp6pnaYSciGEkAt/b6j62mva/IK5\ncy3avRAFgiQFQljQrFnaqXy9e/997bsD33Ht4TW+7Jh5X35zW3T9Oi8cPUqjkiU53Lw5bTIcRZ4m\nPl6bcTdihDYj0hyCgrQdgtauzfJ2GUdH1jZowP95e/N/0dH0OXGCeAvP+OtVpxetq7Rmwo4JpO7R\n4uICo0fDokVpox1CFFmSFAhhIc8+jbe31649ePyAL8O+ZHjT4dRxr2OxvvWqysQLFxh59iyjKlfm\n94YNKefklH2F4GC4f1/bzthcfH2hfXstM8qGoih8WK0aGxo2ZMe9e3Q4coSbSUnmiyGL/qZ3ms7B\n6xZLl/UAACAASURBVAdZdWpV2vW33tKObvjxR4t1LUSBIEmBEBYyb572LXT48L+vzfxjJvHJ8Xzm\n/5nF+k3W6xkcGcnU6Ghm1qjB/Fq1sh4uSKXXax/cffpojzXMKSgI9u7VDk3KQffnniOsaVOuJSXR\n+tAhiy5b9Pf0p1vNbkzcOZFknXZWQ+XK2vLEb7+FlBSLdS2EzZOkQAgLePRIm3z/7NP463HX+eav\nbxjXahwepT1ybsBESXo9r546xcrbt1lZrx7vVa2a++ZDoaEQFaU90jC3l14CLy+Dpvc3LVWKv5o1\nw9nOjvZHjnDKgscxT+s0jXN3z/HDoR/SrgUFweXL2sGPQhRVkhQIYQFLlmjj088+jZ+0cxIuji58\n1O4ji/SZoNPx8okT/H7nDmsbNODV8uUND7ZWLWhrgaWR9vbaBg2rVmnzFnJRzdmZ3U2b4u7oiP+R\nIxyOizN/TEDjio0Z0mQIn+36jAePHwDQtCl06KAdBilEUSVJgRBmlpKifbC88oq25B/g8I3DLD6y\nmH93+Dduzm5m7zNRp6Pn8ePsuX+fTY0a0cPQjYcePYLVq7UPbktts/z661o/2Uw4zKiCkxO7mjTB\n09mZjkePWiwxmBIwhYTkBKaGTU279u678Ndf2knQQhRFkhQIYWYrVsCFC/DR0wcCqqry3rb38HH3\nYZTvKLP390Sv55WTJ/nj4UM2N2pEpzJlDK+8ejUkJMCgQWaPK423t7b98ZIlBlcp6+hISOPG1Cxe\nnBeOHeOkBYYSPEp78GGbD5m9bzYX72mbFPToAQ0awJeWXxgihE2SpEAIM9LrYdo06N5dexwNsOHs\nBkIvhTLzhZk42Jn3YNIUvZ7XTp0i5N491jZogL+bkU8hgoMhIODvRxqWMmQIhIRoSzIM5OrgwNZG\njajs5ETg0aNEWWDy4ftt3sfdxZ2PdmgZnJ0dfPwxbN2qHXEtRFEjSYEQZrR2LZw6BZMmaT8npSTx\n/rb3CfQOpFvNbmbtS1VVRp49y9rYWH6pX58XypY1roHoaG2S4eDBZo0rS6+8AsWKwdKlRlUr6+jI\n9saNcbW3J/DoUa6ZebliCacSTO04lV9O/kJ4dDgAr74KNWvC1Km5VBaiEJKkQAgzUVXtHKGAAHj+\nee3a7L9mc+HeBWZ3mW32o5EnXbzI4ps3WVK3Lr3c3Y1vYOlSbQfDvn3NGleWXF21HZyCg7X/UEYo\n7+RESOPGqEC3Y8e4n5xs1tBeb/w6LSq34O3f30an12Fvrw39rFkDJ0+atSshbJ4kBUKYyZYtcPgw\nTJyo/Xzt4TW+2PMFb7d8m/rl65u1r3lXrzI1OpqvvL0ZWKGC8Q2oqvYB3acPlCpl1tiyNXiw9hgl\nIsLoqlWcndnSqBFXk5LoffIkSXq92cKyU+yY130eR24e4T8R/wG0uZFVq2pDQUIUJZIUCGEGqgqf\nfw6tW0PHjtq1D7Z/QAmnEkzuMNmsfa25fZt3zp0jqEoV3qta1bRGjhyBM2csO8Ewo8BAqFBBm4lp\ngnolSrC+QQP+eviQIZGR6M14VkJLj5YMbzKciTsnEpsQi5MTfPghLF+u/WcSoqiQpEAIM9i0SVvK\n9vnn2sq+PZf3sPzEcqZ3mo6rczbnDZjgSFwcgyIj6VuuHDNr1DB9SGLVKihb9u8MxhocHLQnE6tW\nGT2EkKqdmxvL6tZl5e3bfHH5slnDmxY4Db2qZ9JObULIiBHaToeTJ5u1GyFsmiQFQuSRXg+ffKId\njxwYCMm6ZN7c/CatPFoxpMkQs/UT8+QJPU+cwMfFhSU+PtiZmhCoKvz6K7z8Mjj+P3vnHVdV/f/x\n5+EOLlyWLAeKgAMBcZtaacPSlpqaq7JpWdke3297fTOzvYctf1aO0kq0LC0rtcy9kKECioBshLu4\n8/z++ACKA0UOQzzPx+M87uVyzufzhnvvOa/zeS+dYvadFtddJ8oGNiC0f1xYGP+LiuL5/ftZXFio\nmGnhxnBevORF5myZw6bcTRgM8OyzYmFj507FplFRadGookBFpYEsWSJW4196SawSvPXvW6QUpfDh\n1R/iJSnzFbN7PIxLTsYpyyzt2RPf6g5LZ8KuXaKs8XXXKWJbvRg2DEJDxWpBA3iqc2cmh4dzU1qa\nosWN7hl4D73b9Wb68um4PC5uuQW6dBHiQEXlXEAVBSoqDcDtFheMK66ACy+E/Yf38/yfz/PAoAfo\n176fInPIsszde/awxWTih4QEOhoMDRtw8WIICoLhwxWxr15UuxC+++6MXQgguh1+HhtLvK8vY5KT\nFeusqPXS8sk1n7A9fzvvb3wfnU64D5YuPWVPJxWVVoEqClRUGsA330BaGvzvf+LiPePnGYT4hvDi\nJS8qNsdbOTl8mZ/Pp7GxDA5UID5h8WIYMwbqaqXcmFx3HWRliVSNBuCr0bA0MRGXLDN2924q3W5F\nzDsv4jzuGXgPz/zxDAfLDzJlCsTHH6k9oaLSmlFFgYrKGWKziViCceNgwABYkrqEn/f+zHtXvoef\n3k+ROVaUlPBYRgb/7dSJqe3aNXzAlBRITW0e10E1F18sghwb6EIAiPD25seePdlmMjF9zx5khTIS\nZl46E3+9P/f/cj8ajSh7vGoVrFypyPAqKi0WVRSoqJwh774LeXnwyitQZivjvhX3MTp2NNf2uFaR\n8dOtVianpHBVSAgzY2IUGZPFi0VdgssvV2a8M0GnE4WMGuhCqOa8gAC+6NGDeQUFvJmTo4CBEGgI\n5J0r3uHHtB9ZkrKEMWOEe+ixx4TLSEWltaKKAhWVM6CoSJTBvece0XX44ZUPY3Va+fCqDxUZ3+J2\nMz45mQ7e3nwTF4dGqWqIixfD6NGi5HBzct11sG+fYmH917dty386deK/GRmsOXxYkTGvi7+Oa3tc\ny4yfZ1BWWcrrrwtz581TZHgVlRaJKgpUVM6AF6tCBp55Bn7d9ytzt8/lzRFvEhEQ0eCxZVlmeno6\n+ysrWZKQQIBWoSZKGRki82DcOGXGawiXXipKH59mO+XTYWZ0NEODgpiUksIhBQIPJUniw6s+xO62\n89CvDzFoEEyaJGILGqFpo4pKi0AVBSoq9SQ9HT7+WJQz9vY3cefyO7ks5jJu63ubIuN/kpfHN4WF\nfBobS7zRqMiYACxbJlYIRoxQbswzRa+HK68UYf0KofXyYkFcHBIwJSUFlwKlkNv7t+fNEW8yb8c8\nft77M7NmQXExvPFGw+1VUWmJqKJARaWePPwwRETA/ffDf3/7L8XWYuZcM0eRhkebKyp4YN8+ZnTo\nwJQz6WlQF0uXijREP2WCIBvMmDEiAyE7W7Eh23l7szA+nnXl5TydlaXImLf0uYURXUYwffl02rQ/\nzH33wezZcPCgIsOrqLQoVFGgolIPli+Hn3+GN9+Ev3J+5aPNH/HqZa8S3Sa6wWOXOJ1ct3s3ffz8\neKNrVwWsPYrSUli7VsQTtBSuvFIEHSYlKTrssKAgZsXEMPvgQZKKixs8niRJfDrqUyrsFTzwywM8\n84yI1fzPfxQwVkWlhaGKAhWV06SyEh58UJQyvviKMm5Luo3LYy7nnoH3NHhsjyxzU2oqJrebbxMS\n8PZS+Ku5YoUIm7/mGmXHbQiBgSI9UUEXQjWPdurEtaGh3JSaSqbN1uDxIgMjee/K95i3Yx6/533P\nK6+I8sd//aWAsSoqLQhVFKionCZvvSXK9r/7Ltz3y71YHBa+GPOFIm6DWdnZrCgt5Zu4ODo3tGLh\niUhKEsUUIhoeCKkoY8bAn3+CQhkD1UiSxJexsYTqdFynUGGjqb2mcm2Pa5m+fDojxxUweDDcdx+4\nXAoYrKLSQlBFgYrKaZCTI3ob3H8/7PJ8y/xd8/ngqg/oGNCxwWP/WVbGs1lZPN25M1eEhChg7TE4\nHGKloCW5DqoZNUpcVX/5RfGhg3Q6FickkGq18uC+fQ0eT5IkPrnmEyQkpv90B+++K5OcLIJOVVRa\nC6ooUFE5DR54AAIC4JYHD3DnsjuZED+B6xOvb/C4xQ4HN6SmMiwoiOeiohpu6In46y8wmVqmKIiM\nhL59G8WFANDH3593unblk0OHFOmoGG4M57PRn7FszzI28RHTpokUxUOHFDBWRaUFoIoCFZVTsGwZ\nfP89vPGWi7tX3UCgIZA5oxqebSDLMrekpeGQZWULFB1LUpK4+Pbq1TjjN5QxY0T0psPRKMPf0b49\nE8LCmJaeTpYC8QWjY0dzz4B7eGTlI9z8WDJ6PTz0kAKGqqi0AFRRoKJSB2Yz3Huv6IKY3vYl1ues\nZ/64+QQZgho89ts5OfxUWsr/9ehBh8aqMCjLQhSMHi36OrdExoyBiopGi9qTJIk53bvTRqdjSkoK\nTgXqF7w+4nW6tOnCXb9N4ZXXbSxa1CgeEBWVJkcVBSoqdfD886Kk8S3PreGltf/j2WHPckHkBQ0e\nd3NFBf/NzOSRjh25qjHiCKrZuVPUAWiJroNqeveGzp0bzYUAIr5gYXw8W8xmnlGgfoGPzocF4xew\nt2Qvm0IeYfhwUfLaalXAWBWVZkQVBSoqJ2HbNnj7bXj4mUIe+nsyF0ZeyFPDnmrwuBUuF5NTUujt\n58fLSjU6OhlJSSKp/qKLGneehiBJQrQkJSnSIOlkDAoI4OXoaGYfPMjK0tIGj5fYNpG3Rr7Fx5s/\n4prHF5GXBy+8oIChKirNiCoKVFROgNMJt94K8Qlu1ne4AbfsZsH4BWi9GtaHQJZl7tqzh0Knk4Xx\n8eiVrkdwLElJokiQXt+48zSUMWNEicBt2xp1mkc6dWJkmzZMTU0lX4H+CHcNuIvJPSfzzKZpzHh2\nD6+/Dps2KWCoikozoYoCFZUTMHs2JCfD+Y+/xB/7f2f+uPl08O/Q4HHn5uezoLCQOd2708XHRwFL\n6yA3FzZvbtmug2qGDRPFjBSubngsXpLEvLg4vCSJqWlpeBq4MiFJEnOumUOEfwQrg66jVz8bt93W\naDGTKiqNjioKVFSOITlZdEGc8PhK5ux5gecvfp7hMcMbPG6qxcK9e/cyrX17Jivd1+BELF8OGo1Y\nKWjp6HRw1VWNGldQTbhez9dxcfxeVsZsBfou+Hv7s3jiYjLK9tHx7rtITZOZOVMBQ1VUmgFVFKio\nHIXLBbfdBpG9M/nVfzIju47kqaENjyOwud1MSkmhs8HAO0r3NTgZSUkwdCgEBzfNfA1lzBjYvl2U\njWxkhrdpw5ORkTyTlcU/5eUNHq9neE8+HfUpyw/OY/iTb/Pyy7BjhwKGqqg0MaooUFE5itmzYfNO\nC9LksQT7BDN/3Hw0XpoGj/tIRgZ7bTYWxcfjq2n4eKfEbIbffxcX2rOFRmqQdDKej4picEAAU1JS\nKHM6GzzeDb1u4LHzH+M3zaN0vGglN90ECoQtqKg0KaooUFGpYutWeO55mbj/TOOQPYMfJ/9IG582\nDR53SVERH+Xl8XbXriQ2VdviVavEFWnUqKaZTwkCAuCSS5rEhQCg9fJifnw8JrebaenpyApkPswa\nPouRXUZScukkUvL38txzChiqotKEqKJARQWw2eDGGyF83CukaBby5Zgv6Rnes8Hj7rfZmJaezoSw\nMO5s314BS0+TpUshIQG6dGm6OZVgzBhRxEjhBkknI9Jg4IvYWL4vLuajvLwGj6fx0rBg/AI6BLSl\nzYzRzH6nnHXrFDBURaWJUEWBigrw5JOwV/s9h+Kf5LmLnmNCwoQGj+n0eJiSmkqQVsuc7t0V6aZ4\nWrjdIsjwbMg6OJbRo0Vgx4oVTTbltWFh3BsRwcP79rHDbG7weIGGQJKmJOHQH6LN7Tcw9WY3JpMC\nhqqoNAGqKFA551m5Et5etBWv8VOZlDCJ5y5SZs332f372WwysTA+niCdTpExT4v166Gk5OwUBR07\nQr9+TeZCqOa1mBjijEYm7d6NRYE2y91DurPoukWUh68gN/YpZsxQwEgVlSZAFQUq5zSFhXDD3Qfx\nvmU0vTsk8OWYLxW5o19ZWsor2dm8HB3NoIAABSytB0lJEB4O553XtPMqxZgxYqWgCZP9DRoNC+Pj\nybHbuXfvXkXGHNl1JK9e9irOQbP5auc85s5VZFgVlUbljESBJEkzJEnKkiTJJknSv5IkDaxj34sk\nSfIcs7klSQo/c7NVVBqOxwNTbiul7JqRhIfqWDp5KT66hhcUyrfbmZqaysg2bXikUycFLK0nSUki\nwLCxqyU2FtUNkv78s0mnjfX15cPu3Zmbn8/X+fmKjPnwkIe5ve/tSNfezl2v/0JKiiLDqqg0GvU+\na0iSNAl4A3gO6AvsAH6VJCm0jsNkoBvQrmprL8tyw5ubq6g0gNfetrK67SiM4UX8dvOvtPdveCCg\nR5aZmpZWq3Jek5KeLraz0XVQTa9ejd4g6WTc1K4dU9u25e69e9mrQHcjSZL4+JqPuaLLFTjGjmfU\nXRvVpkkqLZozuZV4CPhEluV5siynAXcBVuC2UxxXJMtyYfV2BvOqqCjGP/+6eGLrJLSdtrPqlp/o\nHtJdkXFnZ2fze1kZX8fFEd4c/QaWLQODAS67rOnnVgpJEqsFjdwg6WR82K0b7fV6JqWkYFegzbLW\nS8viyYvoFd6bzCFXc/Mj6QpYqaLSONRLFEiSpAP6A79XvyaL5N7fgCF1HQpslyQpT5KklZIknX8m\nxqqoKEFxscyI96ZDl1/4ftISzotQxvf+T3k5z2Rl8WRkJMPbNLy+wRmRlASXXw6+vs0zv1KMGQM5\nOaJ4RBPjp9WyKD6e3RYL/83IUGRMX50vq6ctp31AGIsNI/lgXsPTH1VUGoP6rhSEAhqg4JjXCxBu\ngRNxCJgOjAfGAQeBPyVJ6lPPuVVUGozHA4OfehpL9y948+IvGRV3hSLjljqdTElJYXBAAM9HRSky\nZr0pLoa//z67XQfVDB0KQUFNVt3wWPr6+/N6ly68k5tLUnGxImMG+wTz74xf8TW6uX/DlWxJbnh5\nZRUVpWlYH9jTQJblPcCeo176V5KkLgg3xM11HfvQQw8RGBhY67UpU6YwZcoUxe1UOTcY8/J7ZHR4\nmdsjX+PBS25UZExZlpmWno7J7WZ+fDza5grw+/lnoXquuaZ55leSoxskvfBCs5hwb0QEv5eVcWta\nGtsHDKCTwdDgMSODOvHH7b8w5LMLueiTMeTM+oUgv4aPq3LusmDBAhYsWFDrtfIG9POQ6lPas8p9\nYAXGy7KcdNTrc4FAWZbHnuY4rwIXyLJ8wUl+3w/YsmXLFvr163fa9qmo1MWDX83hnczpDPY8wvoX\nXlds3A9yc7l3715+SEjg2rAwxcatN+PHiyX3DRuazwYl+fZbmDQJsrKgmVZfSp1O+mzeTGeDgT96\n91ZM8M1d/Te3rr6cDs5h7HvpR3x0qjBQUY6tW7fSv39/gP6yLNfLB1evT7gsy05gC1DTR1YSSd3D\ngX/qMVQfhFtBRaVJmPmLEASd8+9l7TOvKTbuDrOZR/bt496IiOYVBFYr/PILjD0tXX52cMUVTdog\n6UQE63TMj4tjfXk5LyjYvfGWSy/giahl5OnW0Hf2tVS6KhUbW0WlIZyJ7H0TuEOSpJskSeoBfAz4\nAnMBJEmaJUnS/1XvLEnSA5IkjZYkqYskSQmSJL0NXAK833DzVVROzfvrP+XpDdMJSLuXrTPfRatV\nJk2w3OXiut27iTMaeS0mRpExz5iVK4UwGDeuee1QkoAAuPTSZklNPJoLg4J4ITqamQcOsLK0VLFx\nX542nMnyMtIr13Dh+6owUGkZ1FsUyLL8LfAo8CKwDegFjJRluahql3bA0RVb9Ii6BjuBP4FEYLgs\ny3+esdUqKqfJnC2fct/KO9FuvZd1T75LcLAygkCWZW5PS6PI4WBxQgKGpmiHXBc//ADx8dBdmdTK\nFkN1g6SysmY144nISK4IDub6lBQOVip38f76xeEMzlzGluI1jPh8rCoMVJqdM3KQybL8oSzLUbIs\n+8iyPESW5c1H/e5WWZYvPern12RZ7ibLslGW5TBZlofLsrxGCeNVVOris62fMX35nbBxBgtvepfE\nROUKCb2Tk8OS4mLm9uhBF5+GV0FsEE6nWGJvTasE1YweLRo8/fxzs5rhJUl8FReHr0bDxJQUHArU\nLwDQaODXj4cTuW4Z63L+4pqvVWGg0rycpXVQVVTqZs6WOdyx7A7YOIOXzn+P8eOVEwT/lJfzWGYm\nj3bq1LxxBNVUtxpuTfEE1UREwIABzRpXUE2ITsd3CQlsMZn4j0L1C0B4SX7/bDjGpGX8kfkX1y4Y\ni9Wplj1UaR5UUaDSqpBlmZlrZjJ9+XQ0m+/jxpD3ePJJ5QRBkcPBxN27GRwQwMvR0YqN2yB++EGU\nBe7bt7ktaRxGjxYNkuz25raEQQEBvFlVv+C7QuUKs3btCj++ORzmL2d1xlpGfDWCUpty8QsqKqeL\nKgpUWg0e2cP9K+7n6T+exnfDi5xX+g6ffSqhVPsBtyxzQ2oqTllmUXw8upbQcMjjEaJg7FgU+0Nb\nGmPGgMnU5A2STsaMiAgmhYVxW3o66Qo2Mhg+HN5+4FKcn61me04aQ78cSk5FjmLjq6icDi3grKai\n0nDsLjtTlkzhw80fEvbvHCL2PUPSUglvb+Xm+N/+/fxeVsb8+Hg6KDlwQ9i4EQ4dap2ug2oSE0Wd\ngmbOQqhGkiQ+jY0lQq/nut27sbrdio19771w//jzsL7/N0WHLZz/+fmkFqUqNr6KyqlQRYHKWU+F\nvYKr5l/F0rSldNm8BGnrHfzyC4TW1beznqwoKeHFAwd4MTq6+foanIglSyAsDC44YR2w1kF1g6Sl\nS8XKSAvAX6tlSc+eZNpsTN+zh/oUgasLSYI334Rxw2Ixvf0Pek8QF355IesPrldkfBWVU6GKApWz\nmnxzPhfPvZgteVvot3sleauv5aefQMmyAXusVqakpHB1SAhPREYqN3BDkWVR9W/cOBHG3pq57jrI\ny4N/6lMjrXFJMBr5LDaWrwsKeDtHuWV+jQa+/hoGxHag9I01RBsTGD5vOD/t+UmxOVRUToYqClTO\nWrYd2sZ5n55HvjmfC9LXsnnJMJYsEcHqSlHhcjEmOZl2ej1fx8Xh1ZL89v/+C9nZMHlyc1vS+Jx/\nvshEWLSouS2pxZS2bXmsUycezcjgNwULGxkMYmEkIiSIQ6/+ytD2Ixm9cDRvrn9TsVUJFZUToYoC\nlbOSb3d/ywVfXEBbv7ZclrWJX+Yl8s03MHKkcnN4ZJkbU1M5ZLezNDGRQG2j9w+rH4sWQbt2oqNg\na8fLS/RB+O47UbegBTErJobL27RhYkoKGTabYuMGB4tClQatDxmvLObuXo/yyMpHuHXprWotA5VG\nQxUFKmcVHtnDM6ufYdLiSYztMZZhGWv46oMIPvkEJkxQdq7n9u9neUkJ8+PjifX1VXbwhuJ2C9fB\nxImt33VQzaRJUFAg6jK0IDSSxIL4eEJ0Oq5NTsbscik2dvv28NtvYDVrWPPsbD68/CsWJi/kkv+7\nhEMmtX2MivKookDlrMFkNzFu0Thmrp3JrOGvELXta96c7cObb8K0acrOtbiwkJcOHGBWTAxXhYQo\nO7gSrFsnsg4mTWpuS5qOgQMhOhoWLmxuS46jjU7H0p492V9Zyc1paXgUXOKPjhbCID8f5sy4keXj\n15Jdns3ATweyOW/zqQdQUakHqihQOSvYV7qP8784n9VZq0makoT99//y8kyJ116Dhx5Sdq6dZjM3\np6UxOTyc/3TqdOoDmoNFi6BTJxg8uLktaTokSYigJUtEaecWRrzRyNdxcXxfXMxMBTsqgmhr8dtv\ncPAgPD51IL9N3EREQARDvxzK1zu/VnQulXMbVRSotHgWJi+k3yf9qHRV8s9t69k8/xqefx5mzYJH\nH1V2riKHgzHJycT6+vJ5bCxSSwosrMblgsWLxQWyJRRQakomT4bSUnGFbIGMCQ3lxagont2/nx+K\nik59QD3o1Uv82ZmZcPO4Dvww+i8mJkxk6g9TmZY0TS2NrKII59gZReVswuq0cueyO5myZApXd7+a\nzXds4f9eT+CFF4QgePxxhedzuxm1axc2t5sfevbEt6X66v/4A4qKzi3XQTW9ekFsbIt0IVTzVOfO\nTAgL44bUVDZWVCg6dp8+8PvvQhhcebmBV4fM5YvRXzB/13wGfjqQ3YW7FZ1P5dxDFQUqLZKUohQG\nfTaIr3d+zaejPuXra+fz5CMBvP46vPOO8oLAXZVpsMtiYXliIp0NBmUnUJJvvhHF8vv3b25Lmh5J\ngilT4PvvwWJpbmtOiJck8X89etDHz49rdu0iU8GMBBAtLv76CwoL4eKLJUa2vZVNd2xCQmLgpwP5\nYtsXatqiyhmjigKVFoUsy8zdPpeBnw7EI3vYeMdGbk6cxi23SHz0EcyZA/ffr/y8j2ZksLS4mIXx\n8QwICFB+AqUwmURa3s03t95eB6di6lQwm0XPhxaKj0ZDUs+eBGq1XLVzJyUKx0AkJMCaNUIXXXAB\n6MsT2HjHRm5IvIHbk27nxh9uxGQ3KTqnyrmBKgpUWgwF5gLGfzueW5feyuSEyWyctpEo356MGiXi\n6hYsgDvuUH7etw8e5O2cHN7r1o1RStZGbgwWLwabDW66qbktaT5iYuCii+DLL5vbkjoJ1etZkZhI\nicvFtcnJVCpcX6FbN5GEYjAIYZCyw5dPR3/KN+O+ISk9icSPElmdtVrROVVaP6ooUGl2ZFlmUfIi\nEj5MYF32OhZPWMznYz7HctjIpZeKyrYrVjSOC/37oiIezsjgP506cU9EhPITKM3cuaKdXksqt9wc\n3HorrF4NCkf5K01XX1+W9ezJZpNJ8VRFEB+DdeugSxe4+GL49Ve4PvF6dty1g6igKIbPG869P9+L\n2WFWdF6V1osqClSalSJLERMXT2TykslcEn0Ju+/Zzfj48aSkwKBBoorvn3+K66DSrC8v54bUVCaG\nhTFLyWYJjUVGhlgzvvXW5rak+Rk/HoxGmDevuS05JYMDA5kfF8d3RUU8npmp+PghISIr4ZJLYTuI\nTgAAIABJREFU4Oqr4eOPIaZNDKtvXs27V7zLF9u+oPfHvVlzYI3ic6u0PlRRoNIsyLLMd7u/I+HD\nBP7I+oNF1y3iuwnfEWYMY9UqGDIE/P1FZ+B+/ZSfP8ViYdSuXQzw92dujx4tq6fByfi//4OAALj2\n2ua2pPnx8xMlLOfOFY2hWjhjw8J4q2tXXjt4kLcOHlR8fKNRhFjMmAF33w0PPwyyx4v7Bt3Hzrt3\n0sG/AxfNvYgHVjygrhqo1EkLK+auci6QUZrBfSvuY8W+FYztMZaPrv6Itn5tkWV47z1xQhsxQmSd\nNUbM316rleE7dhDh7c3Snj0xtNTUw6PxeIQomDQJTrPksizLuMpdOIudtTa3yS0281Hb0T9b3MhO\nGdlVxybLSFqpzs3L4IXWX4vGT3Nk86/9XBesQxd6ZNMGa/HSnea9yq23ClGwdi0MG3bm/9sm4oGO\nHcl3OHg4IwNfjYbpHTooOr5WKzJzunWDBx6AtDSYPx+6Bnflz5v/5L2N7/HE70/wQ9oPvH3F24zt\nMbZl1uFQaVZUUaDSZFS6Kpm9bjaz1s2irV9bfpz0I6NjRyNJEjYbTJ8OX30lRMHs2eIkpzQHKisZ\nvmMHQVotq3r3JlinU36SxmD1auFLufVWPHYP9jw79lw7jlwH9hzx3J5rx1HgqLn4u0pcyK7j76K9\nDF4nvkj7adCF6dAYNUj6Oi74uqoLiZtaQsHj9Bz52SnjqfTUCA5nifN4EWJxwwlu8rVBWiEQQrTo\nw/ToI/R4d/TGO6Jqq3quueACpJgYEXB4FogCgJejo7G63dy9Zw8+Xl7c1K6d4nPcey907y7qPA0c\nCD/+CAkJGh4c/CCjY0dz/4r7Gf/teK7seiXvXfkeXYK7KG6DytmL1BLzWSVJ6gds2bJlC/0aY+1Y\npcn5Zd8v3PvzvWSXZ/PIkEd4etjTGPVGQBRimTABUlPhs8/g+usbx4Zcu51h27YBsKZvXyK8vRtn\nogbicXioPFCJLcNGZWYltkwblQv+orJUh92/C86i2ultXkavmgulvp0eXZgOXUjtO/CaLUSHl75l\neA1lj4zr8PErGc6SI88dBQ4ceUL4nPDv9jFjKE3DcNNl+MQFYYgx4BPjgyHGgC6oZQo+WZa5c88e\nvjh0iIXx8UwID2+UeTIyYOxY8f36/PMjgbqyLJOUnsT9v9xPgbmAJy58gv9e+F8M2hZcm0OlXmzd\nupX+oo5Jf1mWt9bnWHWlQKVR2VOyh//+9l9+TPuRS6MvZfn1y+kR2qPm90uWwG23QWioyDLo06dx\n7Ch0OLhsxw6cssyaPn2aXRB4XB4qsyqxplnFlm6lMkMIAPtBe80dtKSVMHTUYMgrxf/CroSOjDj+\njjlAc1YuA0teknAfBOug+6n3P+EKyd4yKudsxbTyIIVLynGbjqT9adto8ekiBIJvd198e4jNJ9YH\nrV/znfokSeLj7t2xud1cn5qKwcurUVJhu3SB9etFGu/kySJg9623wGCQGNNjDJfFXMbMtTOZuXYm\nX+38itmXzWZc3Liz8rOkohzqSoFKo1BoKeTFv17kky2f0MG/A68Mf4XJPSfXnHBsNvjPf+D99+G6\n68QKQWBg49hS6nRy6fbtFDidrOnTh25N2AbZbXNjTbFiSbVgTbXWiADbXhuyU3z3NH4afGJ98Onq\nU3OX69NFPPfu6I309JPw0UeQkyMC7FRqc/PNsGYN8t69uMplsbJSvcKSKVZcrOlWHHmOmkO8O3rX\niITqzdjLiD5M32RmuzweJqeksKykhOWJiVweHNwo88gyfPqpKPoVFyfiDOLijvw+rTiNh359iF/2\n/cKQjkN4fcTrnN/p/EaxRaVpaMhKgSoKVBTF6rTy1vq3mP33bLwkL54a+hT3Dbqv1tLk9u1www1i\nefONN+CeexqvOF+p08nInTvJstn4q29fEozGRplHlmXs2XbMO81YdlpqHq17rOAR++gj9OLiE2es\ndTHSd9Cf/O7MZhPdEKdOFbd5KsezefMR5/mYMSfdzVXhwpp+RJidSKDp2+kx9jLi18uv5tG3hy9e\n3o3jcnF4PIxNTuaPw4dJ6tmTyxpJGADs2CFcCPv3w8svi2DEo2Nsf8v8jcdWPcb2/O2MixvHrOGz\n6B5yGks4Ki0OVRSoNDsOt4P/2/5/vPDXCxRaCpkxcAZPD3uaEN+Qmn1cLiECnnlGtIL9+mvo2bPx\nbCpwOLh8xw4OORys6tWLPv7+iozrtrgx76p98TfvNOMuF0vX2iBtrQuLMdGIMcGI1v8MlqznzhVR\n9nv3in4HKidmyBCRl3cG3RM9Lg+VGZXHvaeVWZWAcOFUryTUiIU+fni3V8YFZXO7Gb97N6vLylic\nkMA1jVhV02aDp56Ct9+GCy8UMZpdjooz9Mgevtn5DU+tfopD5kPc2e9Onhz6JBEBZ0FhL5UaVFGg\n0mw43A6+3PYlL697mYPlB5mYMJGZl848LqJ5yxaYNg127oRHHoH//Q8a062fU1nJZTt2UOF281vv\n3sSf4QqB2+LGvN2MaYupZrOmVt39a8A31rfWXaWxl1Es+Sux9CHLMGAAtG0LP//c8PFaM/Pni+Wn\n3buF4lQAV4ULS7LluNWf6rgFfXs9/v398evvh39/f/z7++Pd4cw+1HaPhylVroT5cXGNFnxYzV9/\nCa1ZWAivvQZ33VV7ta7SVcm7G97llXWvYHFamNZ3Go9f+DidAjs1ql0qyqCKApUmx+6y88W2L5i1\nbhY5FTlMTJjIM8OeISE8odZ+ZjM895y4M0lMFL7NgQMb17Z0q5WRO3YgA7/37k3X04whqEsASHoJ\nv95+Ry4C/fzxjfdFY2jEGgfr1sHQoUIQXHll483TGnA4oHNn4T74+ONGm0aWZSr3Vx75nGw2Yd5i\nxlksMiMaIhScHg83p6WxqLCQj7t35w6F6xgci8kEjz0Gn3wCl10GH34oahwcTYW9gvc3vs8b69/A\n7DBze9/beeLCJ1Rx0MJRRYFKk2Gym/h82+e8sf4NcitymdxzMk8Pe5r4sOPvzlasENXVCgvhhRfg\nwQehscsCbKio4OqdOwnX6/m1Vy86naQFsuyWsaRaqPi3gop/KzBtMGFJsRwRAL388B9w5ORuTDA2\nfSrflVfCwYNiecWrZaQRtmhmzYLnnxc5eE3Yx6I6nuRoMWnabMJV4gKqhMJAfwIGBxAwKAD/gf4n\ndSW5ZZkH9u7lg7w8XoyK4unOnRs9G+DXX+HOOyE/Hx59FJ58UnhijsZkN9WIgwp7Bbf0uYWHhzxc\nK5NIpeWgigKVRienIod3N7zLnC1zsDgtTO45mScvfJK4sLjj9s3IgMcfFw39LrtM3Lh1aYL6KD+X\nlDBh9276+PmxLDGxVmEiR6GDig0VR0TARhNusxu8wNjTWHOy9u/vj7FnMwiAY9m4UTR/WLiwcTpB\ntUYqKiAqSgRlvvNOs5pSSyhsMlGxsQLTJpNwPUhgTDDiP+iIUDDGG5E0Us2xL2dn83RWFnd36MC7\nXbuibWRRaLWKgmGzZ0NYGLz+OkyceHwAsMlu4sNNH/L2hrfJN+dzTfdreGTII1zU+SI1lbEFoYoC\nlUZj26FtvLH+DRbtXoRRZ+TO/ndy/6D76RjQ8bh9S0rgpZfggw8gPBxeeUW4eZviXPFeTg4P7tvH\nNSEhfNO1B55kW40AqPi3gspMETSmC9cRMCTgyF3bgJPftTUro0bBvn2QnFw7RFylbl58UawYZGZC\n+/bNbU0tZLeMNc0qPpMbxGZJFqtTGj+NWE0YJD6b/oP8+VouZXp6OiOCg1kYH09AY5T4PIbMTHjo\nIUhKEl0X33vvxMHAdped+bvm8+a/b5JcmEy/9v14ZMgjTIifgE7TMotGnUuookBFUSpdlXy3+zs+\n3vIx/xz8h86BnXlw8IPc3vd2/L2Pj+CvrBQnj5kzRYn+xx8XroKmKAfgdLt5cm06m/8o5MaDfvRJ\n88K81YRsl5H0Ev79/I/ckQ0OwNDZ0PLvaLZuhf79RXrGDTc0tzVnF4cPi9iCadNEqksLx2V2Ydps\nwrTBVLOS5Tgk6il4R3pj62dgbkQFh3vreXt8L2KCGyel9lh++UXUNcjIEB/B55478WqfLMuszFjJ\nG+vfYFXmKtr7tWdav2nc0e8ONe6gGVFFgYoi7CnZwyebP2HujrmU2kq5LOYy7up/F2N6jEHrdfxd\nisMhrlsvvijq6tx1Fzz7rFglaCzcFjemzeIEWvzPYXL/KcO/SHyGDdGGmjutgMEB+PXxa7T88kZl\n7FgRRZ+S0jgNIFo7zz4rBEFWVuN+GBsBWZax59hrBIJpg4nyLSaweXBpQNPLh44XBIvP+JAADNGN\nJ3LtdlFUbOZMERd0660inTgy8sT77yrYxUebP+KrnV9hdVq5utvV3D3gbkZ0GYHGS13takpUUaBy\nxpgdZn5I/YG5O+ayOms1wT7B3NrnVqb3n063kG4nPMZmE7XUX31VxMGNGyeKocTGKmub7JGx7a3t\nBjDvMoMbMHqREiuTliAx9opILhjeAX3bpqtG12hs2SLSEOfOFZX6VOpPaalYLbjzzrNiteBUeJwe\ncrce5v2le9FssTFsnw5Dpsh20IXpakRwwOC6gxjPFJtNFNScNUuEbdx5JzzxBJwsOcJkNzF/13w+\n2vwROwp2EBUUxS29b2Fq76nEtIlR1DaVE6OKApV64ZE9/Ln/T+btmMfilMVYnBaGdR7GHf3u4Lr4\n607aGKWiQpwc3nwTiothyhRxckhIOOHu9cZZ5sS00XREBGyowFUmIrh9431r4gB+62Lnbq9s4gP8\n+D4h4aQZBmcdsixSEA8fFmUf1VWCM+ell8QSVnKyaBnYCnB5PDyVlcWrBw8yVRfCzLK2uDZZar4r\n7vKqIMaexlpCwbeHL5JXw1cTzGbhJnztNfH8+utF/EHv3ifeX5ZlNuZu5JMtn/BdyneYHWaGRg7l\npt43MSF+AoGGRqprrqKKApVTI8sy2/O38+3ub/lm1zccrDhI1+Cu3NTrJqb2nkpUUNRJj923T2QQ\nfP65iFK+5RbRt6AhGQUelwfrbmutVQBrmhUAbbD2yEltUAD+5/mjC9JR7nIxPT2dRUVF3Nm+Pe90\n7YqhNQXhVRfg+e03GD68ua05u7HZRBGj+Hj46afmtkZRFhcWcnt6OqE6HYvi4xkQEIDsOSqIsWqz\nJFtABk2AppZbLWBQALqQMw8GLC8XboV33hErhZddJtqdjxx58sxZi8PCj2k/Mm/nPFZlrMJb6821\nPa5lcsJkRnYdqXZoVBhVFKickGoh8F3Kd3y7+1syyjII9glmQvwEbu59M4M7Dj6pP9LtFjVzPvhA\n5DEHB4tuhg8+eGYp4PZ8uwimqj5pbarAY/GABvx6+9W6s/Hp6nOcXevLy7kxNZVip5M5sbFMOst8\nxafEYhH+l/POg++/b25rWgdLlohuWz/9BFdd1dzWKEqmzcbklBS2m83MjI7m4U6d0BzznXGZXCId\n8iihUN1+2qebz5EsnMEBGBONeGnrF3/jdIqP6htvwKZNosnSjBliBaFNm5Mfl1uRyze7vuHrnV+z\nq3AX/np/RsWOYmL8RFUgKIQqClRqcHvcrM9Zz7L0ZSxJXVIjBMb2GMvEhIlcEnVJnSlDOTnw1Vei\nytmBA8K9PWOGSJX38Tk9Gzx2D+bt5topgftFSqC+vb7Wyci/vz8a35Pf7ds9Hp7fv59Xs7MZ6O/P\n/Ph4Yk7XkLOJZ54R67KpqRAd3dzWtA5kWay45ObCrl2gbwUxJ0fh8Hh4OiuL1w8e5ILAQOb26EGX\nOr4bsixTmVVZO0ZnmxnZJePl64X/AP9a4vx0ezvIMvz9t3ArJiWJAmXjxsHtt4u0xrpKLKQVp/Hd\n7u/4NuVbkguT8df7c033axgdO5orul5BkCGonv8VFVBFwTlPhb2CX/f9yrI9y/h578+U2EoIN4Yz\nqvuo0xICFRXipurrr+GPP0RPgsmTRffCU5Uklj0y1j1WkVa1yYRpownTVhOyQ0bylvDvf8yJph59\nATZVVHB7ejppVisvREXxWKdOjV7EpVlIS4M+fUQ5uZdeam5rWhe7don/7f/+J0r1tULWHD7MLWlp\nFDgcvBITwz0REcetGpwMt82NeWttAW/PsQMiJfLo765fX79TlvXOz4d584Srcc8eiIkRWQs33HBq\nrVstEJakLmFHwQ60XlqGdR7GqO6jGNV91HH9VFROjioKzjE8soft+dtZmbGSlRkrWZe9DqfHSWJ4\novgCxY7ivIjz8JJOfgG12WDVKuHGXrpUpB9dcgnceKNQ+YEniAGqrvtu2mQ6IgK2mGoaxPh09REF\nWKpWAvx6+51RZcAKl4uns7J4PzeX3n5+zO3Rg95+fvUe56zA4RAd/sxmUZ+gkVo7n9M88YQo0bdh\nA7TS84nZ5eI/mZl8lJfHef7+zImNPePvTGVOZS1Xn2mzCU+lB0krYUw04j/AX2z9/YXb4QTf8erV\ng88/h2+/FbFIAweKKokTJojkkLrILs9m+Z7lLNuzjNVZq3G4HcSGxDKyy0hGdBnBRVEX4advpecE\nBVBFQStHlmX2H97Pn/v/ZFXmKlZlrqLYWoxRZ+SS6EsYETOCUbGj6gwWBCgrE+7VH38UxUksFlGt\nbOpUkUnQ6ZhaI/Y8u/BJbhInhqPruXtHeosTw8AjJwhdm4ZVMpNlmYWFhTyWkUGZy8WL0dE8EBHR\nOlcHqnnqKZHbuX698NWoKI/DAYMHCyW8ZUvTVNVqJv4uL+fO9HTSrVbu79iRZzt3JqiBDUc8Tg+W\nnRYqNlTU9HWw7LaAu3afkOpeIcYEI166I99Zs1mcd779VsQpVVaKCt4TJsDo0cc3YToWs8PMqoxV\nrNi3gpUZKzlQfgCdl44LIi/g8pjLGR49nH7t+6mVFI9CFQWtDI/sIbUolbXZa1lzYA1rDqwh15SL\nhMSADgMY0WUEI7qMYHDHweg1J/eTyrJYwlu5Uvj6/vwTXC4Ry3bttWKLizuqTvs2E+ZtZszbzJg2\nm2oqq+nCdWIFYGBAzZdf6ZoAW0wmHti7l78rKhgTEsLbXbsS1RpjB45m7Vq46CLhMmilS9sthtRU\nsUpw++3w/vvNbU2j4vB4eOPgQWYeOICPRsPM6Ghub9/+tF0Kp4Pb5sa8w1xzs2DafFRHUW8Jvz5+\nNTcL/gP88e3hi5fOC5MJli8XAmHFCrFC2a0bXHMNXH21yMitK/RDlmX2le4Tq6SZK1mdtRqzw4yv\nzpchHYcwrPMwhnUexqCIQfjoWvn5ow5UUXCWY3fZ2VGwg3XZ61ibvZa1B9ZSYitB66Wlf/v+DOs8\njKGRQ7kg8gKCfYLrHKu0FH7/XQiBlSshO1t8yYYNE4XyRl/lIchkxbxdXPzN28VWXQ9AF6rDr69f\nrVWA+sQB1Je9VivP7t/PwsJCehqNvN21K8PrCl1uLRQVifXUjh1Fc/vWlFrZUvngA7j3XvjhB6GI\nWzm5djtPZGbyVUEBCb6+vBQdzZjQ0Eb7Lte0Hj9aKKRbQRZCwZhgxK+PH369/fDr44fUxchfW3T8\n9JNYScjNBX9/keI4fLhwZ8bF1d07xel2si1/W83N09rstRyuPIzOS8fAiIEMjRzK0MihDO44mBDf\nkEb5u1siqig4i3B73KQVp7EpbxObcjexMW8jOwt24nA7MGgNDO44mGGRQu0O7jgYo75uH3NRkbjh\nXLsW1qwRNW88HvFlunqok+ExFrprzTjTLEIA7DIj26vKAscY8OsrvqD+ff3x6+OHvoO+SXoDZNls\nzMrO5otDh2jv7c1znTtzS7t2rdtVUI3NJs56GRmiG+KpHKwqyiDLwqn9009CiJ0qiraVsKmigiez\nsvitrIxB/v48HxXFyODgJvmeu0yuIzcfO8SjJdmC7Kg6B0UZMPY24tfbj+IgP9bk+rF0g4F/N0g4\nndC2rRAHl14qMhm6dq1bJHhkD7sLdwuRkC2EQr45H4DooGgGRgzkvA7nMTBiIP3a92u1cQmqKGih\nON1O9pTsYWfBTrYc2sKmvE1sPbQVs8OMhESP0B4MjBjIwA5i69OuD97ak6cBuVxiFXTTJhEztWaN\nCFw34GJIeyuXRllI9LPQzmbBvc+CI18s/0s6Cd9435oLv19foda1gU1fMS/NYmFWdjbfFBQQrNPx\n38hI7unQAZ9z5U7Z4xH5nefYhanFYLOJK0xmJvz77zmV/rm6rIynsrL4t6KC/n5+PNm5M2NCQxV1\nK5wOHqcHa7oVyw5LLbFQXUNB46/B0MMXU7CRPU4j63KNrNpjpFDWExoqcf75cP75Ij53wIC6Q0Rk\nWSazLJONuRvFjVjVOdjqtOIleREXGsfAiIEMaD+A3u16kxie2CoqLaqioJmRZZl8cz47C3ayq3AX\nOwt2srNgJ6nFqTjc4sLcObBzLQHQv0N/ArwD6hhTVBLctOnIlr7VRajNSiRW+oVY6elnoa3VgrZI\n1ABAEhkAxp7GWptPN59agT9NjSzLrCor4+2cHFaUlhKh1/NYZCR3tG+P77kiBkC8qY8+Cm+9JZaw\nx4xpbovOTYqKxBVFpxNLbKGhzW1RkyHLMn8cPsxLBw7wx+HDRBkM3NWhA7e1a0dYM9ZxkGUZR74D\n8w4zlh0WLLvFZk214rF5xD5GDeVtjGR6jGwuNpLu8CVHY6RjLz19+0n06yfCRnr1qlsouDwuUopS\n2JS7qUYo7CzYicsjXKidAzvTq20vEsMT6dW2F73a9qJbSLcTNoVrqaiioIlweVxklWWRXpJOenE6\nacVppJekk1qcSrG1GAA/vR89w3vSK1x8mBLbJpIYnkgbn5P7yc1m0RBv1y5I3imTvamSwzutBFus\ndMJKN4ONSKwYKx01x3h39D7u4u8b51tnIaCmJtdu56v8fL7Mz2ePzUYfPz8e7NiRyeHheJ8LboKj\ncbtFL9oPPxQF5O+9t7ktOrfZswcuvBBCQkTJzpO1/mvFbKqo4IPcXBYWFgIwMTycGRERnOfv32La\ni8tukQZdLRIsuy1Yd1uxpFpq3KAOvYZDWl/22nzJln3IlXzRRvnQYZAviQM19Osnsqzq0n4Ot4P0\n4vSaG7rqm7tcUy4A3hpveoT2IDY0lh4h4jE2JJbY0NgW6YJQRYGCONwOssuzySrLIutwFpllmTUi\nYF/pPpwescRl1BlrPhg9QnvUKMqooKiT1gewWMTdf8ouD/vW2ynYasOytxJ9iY0O2OiEjU6SFV3V\neyLrvfDp5oN/gi++PXzxjRWbT3cfxTuhKUWl283SkhLm5uezsrQUvZcX40NDuaNDB4YFBraYk02T\nYreLvM8lS0SpyGnTmtsiFRDCYMQIIdhWrhSBOOcgxQ4HX+bn81FeHlmVlfT18+P68HAmhocT2UKb\njcluGVumTYiENCu2PTYsaVbMaTbkMmfNfsXoycaXHHwoMwqxEJzgQ6fzDPTorSE+Htq3P3mcQqmt\nlF0Fu9hRsIOUopSaa8Eh86GafSL8I2quBd2CuxHTJoaYNjFEt4luNsGgioJ64JE95Jnyai76NY9V\nz3MqcpAR/xMvyYtOAZ2OqMIqARAbGkuEf8QJL3B2u4gf27vdRc7mSkp22ajMtOFVUEmgRVz822JH\nWzWHRwJ3qAFDtIGQ3r4E9Ky6+PfwxbuTtyLdzRobi9vNb2VlJBUX831xMYddLoYEBHBru3ZMDA8n\n8Fzu9peTI8q5bdgACxeeE1HvZxV5eaKTT16eKMV39dXNbVGz4ZZlfikt5fNDh/i5pAS7LDMkIIBJ\n4eFMCAujg/fplT1ubpylTmx7bVj3WLGk2SjcasWSakXKs6Fxemr2K0HPIQwU6wx4wn3QdjIQEGug\nbR8fOg/wplusRGjoiQVDeWU5e0r21KwWp5eIlePMskysTmvNfmG+YTUiIaZNDNFB0UQGRtIxoCOd\nAjs1mmhQRUEVNqeNXFMuORU55FbkkmvKPfJY9fyQ+VCN7wgg3Bhe82ZFB0UT3ebIY6eATscVxHA6\nITvTw4EtdvJ32ClLr8R2wI6cX4nusJ2ASjvhVOKHu+YYh1aDPdiAppMP/rE+tOtnICTRB58uPnhH\nejerv/9MOWS3s7ykhKSSEn4rK6PS46GHry9jQ0O5uV07YltxgZjTZvFi0Xzex0cIgqFDm9silRNR\nViZKef78s2j08dprp9/oo5VS4XKRVFzMoqIifi0txSXLnB8QwIjgYEa0acMAf/+zLlNI9sjY8+xU\nZlVi2VtJ/jYbJcmVVGZVoimyYbQecc+6kCjAmyKNAUegASncG59Ib/xivAmN86ZDH28i47WEhEi1\nRIMsyxRaCsksyzyyHT7yPLcit+amEyDQO7BGIHT0r3oM6CheC+hEREAE/vr6u3NatSiI7RlLoaXw\nxJtVPOab88mtyKWssqzWOAHeAUT4RxARECEeq55HBkYSHRRNVFBUrZS/inKZg8lOClIclKTbMe13\nYMt14C6041XmwGCxE+SwE4yDo78OVq0Wq78BT6g3+o4GgmK9ad/Lm7a9Dfh29UEXpjvpm7pgwQKm\nTJnSCP9FZcmz21lbXs6aw4dZU15OssWCF3BhYCCjQ0MZFRJC9zqEwNnydzaUBQsWMOX880Uxovnz\nYfx44TIIaV050q3u/ZRlUcfgsccgKkp097niChYsXNi6/s6TUNf7edjp5MfiYpaVlPB7WRnlbjeB\nGg2XtmnDiDZtuDgoiO6+vnidJa7Bk/2t7ko39gN2ylJs5GyqpHR3JbasSjwFdvTllfjZjznvo6FY\n8sbs440jyBsp3BtDR28CovSEddcTHqenfbye0PZeNcLB4XaQWyFuXA9WHBSP5QfJMVU9VuRQYCmo\nZZdBa6CtsS1t/drSzq+deH7sz37itQDvACRJanpRIEnSDOBRoB2wA7hPluVNdex/MfAGkABkAzNl\nWf6/OvbvB2wx3GOgMryy9u+QCPENIdwYfmTzDT9y4T/q0ajzo6LUQ0G6k+J9Dg7vd2LcuTrUAAAN\nuUlEQVTKcWLNc+IscuAudeJ12IG32YGf3U6Q7KxZ1q/GpNFh89HjCtBDqB5DpIHAbt60TTTQsb83\n/l0MaIxnHtw3evRokpKSzvj4xsDscrHTYmGb2cwWk4m15eXss9kA6O7jw7CgIC4KDOTKkBBCTrOE\nakv8OxWnuJjRgwaRlJMjese+8grcfHPdidVnKa32/dy9G+6+W2QlXHIJo51OktaubW6rGp3TfT9d\nHg+bTCZWlZWxsrSUfysqcAOBGg3nBQQwKCCAQf7+DAoIaNZshro408+ux+WhMtdB/k47BbvslKbb\nse6vxJlnR1Nix8dsx89ZWzgAVKClQqsX1xF/PQTr0bbV4xOhJ6CjjjbROsJidIR11dGmowaX7KwR\nDrmmXArMBRRYCigwF5Bvya/1c3WMWzV6jZ4w3zB8i33ZO2svnIEoqLezV5KkSYgL/J3ARuAh4FdJ\nkrrLslx8gv2jgOXAh8D1wGXAZ5Ik5cmyvKquue457x769+tPmG84/s4wtMUhUBiA5ZCMOc+JtdBF\nZbELZ6kL92EXNpOTTLOLbIsFH3syfi4nxqOW8X2rNjtemLU6Kg068SZ1NeJoF4ylk56gLnpCY73p\nkKgnMFp/Rg19zhZsbjf7bDb22GzssVrZYTazzWxmr82GDOgkiQSjkSuCgxkWGMjQwEDanSV+xSZD\nlmHdOrEasHixKCbx3HPw0EPQWps4tWYSEkT9iJ9+gscfFyJh0CCYPl3UlzjHG1ZpvbwYEhjIkMBA\nno2KosLlYkNFBRsqKthoMjEnL4+XnOJC1UGvJ95oJMHXlwSjsWY7W2OMvLRe+HY2ENPZQMyoE+8j\nu2XMuU5ydzkoTrdzONOBLceBfMiBvtiBodyO9z4Tfrsd+Mri2uQG8qs2JxJmLx02nRa7wQ+XsSe+\n/n3pEqyjR4gOn7Y6jB10+HfQ4hejQRNuxR5cjFlXSKG1gCJrEUWWIlJ2prCXvWf0d57Ju/MQ8Iks\ny/MAJEm6C7gauA149QT73w1kyrL8n6qf0yVJurBqnDpFQci03vh4opA9TiopA464B4xVmwMJq5eO\nSp0Wh7cWt48WT4QvtuBAPOE6XB10+HXU0yZaR2hXHW1j9XgHtpy0vcbE4naTY7dzsLKSg3a7eG63\nk2Gzsddm46DdXrNvkFZLYpUAeNzPj75+fsQZjede6uDpYDKJO8mffhKF3LOzRam1//1P1Jh+5pnm\ntlClIUiSKMZ/5ZWinkFwsMgYue8+UYP3mmtE1kJkZKtcBaoPAVotlwcHc3mwKL8uyzIHKivZYDKx\ny2wmxWplRWkp7+XmUh3iF6LVEu3jQ5TBQHTVFmUw0MHbm3Z6PWE63VnjijgWSSPhH6mnR6Qerq77\npsBldVOc4aRgr4uSTCfmXCfWfCf2IieuEidyuRONyYn+YCXeGS78XE58OCIkyo8ay4mEt9SRUG0U\nATotsjYOWHxGf0O9RIEkSTqgP/By9WuyLMuSJP0GDDnJYYOB34557VfgrVPNZ+8aSEWHDlQGazGE\nafEN1+LfXkdARy1tIrUER2rxaeUXeI8sY3a7KXe5KHe5qDjqeXnV8zKXiyKHgyKnk0Kns+Z5hdtd\na6xwnY6O3t7E+PgwNSCAbj4+dPf1pbuPDyG6k8c9nJPIsrj4HzoEe/eK9LXkZFGWOCVF/D46WhQg\nuvZaUYPVy0uIBZXWgUYD7dqJbmJZWaKLz08/wV13icqUbduKVYQ+fSA2Frp3F0IhJOSc7WUhSRJR\nPj5E+fgwKTy85vVKt5t0m40Ui4Wsykqx2WxsNZnItttxHeXG1gBt9Xra6/W01etpo9USrNPRRqsV\nm05H8FHPAzQajBoNvl5eeHt5nTXnMa2vhnaJGtolnv4xdrOHogwnZdkuTHkuzAWumhVzR5lYMfdU\nuKBUAxVnaFc99w9FvGcFx7xeAMSe5Jh2J9k/QJIkb1mW7Sc4xgDgGLcRS2Q+5qoXZQBZRj4I8kHg\nb/GafPTvOfLa0dEB8olel+Xjjj/lWEcd467aXJIknsuy+BnwVL9W9bP7qM0B2Ku2Xfv30/fjj2t+\nrtkkqWa/k96NyDJ+gL8s0wYIlmXCZJnuskwbWSZUlmkry4R7PITLMidc+JdlshGBHse+frI5640s\nU75vH1tfe+30x6trnvraIMtiWd9uF2107fbjn9vtcPgwlJSIaPTSUvH7agwGEYAWHy+CB3v1Ej9X\nvzfbtwNQXl7O1q31cuGdlZyTf+fll4utvFy838nJYvvgA/G5qUaSIChIrDBUb0aj+AwdvXl7C/Gg\n0YhjNBohLL28aj+v/lmSTm9l4gz2Kd+/n60ff1yP/8yZzRfL8RcKF1B8kq0UOASYEHfGJqCSkyMB\nPkdthqpHHaCvekw5cICRc+agq/q5etMe9eh1ik1TNdepHo+1ra7Ho/c7nX2QQGoHtDt+P1d2NsyE\nqn9BvahXoKEkSe2BXGCILMsbjnp9NjBMluXjVgskSUoHvpBlefZRr12JiDPwPZEokCTpeuCb+vwh\nKioqKioqKrW4QZbl+fU5oL4rBcWIm922x7zeFhEncSLyT7J/xUlWCUC4F24A9lO3MFRRUVFRUVGp\njQGIQlxL60W9RIEsy05JkrYAw4EkAEk4cIYD757ksPXAlce8NqLq9ZPNUwLUS92oqKioqKio1PDP\nmRx0JqHlbwJ3SJJ0kyRJPYCPEZl+cwEkSZolSdLRNQg+BmIkSZotSVKsJEn3ANdVjaOioqKioqLS\nQqh3SqIsy99KkhQKvIhwA2wHRsqyXFS1Szug01H775ck6WpEtsH9QA5wuyzLx2YkqKioqKioqDQj\nLbLMsYqKioqKikrTo1amUVFRUVFRUQFUUaCioqKioqJSxVkjCiRJ0kuStF2SJI8kSb2a2x6lkSRp\nqSRJByRJskmSlCdJ0ryquhCtBkmSOkuS9JkkSZmSJFklSdorSdLzVZUyWxWSJD0pSdLfkiRZJEkq\nbW57lESSpBmSJGVVfVb/lSRpYHPbpCSSJA2VJClJkqTcqvPN6Oa2qTGQJOkJSZI2SpJUIUlSgSRJ\nP0iS1L257VIaSZLukiRphyRJ5VXbP5IkXdHcdv1/e3cXKlUVhnH8/3gQCY1EwUqLyERCCjXKbrLo\nwrvMSjKjKAwSiT4gpA8iipCQ6Isyoy7SgqCCiuoi+jpKSJRRghcaIQopBcEppLSL0reLtU6MJz1n\nZtxz1uzt84NhmD2zh3ex9+z9zl5rr7fXJD2U99+OBvXXJikg1VU4ADR1EMQgcCMwF7gBuAB4t2hE\n1buQNOnWncA8Uv2LNQzPvdUsE4F3gJdLB1KlloJojwELSVVSP8mDj5tiMmkA9V0093gDsBh4Ebic\nVKhuIvCppNOKRlW9/cCDwCWkafoHgQ8lzSsaVQ/lRH016ffZ2bp1GGiYZ0B8GlgO7AIWRMTOslH1\nlqSlwPvApIg4Mtbn60rSWmBNRMwpHUsvSLodeC4ippWOpQqSvga+iYj78muRDrovRMTxCqLVmqSj\nwHUR0cA60cfKid2vpNlpt5WOp5ckDQFrI2JT6ViqJmkK8B2pGOGjwI6IuL/d9fv+SoGkM4FXgVuB\nvwqHMy4kTSPN6LilyQlBNpU0xbn1uZaCaF8ML4v0r2K0gmhWH1NJV0Ya+3uUNEHSSmAS0NTqZS8B\nH0XEYDcr931SAGwCNkbEjtKB9Jqk9ZL+JE0nfT5wU+GQekrSHOBu0gRX1v9GK4h21viHY1XJV3ye\nB7ZFxK7S8VRN0kWS/iDVmHsFWBERewqHVbmc8CwAHu72O4okBXnWw6OjPI5ImivpXmAKMFxMqR41\nMbN229myylOkDbqEtPN+UCTwDnXRTiTNAj4G3o6I18pE3plu2mlWExtJ43xWlg6kR34A5gOLgA3A\nW5IWlg2pWpLOISV2t0TE311/T4kxBZKmA9PH+Ng+0kCta0YsHyBV23wzIlb1ILzKtNnOvRHxz3HW\nnUXqq70iIrqaw3q8dNpOSTOBLcBX/b4NW3WzPZs0piB3HxwGlrf2sUvaDJwREdeXiq1XToUxBZI2\nAEuBxRHxvyrqTSTpM2BfRKwuHUtVJC0D3iMVLRz+Az1A6hI6QhqfNuYJv+NpjquQCx4NjfU5SfcA\nj7Qsmkmq+rQC2N6b6KrTbjtPYGDEc9/qpJ052RkEvgXu6GVcVTvJ7Vl7XRZEsz6WE4JlwFWnSkKQ\nTaAGx9YOfQ5cPGLZZmA3sL6dhAAKJQXtiogDra8lHSJlQHsj4ucyUVVP0iLgMmAb8Dswh1Rb4kdG\nqSZZN/kKwVbSVaAHgBnpnAIRMbKfutYknQtMA84DBiTNz2/tiYhD5SI7ac8Cm3NysJ10W+l/BdGa\nQNJk0m9w+N/W7Lz9fouI/eUiq5akjcDNwLXAoTyoG+BgRDSmZL2kJ0ldlT8Bp5MGcV8JrCsZV9Xy\nceWY8SD5nDkUEbvb/Z6+TgpOoP/voezcYdLcBI+T7pH+hbQTrzte10KNLQFm58fwwVWkbdq0rP0J\n4LaW19/n56uBL8c/nGq0URCtCS4ldW9FfjyTl79Oza5ujWENqX1bRyxfBbwx7tH0zgzStjsbOAjs\nJO2zW4pGNT46Pl/WYp4CMzMz67063JJoZmZm48BJgZmZmQFOCszMzCxzUmBmZmaAkwIzMzPLnBSY\nmZkZ4KTAzMzMMicFZmZmBjgpMDMzs8xJgZmZmQFOCszMzCz7F30nyJBjDCnRAAAAAElFTkSuQmCC\n", + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot the primitive gaussians, the overall GTO, and the STO it mimics\n", + "\n", + "nplot = 161\n", + "delta_plot = .05\n", + "rvals = np.zeros(nplot)\n", + "yvals = np.zeros(nplot)\n", + "indvals = np.zeros((3,nplot))\n", + "stovals = np.zeros(nplot)\n", + "for i in range(nplot):\n", + " rval = i*delta_plot - 4.0\n", + " y = cg_unroll.subs(r,rval).evalf()\n", + " rvals[i] = rval\n", + " yvals[i] = y\n", + " stovals[i] = 2.0*math.exp(-abs(rval))\n", + " for j in range(3):\n", + " val = h_coeff[j]*gto000.subs({alpha:h_alpha[j], r:rval})*radial_norm.subs({alpha:h_alpha[j]}).evalf()\n", + " #print(val)\n", + " indvals[j,i] = val\n", + " \n", + "plt.plot(rvals, yvals, label='H GTO')\n", + "plt.plot(rvals, stovals, label='H STO')\n", + "plt.plot(rvals, indvals[0,:], label= r'$\\alpha=%.2f$'%h_alpha[0])\n", + "plt.plot(rvals, indvals[1,:], label= r'$\\alpha=%.2f$'%h_alpha[1])\n", + "plt.plot(rvals, indvals[2,:], label= r'$\\alpha=%.2f$'%h_alpha[2])\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAABkAAAAOBAMAAAAoFKpzAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIpm7MhCriUTv3c12\nVGZoascqAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAoklEQVQIHU3MPw4BURDH8e8mXsjaQq0SCoVG\nIlG/G3iV1jYaDdEouINGQycOwA1WFBLdxgW2VkkUiGb9XjSmmJlP5g9BvWNRmOZJtU/x7VVOzRaG\nsPa6wRw2MLZSFzJH4n56wjL1S3unJw/pKJReSoFSb6ImbP+rIWA085tR7IXuMn25ElSkHdwdhRqR\n1wVaUJ1NRwwOhKlZQZLnH8IYszhbvj+3KLWMktYcAAAAAElFTkSuQmCC\n", + "text/latex": [ + "$$2.0$$" + ], + "text/plain": [ + "2.0" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Compute normalization for STO for previous plot\n", + "sto00 = exp(-r)\n", + "val = integrate(sto00*sto00*r*r,(r,0,oo))\n", + "1/math.sqrt(val)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 2", + "language": "python", + "name": "python2" + }, + "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": 1 +}