diff --git a/notebooks/gto_integrals.ipynb b/notebooks/gto_integrals.ipynb index 38d7894..473e58b 100644 --- a/notebooks/gto_integrals.ipynb +++ b/notebooks/gto_integrals.ipynb @@ -6,268 +6,451 @@ "source": [ "*Copyright (c) 2023 Graphcore Ltd. All rights reserved.*\n", "\n", + "$$\n", + "\\def\\br{\\mathbf{r}}\n", + "\\def\\bA{\\mathbf{A}}\n", + "$$\n", + "\n", + "# Integrals over Gaussian Type Orbitals (GTO)\n", + "\n", + "We have a molecule comprising $M$ atoms (Mnemonic: $M$ for $\\text{ato}M$),\n", + "each with atomic number $Z_m$ and position $\\br_m \\in \\mathbb{R}^3$\n", + "$$\n", + "\\mathcal{A} = \\{ (Z_m, \\br_m) \\}_{m=1}^{M}\n", + "$$\n", + "Each atom in the molecule contributes both negatively charged electrons that are bound \n", + "to a positively charged nucleus. Both the number of electrons and the charge on the\n", + "nucleus can be found from the periodic table. For a concrete example, a single Oxygen\n", + "atom has atomic number 8 and this implies a charge of $Z = 8$ located at the position\n", + "$\\br$. This charge on the nucleus is balanced by eight negatively charged electrons.\n", + "\n", + "\n", + "We will represent the electron density of the molecule from the eigenstates:\n", + "$$\n", + "\\rho(\\br) = \\sum_{m=1}^M f_m | \\psi_m(\\br) |^2\n", + "$$\n", + "Atom $m$ has $I(Z_m)$ basis functions (determined by $Z_m$ and the basis set in use),\n", + "with coefficients $C = [c_{mi}]_{i=1}^{I(Z_m)}$:\n", + "$$\n", + "\\psi_m(\\br) = \\sum_{i=1}^{I(Z_m)} c_{mi}~ \\phi_{Z_m,i}(\\br - \\br_m)\n", + "$$\n", + "We can think of $C$ as a jagged array of coefficients, in terms of which the overall system is described by\n", + "$$\n", + "\\psi(\\br) = \\sum_{m=1}^M \\sum_{i=1}^{I(Z_m)} c_{mi} \\phi_{Z_m,i}(\\br - \\br_m),\n", + "$$\n", + "and we can think of the task of DFT as being to determine the values $C$ which minimize\n", + "total energy. That is: the task of DFT is to determine $\\psi$, and here $\\psi$ is \n", + "specified by $C$, so the task is to determine $C$, for example by iteratively solving\n", + "Kohn-Sham equations using the self-consistent field (SCF) method.\n", "\n", - "# Integrals over Gaussian Type Orbitals (GTO)" + "\n", + "\n", + "Each function $\\phi_{Z,i}$ is defined by the basis set as a fixed linear combination \n", + "of primitive functions $p(\\br; \\nu)$, known as a \"contraction\":\n", + "$$\n", + "\\forall_{i=1}^{I(Z)}: \\phi_{Z,i}(\\br) = \\sum_{k=1}^{K(Z,i)} d_{Z,i,k} ~ p_{f(Z,i,k)}(\\br; \\nu_{Z,i,k})\n", + "$$\n", + "The values of $d,\\nu$ and the function type $f$ are read from standard tables of basis sets.\n", + "In general the lengths $K(z,i)$ of the contractions have been pre-optimised alongside the contraction coefficients $d_{z,i,k}$ and the primitive exponents $\\alpha_\\mu$ to best approximate atomic orbitals. A range of Gaussian basis sets covering the periodic table with different compute-time-vs-accuracy tradeoffs are available from the [basis set exchange](https://www.basissetexchange.org/). We use their python API in this project to provide programmatic access to these basis sets.\n", + "As these functions $\\phi$ are often taken as approximations of atomic orbitals, \n", + "this is referred to as the linear combination of atomic orbitals (LCAO) method in the literature.\n", + "\n", + "In this, we will consider only the function type \"GTO\", for Gaussian-type Orbital,\n", + "where the parameter packet $\\nu$ comprises three non-negative integers $(l,m,n)$, \n", + "called \"angular momentum quantum numbers\",\n", + "and a real value $\\alpha$, the \"exponent\".\n", + "The primitive is\n", + "$$\n", + "\\def\\normp{N}\n", + "\\def\\pgto{p_{\\text{GTO}}}\n", + "\\pgto(\\br; \\nu) = \\pgto(\\br; l,m,n,\\alpha) = \\normp(l,m,n,\\alpha) ~ x^l y^m z^n \\exp(-\\alpha\\|\\br\\|^2)\n", + "$$\n", + "where the normalizing constant $\\normp(l,m,n,\\alpha)$ is a function of $\\alpha$ and $(l, m, n)$\n", + "and is chosen so that the function integrates to 1, as derived later in this notebook.\n", + "In the following, as we are dealing only with GTO, we will simply write $p$ instead of $\\pgto$.\n", + "\n", + "Noting that a linear combination (LC) of LCs is just another LC, we will often \n", + "contract the two sets of coefficients, writing\n", + "$$\n", + "\\psi_m(\\br) = \\sum_{i=1}^{I(Z_m)} c_{mi} \\sum_{k=1}^{K(z,i)} d_{Z_m,i,k} ~ p(\\br - \\br_m; \\nu_{Z_m,i,k})\n", + "$$\n", + "as\n", + "$$\n", + "\\psi_m(\\br) = \\sum_{j=1}^{I_m} a_{mj}~ p(\\br - \\br_m; \\nu_{mj})\n", + "$$\n", + "where $I_m$ is just the total number of different $\\nu$ values in the basis set for atom $Z_m$.\n", + "For most basis sets, this will mean $I_m = \\sum_i K(Z_m,i)$.\n", + "\n", + "### This notebook\n", + "\n", + "This notebook derives several key computations involved in DFT.\n", + "In particular, DFT involves integrals of the form\n", + "$$\n", + "\\def\\op{\\mathcal{Q}}\n", + "\\langle\\psi|\\op|\\psi\\rangle = \\int \\int \\psi(\\br_1) ~\\op \\psi(\\br_2) ~g(\\br_1, \\br_2) d \\br_1 d \\br_2\n", + "$$\n", + "where $\\op\\psi$ is a transformation of function $\\psi$ by an operator $\\op$, such as gradient $\\nabla$, and the function $g$ is some function of a pair of points, e.g. $g(\\mathbf r, \\mathbf s) = \\|\\mathbf r - \\mathbf s\\|^{-1}$.\n", + "Such integrals can be written in terms of the per-atom (or \"per-center\") functions $\\psi_m$\n", + "$$\n", + "\\langle\\psi|\\op|\\psi\\rangle = \\sum_{m_1=1}^M \\sum_{m_2=1}^M \\int \\int \\psi_{m_1}(\\br_1) ~\\op \\psi_{m_2}(\\br_2) ~g(\\br_1, \\br_2) d \\br_1 d \\br_2\n", + "$$\n", + "and then in terms of the primitives $p$, where again $\\op p$ is the operator $\\op$ applied to $p(\\br;...)$.\n", + "$$\n", + "\\langle\\psi|Q|\\psi\\rangle = \\sum_{m_1=1}^M \\sum_{m_2=1}^M \n", + " \\sum_{j_1=1}^{I_{m_1}} \\sum_{j_2=1}^{I_{m_2}} \n", + " a_{m_1,j_1} a_{m_2,j_2} \n", + " \\int \\int p(\\br_1 - \\br_{m_1}; \\nu_{m_1,j_1}) ~\\op p(\\br_2 - \\br_{m_2}; \\nu_{m_2,j_2}) ~ g(\\br_1, \\br_2) d \\br_1 d \\br_2\n", + "$$\n", + "all of which depends on the ability to compute integrals of the form\n", + "$$\n", + "\\def\\bB{\\mathbf B}\n", + "\\int \\int p(\\br_1 - \\bA; \\nu_A) ~\\op p(\\br_2 - \\bB; \\nu_B) ~g(\\br_1, \\br_2) ~d \\br_1 d \\br_2\n", + "$$\n", + "We will describe how to do so in the following, but first let's briefly look at the basis functions for a simple example." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Basis set for a single-atom molecule\n", + "\n", + "Before deriving the integrals, we use the `basisset` function to build the atomic orbitals of a single oxygen atom $M=1$, $Z_1 = 8$, illustraiting typical values for the numbers of orbitals and consequent numbers of primitives.\n", + "The `Basis` object built by `basisset` consists of a list of `Orbital` objects which are defined by a set of `coefficients` and corresponding `Primitive` objects.\n", + "\n", + "A Basis set should define the following quantities:\n", + "\n", + "For each atomic number $Z$, the number of basis functions, or \"orbitals\" is $I(Z)$.\n", + "\n", + "For each basis function $i \\in \\{1..I(Z)\\}$, the number of primitives is $K(Z, i)$." ] }, { "cell_type": "code", - "execution_count": 41, + "execution_count": 16, "metadata": {}, - "outputs": [], + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "The sto-3g basis for Z=8 has I(Z)=5\n" + ] + } + ], "source": [ "import numpy as np\n", - "import IPython\n", - "from sympy import *\n", + "from pyscf_ipu.experimental.structure import Structure\n", + "from pyscf_ipu.experimental.basis import basisset\n", "\n", - "init_printing(use_unicode=True)" + "Z = 8\n", + "basisname = \"sto-3g\"\n", + "oxygen = Structure(atomic_number=np.array([Z]), position=np.zeros(3))\n", + "basis = basisset(oxygen, basisname)\n", + "print(f'The {basisname} basis for {Z=} has I(Z)={basis.num_orbitals}')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Primitive GTO\n", - "\n", - "The unnormalized Cartesian GTO (Gaussian-type Orbital) primitive is defined as follows:\n", - "\\begin{align*}\n", - "\\hat p(\\mathbf{r}; l,m,n, \\alpha) =\n", - " p(\\mathbf r; \\nu) = ~ & x^l y^m z^n e^{-\\alpha |\\mathbf{r}|^2} \\\\\n", - "& \\text{where} ~ \\mathbf{r} = (x, y, z) ~ \\text{and} ~ \\nu = (l, m, n, \\alpha)\n", - "\\end{align*}\n", - "\n", - "From which the normalized primitive is\n", - "\\begin{align*}\n", - "p(\\mathbf{r}; \\nu) = & \\frac 1 {Z_\\nu} ~ \\hat p(\\mathbf{r}; \\nu)\n", - "\\end{align*}\n", - "where the normalizing constant $Z_\\nu$ is defined so that\n", - "$$\n", - "\\int_{-\\infty}^\\infty p(\\mathbf r; \\nu)^2 \\mathrm{d}\\mathbf r = 1,\n", - "$$\n", - "that is:\n", - "\\begin{aligned}\n", - "Z_\\nu^2 = \\int \\hat p(\\mathbf{r}; \\nu)^2 d\\mathbf{r}\n", - "\\end{aligned}\n", - "\n", - "So, squaring all terms in the integrand:\n", - "\\begin{aligned}\n", - "Z_\\nu^2 & = \\iiint x^{2l} y^{2m} z^{2n} e^{-2\\alpha (x^2 + y^2 + z^2)} \\text{d}z \\text{d}y \\text{d}x\\\\\n", - " & = \\iiint x^{2l} y^{2m} z^{2n} e^{-2\\alpha x^2}e^{-2\\alpha y^2}e^{-2\\alpha z^2} \\text{d}z \\text{d}y \\text{d}x\\\\\n", - " & = \\iiint x^{2l} e^{-2\\alpha x^2} y^{2m} e^{-2\\alpha y^2} z^{2n} e^{-2\\alpha z^2} \\text{d}z \\text{d}y \\text{d}x\\\\\n", - " & = \\int x^{2l} e^{-2\\alpha x^2} \\int y^{2m} e^{-2\\alpha y^2} \\int z^{2n} e^{-2\\alpha z^2} \\text{d}z \\text{d}y \\text{d}x\\\\\n", - " & = \\int x^{2l} e^{-2\\alpha x^2} \\text{d}x \n", - " \\int y^{2m} e^{-2\\alpha y^2} \\text{d}y\n", - " \\int z^{2n} e^{-2\\alpha z^2} \\text{d}z\\\\\n", - " & = ZZ(l,\\alpha) ZZ(m,\\alpha) ZZ(n,\\alpha)\n", - "\\end{aligned}\n", - "where \n", - "$$\n", - "ZZ(k, a) = \\int_{-\\infty}^{\\infty} t^{2k} e^{- 2 a t^2} dt\n", - "$$\n" + "The number of atomic orbitals follows from the [Aufbau principle](https://en.wikipedia.org/wiki/Aufbau_principle). Applying this to a single Oxygen atom predicts the electron configuration as $1s^2 2s^2 2p^4$ which keeping in mind that the $2p$ atomic orbitals will consist of three orbitals we arrive at five atomic orbitals for a single oxygen atom.\n", + "\n", + "This rule applies when using a minimal basis set such as `\"sto-3g\"` which uses a\n", + "contraction of 3 Gaussians to approximate each atomic orbital. From this we expect a total of 15 primitives in the basis set for a single oxygen atom:" ] }, { "cell_type": "code", - "execution_count": 42, + "execution_count": 17, "metadata": {}, "outputs": [ { - "data": { - "text/latex": [ - "$\\displaystyle ZZ(k,\\alpha) = \\int\\limits_{-\\infty}^{\\infty} t^{2 k} e^{- 2 \\alpha t^{2}}\\, dt = \\left(\\frac{1}{2 \\alpha}\\right)^{k + \\frac{1}{2}} \\Gamma\\left(k + \\frac{1}{2}\\right)$" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 42, - "metadata": {}, - "output_type": "execute_result" + "name": "stdout", + "output_type": "stream", + "text": [ + "sto-3g for Z=8 has 5 basis functions/orbitals\n", + "phi_(8,0)= 0.154 * r**[0 0 0] * exp(-130.71 r**2) + 0.535 * r**[0 0 0] * exp(-23.81 r**2) + 0.445 * r**[0 0 0] * exp(-6.44 r**2)\n", + "phi_(8,1)= -0.100 * r**[0 0 0] * exp(-5.03 r**2) + 0.400 * r**[0 0 0] * exp(-1.17 r**2) + 0.700 * r**[0 0 0] * exp(-0.38 r**2)\n", + "phi_(8,2)= 0.156 * r**[1 0 0] * exp(-5.03 r**2) + 0.608 * r**[1 0 0] * exp(-1.17 r**2) + 0.392 * r**[1 0 0] * exp(-0.38 r**2)\n", + "phi_(8,3)= 0.156 * r**[0 1 0] * exp(-5.03 r**2) + 0.608 * r**[0 1 0] * exp(-1.17 r**2) + 0.392 * r**[0 1 0] * exp(-0.38 r**2)\n", + "phi_(8,4)= 0.156 * r**[0 0 1] * exp(-5.03 r**2) + 0.608 * r**[0 0 1] * exp(-1.17 r**2) + 0.392 * r**[0 0 1] * exp(-0.38 r**2)\n", + "Found 15 unique primitives for Z=8 in sto-3g\n" + ] } ], "source": [ - "a = Symbol(\"alpha\", positive=True, real=True)\n", - "k = Symbol(\"k\", integer=True, nonnegative=True)\n", - "t, x, y, z = symbols(\"t x y z\", real=True)\n", + "print(f'{basisname} for {Z=} has {len(basis.orbitals)} basis functions/orbitals')\n", + "for i,o in enumerate(basis.orbitals):\n", + " terms = (f'{coef:.3f} * r**{p.lmn} * exp(-{p.alpha:.2f} r**2)' for coef,p in zip(o.coefficients,o.primitives))\n", + " print(f'phi_({Z},{i})=', ' + '.join(terms))\n", "\n", - "ZZ = Integral(t**(2*k) * exp(-2 * a * t**2), (t, -oo, oo))\n", - "IPython.display.Math(r'ZZ(k,\\alpha) = ' + latex(ZZ) + ' = ' + latex(simplify(ZZ.doit())))" + "# Now count unique prims\n", + "prims = set()\n", + "for o in basis.orbitals:\n", + " for p in o.primitives:\n", + " prims = set.union(prims, {(*p.lmn.tolist(), float(p.alpha))})\n", + "print(f'Found {len(prims)} unique primitives for {Z=} in {basisname}')" ] }, { - "cell_type": "code", - "execution_count": 43, + "cell_type": "markdown", "metadata": {}, - "outputs": [ - { - "data": { - "text/latex": [ - "$\\displaystyle Z^2(l,m,n,\\alpha) = \\left(\\int\\limits_{-\\infty}^{\\infty} t^{2 l} e^{- 2 \\alpha t^{2}}\\, dt\\right) \\left(\\int\\limits_{-\\infty}^{\\infty} t^{2 m} e^{- 2 \\alpha t^{2}}\\, dt\\right) \\int\\limits_{-\\infty}^{\\infty} t^{2 n} e^{- 2 \\alpha t^{2}}\\, dt \\quad = \\quad \\Large \\left(\\frac{1}{2 \\alpha}\\right)^{L + \\frac{3}{2}} \\Gamma\\left(l + \\frac{1}{2}\\right) \\Gamma\\left(m + \\frac{1}{2}\\right) \\Gamma\\left(n + \\frac{1}{2}\\right)$" - ], - "text/plain": [ - "" - ] - }, - "execution_count": 43, - "metadata": {}, - "output_type": "execute_result" - } - ], "source": [ - "l,m,n,L = symbols(\"l m n L\", integer=True, nonnegative=True)\n", - "Z2 = (ZZ.subs(k, l) * ZZ.subs(k, m) * ZZ.subs(k, n))\n", - "Z2_expanded = simplify(Z2.doit()).subs(l+m+n, L)\n", - "IPython.display.Math(r'Z^2(l,m,n,\\alpha) = ' + latex(Z2) + ' \\quad = \\quad \\Large ' + latex(Z2_expanded))" + "For a basis set such as `6-31+G` we have a different number of primitives in each basis function: " ] }, { "cell_type": "code", - "execution_count": 44, + "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "Help on function _lambdifygenerated:\n", - "\n", - "_lambdifygenerated(l, m, n, alpha)\n", - " Created with lambdify. Signature:\n", - " \n", - " func(l, m, n, alpha)\n", - " \n", - " Expression:\n", - " \n", - " (1/(2*alpha))**(L + 3/2)*gamma(l + 1/2)*gamma(m + 1/2)*gamma(n + 1/2)\n", - " \n", - " Source code:\n", - " \n", - " def _lambdifygenerated(l, m, n, alpha):\n", - " return ((1/2)/alpha)**(L + 3/2)*gamma(l + 1/2)*gamma(m + 1/2)*gamma(n + 1/2)\n", - " \n", - " \n", - " Imported modules:\n", - "\n" + "The 6-31+G basis for Z=8 has I(Z)=13 basis functions/orbitals\n", + "phi_(8,0)= 0.002 * r**[0 0 0] * exp(-5484.67 r**2) + 0.014 * r**[0 0 0] * exp(-825.23 r**2) + 0.068 * r**[0 0 0] * exp(-188.05 r**2) + 0.233 * r**[0 0 0] * exp(-52.96 r**2) + 0.470 * r**[0 0 0] * exp(-16.90 r**2) + 0.359 * r**[0 0 0] * exp(-5.80 r**2)\n", + "phi_(8,1)= -0.111 * r**[0 0 0] * exp(-15.54 r**2) + -0.148 * r**[0 0 0] * exp(-3.60 r**2) + 1.131 * r**[0 0 0] * exp(-1.01 r**2)\n", + "phi_(8,2)= 1.000 * r**[0 0 0] * exp(-0.27 r**2)\n", + "phi_(8,3)= 1.000 * r**[0 0 0] * exp(-0.08 r**2)\n", + "phi_(8,4)= 0.071 * r**[1 0 0] * exp(-15.54 r**2) + 0.340 * r**[1 0 0] * exp(-3.60 r**2) + 0.727 * r**[1 0 0] * exp(-1.01 r**2)\n", + "phi_(8,5)= 0.071 * r**[0 1 0] * exp(-15.54 r**2) + 0.340 * r**[0 1 0] * exp(-3.60 r**2) + 0.727 * r**[0 1 0] * exp(-1.01 r**2)\n", + "phi_(8,6)= 0.071 * r**[0 0 1] * exp(-15.54 r**2) + 0.340 * r**[0 0 1] * exp(-3.60 r**2) + 0.727 * r**[0 0 1] * exp(-1.01 r**2)\n", + "phi_(8,7)= 1.000 * r**[1 0 0] * exp(-0.27 r**2)\n", + "phi_(8,8)= 1.000 * r**[0 1 0] * exp(-0.27 r**2)\n", + "phi_(8,9)= 1.000 * r**[0 0 1] * exp(-0.27 r**2)\n", + "phi_(8,10)= 1.000 * r**[1 0 0] * exp(-0.08 r**2)\n", + "phi_(8,11)= 1.000 * r**[0 1 0] * exp(-0.08 r**2)\n", + "phi_(8,12)= 1.000 * r**[0 0 1] * exp(-0.08 r**2)\n", + "Found 26 unique primitives for Z=8 in 6-31+G\n" ] } ], "source": [ - "f = lambdify((l,m,n, a), Z2_expanded, modules=\"scipy\")\n", - "help(f)" + "basisname = \"6-31+G\"\n", + "oxygen = Structure(atomic_number=np.array([Z]), position=np.zeros(3))\n", + "basis = basisset(oxygen, basisname)\n", + "print(f'The {basisname} basis for {Z=} has I(Z)={basis.num_orbitals} basis functions/orbitals')\n", + "\n", + "for i,o in enumerate(basis.orbitals):\n", + " terms = (f'{coef:.3f} * r**{p.lmn} * exp(-{p.alpha:.2f} r**2)' for coef,p in zip(o.coefficients,o.primitives))\n", + " print(f'phi_({Z},{i})=', ' + '.join(terms))\n", + "\n", + "# Now count unique prims\n", + "prims = set()\n", + "for o in basis.orbitals:\n", + " for p in o.primitives:\n", + " prims = set.union(prims, {(*p.lmn.tolist(), float(p.alpha))})\n", + "print(f'Found {len(prims)} unique primitives for {Z=} in {basisname}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can plot these atomic orbitals by first evaluating the basis functions on a regular grid and using py3DMol to render isosurfaces. The $1s$ and $2s$ orbitals have spherical symmetry so we plot the $2px$ orbital below." ] }, { "cell_type": "code", - "execution_count": 45, + "execution_count": 33, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAABAAAAAUCAYAAACEYr13AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAABZklEQVR4nKXTv2tVQRDF8c+LFkaLCKYQBFFQ0yqipLLyB3b+A0khaJkqjdUwAUkT1E6QNHZaBBHsBVFQTGGX1FoEIihKgiEQn8Xb4ua+vc8i0xwYZr57Drvb6/f7DlJjB9rG4VozM5/jNs5GxHZmXsYq7kXEcnO2146QmVfwCfMR8ajRf4VpnI+IrVERHuI3nrb6iziJuU4HmXkB61iOiPuVaGs4WqL9rTm4ix5eVpzBC5zGja4I17GHjx2AD0WHAZl5DBexFhHbHYDPRa/VHJzCIWx0LIuIX9gxiDEEOFH0Zxeg1A9M1gB/ih75D2C8MbsPsNlyMlSZOYbjjdl9gA18x9SI06cMrvnLECAi+niHycw81wGYLvq25gBWit7qANw0eCevRwE2MdvezMwJ3MGbiPhWBUTELp7gamZeajFmDW5oqdms/cbH+IqFxunjeICViHg/EhARO5jBannecAbPMN+e/wf9OWdT51stdwAAAABJRU5ErkJggg==", "text/latex": [ - "$\\displaystyle \\Large{z_{000}= \\frac{\\sqrt[4]{2} \\pi^{\\frac{3}{4}}}{2 \\alpha^{\\frac{3}{4}}}}$" + "$\\displaystyle \\left( \\right)$" ], "text/plain": [ - "" + "()" ] }, - "execution_count": 45, + "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "z000 = sqrt(Z2).subs(((l,0), (m,0), (n,0))).doit()\n", - "IPython.display.Math(r'\\Large{' + latex(Symbol('z000')) + '= ' + latex(z000) + '}')" + "from pyscf_ipu.experimental.mesh import uniform_mesh, molecular_orbitals\n", + "from pyscf_ipu.experimental.plot import plot_volume\n", + "\n", + "\n", + "# Evaluate the molecular orbitals on a mesh -> [num_mesh_points, num_orbitals]\n", + "mesh, axes = uniform_mesh(n=32, b=3.0)\n", + "orbitals = molecular_orbitals(basis, mesh)\n", + "orbitals = np.asarray(orbitals)\n", + "\n", + "# Mapping between orbital label -> index in basis set\n", + "# TODO: what is this mapping for 6-31+G?\n", + "#orbital_idx = dict(zip([\"1s\", \"2s\", \"2px\", \"2py\", \"2pz\"], range(5)))\n", + "plot_volume(oxygen, orbitals[:, 5], axes)\n", + "()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Integrals over a single primitive\n", + "\n", + "The GTO primitive\n", + "$$ \n", + "p(\\br; l,m,n,\\alpha) = \\normp(l,m,n,\\alpha) ~ x^l y^m z^n e^{-\\alpha\\|\\br\\|^2}\n", + "$$\n", + "contains a normalization $\\normp$, which should be chosen so that $p^2$ integrates to 1:\n", + "$$\n", + "\\int p(\\br; l,m,n,\\alpha)^2 d\\br = 1\n", + "$$\n", + "That is, \n", + "$$\n", + "\\normp(l,m,n,\\alpha)^2 \\int x^{2l} y^{2m} z^{2n} e^{-2\\alpha\\|\\br\\|^2} d\\br = 1\n", + "$$\n", + "We note that $x^{2l} y^{2m} z^{2n} e^{-2\\alpha\\|\\br\\|^2}$ can be written\n", + "$$\n", + "\\def\\ldp{h} % Read this \"1-D primitive\"\n", + "\\def\\ldpi{H} % read this \"1-D primitive's integral\"\n", + "x^{2l} e^{-2\\alpha x^2} \\cdot\n", + "y^{2m} e^{-2\\alpha y^2} \\cdot\n", + "z^{2n} e^{-2\\alpha z^2}\n", + "= \n", + "\\ldp(x;2l,2\\alpha) ~ \\ldp(y;2m,2\\alpha) ~ \\ldp(z;2n,2\\alpha)\n", + "$$\n", + "for\n", + "$$\n", + "\\tag{defh}\n", + "\\ldp(t; 2k, \\eta) = t^{2k} e^{-\\eta t^2}\n", + "$$\n", + "So the integral\n", + "$$\n", + "\\int x^{2l} y^{2m} z^{2n} \\exp(-2\\alpha\\|\\br\\|^2) d\\br = \\int \\ldp(x;2l,2\\alpha) ~ \\ldp(y;2m,2\\alpha) ~ \\ldp(z;2n,2\\alpha) dx dy dz\n", + "$$\n", + "is the product of three independent integrals\n", + "$$\n", + "\\int \\ldp(x;2l,2\\alpha) dx \\cdot \n", + "\\int \\ldp(y;2m,2\\alpha) dy \\cdot\n", + "\\int \\ldp(z;2n,2\\alpha) dz\n", + "$$\n", + "Where each is of the form\n", + "$$\n", + "\\tag{defH}\n", + "H(2k, \\eta) = \\int t^{2k} e^{-\\eta t^2} dt\n", + "$$\n", + "These one dimensional integrals have a known analytic solution [see equation 42 of gaussian integral](https://mathworld.wolfram.com/GaussianIntegral.html), but we will lean on [SymPy](https://docs.sympy.org/latest/index.html) to help us derive a formula for normalising our primitive Gaussian functions." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, let's define the integral:" ] }, { "cell_type": "code", - "execution_count": 55, + "execution_count": 20, "metadata": {}, "outputs": [ { "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAO4AAAAkCAYAAACOnQbwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAAMKElEQVR4nO2deZBcRR3HP5FEIocG0KACCpQIpAK7mwiIHHIGUUQMoiIEvCgUBWMQJYXw5SsCAZSzVERUImAhEAHD4YGI8QCBkOUGsdQACoZLMRCuEP/49Ute3r6Z2dnM7OzG962a6pn3erp/r6d//Tu7Z8SSJUuoUKHC8MLIThNQYXBg+3hgf+B24LPAO4GTgUXANEk3d466Cs3iVZ0moEL7YXs34I3ARIJxfwqcAnwsvU61PaJzFFZoFhXjDnPYHmF7lQbVeoCZkhZKOg1YFThb0n2S/g7cD7y+zaRWaCEqVXmYwvZoQMBngFG2vw1Ml7S4pPoDwHuAm2y/m2DcabavBZ4BNgOeHBzKVy7Y3gC4EBgLvAycIOmydvdbSdzhi+8BRwNjgNWBo4Cv1Kg7G1jX9kPAN4APATOAW4C7gTMkvdJugldSvAxMlTQOmAScaXv1dnc6ovIqDz/YHgs8BowAjkjvLwUekbRBJ2kbCGzPACZK2r3TtKwobN8B7CXp4Xb2U6nKwxPdBNMCXAz8B7gWeNz2apKe6xRhedieDkwGNgVeAG4m1Pm7C1W7gd5BJW4AsH0asKWkPWrcnwis0m6mhYpxhysyR9Lzkp5K79/XKWLqYCfg28CtxELzNeB62+NydEMw7o8GnbrmsTUwp+yG7bWJZzhkMAipGHd4Ys1UPtNRKhqgKJlsTyG0g+0IuxvbbwTWJUncZB9+H3g7MDl5vdsO273AVcDaRLz7FeAcSSfYfjWwEBgF7Gj7q8B9ya7F9qrAlcAMSX8cDHpbzri2ZwJ7AhtJenYF2/oicDpwgKQft4K+AdAwEbgNOETS+Z2goQRrpHJhR6loHmsSDtGnc9e6iSSQB2xvSsSY5wHbSVo0GETZHgVsnug7HjgLmBK3/D1gAbAtMQ+2AR4iVH9S/PsC4AZJF7aAlg1S+1dImlyrXh/Gtf1LYHdgZ0k31mj8XOBQ4FBJ5+Wub0U88JcyprW9DvBBQpXbAlgPeBG4C/gh8MM6Hs2Jqbyt5pO2GZLm2r4SOMH2JZKGArNkEve/HaWCphfXswjJelPuWjcxF/YBziPCKWe2ms4GGAe8GjhS0pUAts8HjgPGSHrM9puI8b5VUt6jux3wEeBO2/uka1Mk3TVAWrI5f3v+YnGcyyTuRGAxYZfUwjapLDLUiYT69p3ctf3S50eB3xCrybqE0+J8YE/b+xUGI8OE1N6DdWgZDJwM/Inw4J7UYVpgaEncbKLNrVfJ9unA9sD2hVhzN7AJ8ANgb0m/bQVRtr8OHNOgWiacuoCnSOp7whtS+Vgqe4A7ivNU0u9pbVh1QiqL47ncOC/Xoe2NCR3/nlpqru3VgPGEqnBX7vrbgd2ASwsqzp+BvYH1JR0gabqkTxJB/4eBfQkmLvazOuGNnFeDqQcNkm4hsosOtT0UYt9DRuISE2oh8TuXwvYZhN24i6S/Fm53E+rxKGLutQpnEupvvdctORrmliwo8yX9O/d5Xgvpq4VSiUthnIsS9x2p/FODhkcSDPVS7vonCc/hT/KVJd1Q1khSP84lpPROwKxClW5iYVlu5bG9FmFT7E2oXkcV6GgI2/cAGxHe2S8Q6v3GhO11EXB0SQbSJYT9szvwi2b6awM6LnFT7DWf8PGK7ez9QZm9Z/ssQpXcWdL9hTZWI6TtFEIbu9D2jpKKkzarvy/waWAr4LXAfMLcOqX4e0l6Aniin4/TRV/tsYflQ1RdwHX9bK8ubI8EPgd8inDC/Yvwvp9KSNx/SvpXqls6zgNh3HemsviguxEqdjO7TDKGe7nkXh+VwfY2xMKwFrCfpMub6Ctr4zWEJJ9PTJbRwA3A74mE+y8BjxCLQh5/SOVQYNyhIHFvB2YCBwN/BH6Vu3cjgO1vEUy5D/B08iADLEy+gi2BJcDdkm61vRkw2/bWkv6RNZZysS8CPgr8BbiM0Pj2JBb+TRMdA0UXkYmWR0/hmUYCm9l+M/BcThI3heShvpqYR73AOcA6hFB4G7EZJK+yl45zLcb9iO2da/Tdx75Nam034SLvlyc5rToHpY8/L6mynE5vexqRpnc/sLukgdq9XcAqhIQ9Dzg1U8Vt30gkNOxCX8bNbP4d6zVueyqRhthf9GYOkSbQcYkr6VLbY4gJNTPvpMzhsFT+unDdxETtBh7MmVbHEUz4M9s75BJJziKYdgZwrKSXAWwfRSwSB9k+RdK9zT6H7fUJxunNXRtBzJNTc1WPIXZUTQO+S2yNHAi+RTDtccDXc3PvAiCz75dqHLXGeSnjJmIzKdef1LO8xF2PYIZHm3iAGYStfK2kMgk2gZAoT9i+ilCNZwKfXcEwQU8qr5B0SuFeNmB97FhJ/7H9PPCWBu1PBd7aBD0ziRhgMxgKEheWzZdS1VZS3a2Cks4Fzs19XgJ8OF8naVmHAVdJml74/ksp/PguQqA0zbiSHmFZFlqGTYjFcV6u3sXEoj5g2N6aUPWvlnRCgY45tu8jbO+iY6rPOOcl7ibA64CbJW1bo+O1CO/bIuCe3K11Uvl0ny+Vt3MEcCQhPaeU3B+dHmBBInYsrYujZoNQJiHGprLoQMnwFOERrwlJGw6MrKbQcYmbMIEwdwYa+ugPDicY67l0GEAR41PZSqdhD/BkG1IXD0/liTXuZzu0igthn3HOM26mJpeunrkGINS7vDMgk4Cj63wXANufJ1Sfe4FdC6lvGboSbasCbwYuamHyQw9hi5eFHbpSWct7+BqWPWsn0XGJm0ydLYB7Jb3Qxq4mpXL/BvXmt7DPomOqVZhEMGctH9LGwIKCfV86zgNl3KJjakEq16EOkv13BrGVbFdJC2pUzfo5ktiCdqDtuSsamE8ZMuMJW7yMAWvGJFMYaAzwtwZ9TOX/wMYlkhZGU3++rBCS5vUGYI6kd7ernyIkHd3qNtOzjKVGeNP2BEJIFT3XpePcLONm9mGRcR8FHiccC7UI/wph1/YSzqV6rvp8xtTlhHT8pu2HJP20zvcaYRwhxWslC0wkJGqZrbQpobL1NuhjKu23cTPG7aSN253KdsY2M9tzZTidY3F6ja1xP0sWKfJfdyqXG+eRsFSa9BCpiPdQG6USV9IS23OAfW2/TdJf8vdtH0vsDJkLTKqhHhf7WQTcL2mx7b2IMNNFtndRycFmySt3MPAJSRc0oL/P4pRiipsTKW1lp0hkYbDf1CO83TZuCmdlR9V0UuJm2lXbNjpIWmT7TmBL25PLFm3b2wM31fjNhgySI+1BIqT0fklLQz5JqGVJSEWhUjrOmcTdjORFk/RiWce21yAcWAsJp1IRs4gsqD2IWFv2vYMJpl0M/A44Iheoz/D3jNlSnGt8omVxeujHbL+XiKXOtr1tcXFgmXOiLCacIdMYyiRuFiaqJY0npWe4qk77g4E1c+87KXGzcTrR9njgWSLjrtXHthwFXAPMsn09cCfxW69HaEijJDXy9A8VnExoWLNsX0KkU+5E2LAPAxvQV6iUjnM22fujJncRAzavxqaAWYSte1Dh+kapXIVQI1Xy+niu/hZE+ttyqkGK0U0mPN/X2S6qT1sQE/maOs/QQ2zX6i25VytHFNuvI5IIrm6Dp7FZrJF738k47hzCS/psKsUyta6V/fySCPfMIn7jI4hIxObA9cAnWt1nuyDpRwQPPEI42w4mGHY7IhHlSUnzC98pHeeWHl2TTjw4CZggaTDyOrN+xxDeum9K+vIA2/gBMQm6JN1ZuHc4cDawQ0oq7xhsd7Fs4dlYUl1nWYWVE61OmD+D2P3ztRa32wg7EHGu01egjYnA8xQcU8mmnA7M6jTTJgwVVblCB9HSjfSSnneccrCz7dX7m/7Ygn5n048Yci0kV/04YodI0UbekEjWuGCg7a8IUrbN8YQf4DDCDoKw5Rs5+SqspGj5CRhJJy89l2cII0v46GPfSrqPYJxO4QkimX5ROo3hwHS9t4avocL/AarjWYcBbM8G9ipc/pCk4lbIFe3neKr/FxoWGAqbwis0xhQiwf0Z4A7gA21g2ur/hYYRqlMehwHS3s8DG9XLo8mjWyD3/0LAaen8pLOTqYDt7P+FHm+GjgrtQcW4Ky/OJDaf18NDuffV/wsNI1Q2bgVg6X7s7wDvJY5S2Y8Is51EeLCnSbqicxRWyKNi3AoVhiEq51SFCsMQ/wMZ4YokX4+DYgAAAABJRU5ErkJggg==", "text/latex": [ - "\\begin{align*}~\\\\z_{000}= z_{000}\\\\z_{100}= \\frac{z_{000}}{2 \\sqrt{\\alpha}}\\\\z_{010}= \\frac{z_{000}}{2 \\sqrt{\\alpha}}\\\\z_{001}= \\frac{z_{000}}{2 \\sqrt{\\alpha}}\\\\z_{110}= \\frac{z_{000}}{4 \\alpha}\\\\z_{101}= \\frac{z_{000}}{4 \\alpha}\\\\z_{011}= \\frac{z_{000}}{4 \\alpha}\\\\z_{111}= \\frac{z_{000}}{8 \\alpha^{\\frac{3}{2}}}\\\\z_{200}= \\frac{\\sqrt{3} z_{000}}{4 \\alpha}\\\\z_{220}= \\frac{3 z_{000}}{16 \\alpha^{2}}\\end{align*}" + "$\\displaystyle H{\\left(2 k,\\eta \\right)} = \\int\\limits_{-\\infty}^{\\infty} t^{2 k} e^{- \\eta t^{2}}\\, dt$" ], "text/plain": [ - "" + " ∞ \n", + " ⌠ \n", + " ⎮ 2 \n", + " ⎮ 2⋅k -η⋅t \n", + "H(2⋅k, η) = ⎮ t ⋅ℯ dt\n", + " ⌡ \n", + " -∞ " ] }, - "execution_count": 55, "metadata": {}, - "output_type": "execute_result" + "output_type": "display_data" } ], "source": [ - "out= r'\\begin{align*}~'\n", - "for vl,vm,vn in ((0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,1,0),(1,0,1),(0,1,1),(1,1,1),(2,0,0),(2,2,0)):\n", - " zllm = sqrt(Z2).subs(((l,vl), (m,vm), (n,vn))).doit() / z000 * Symbol('z_000')\n", - " out += r'\\\\' + latex(Symbol(f'z{vl}{vm}{vn}')) + '= ' + latex(zllm)\n", - "out += r'\\end{align*}'\n", - "IPython.display.Latex(out)" + "import numpy as np\n", + "import IPython\n", + "from sympy import *\n", + "\n", + "init_printing(use_unicode=True)\n", + "\n", + "eta = Symbol(\"eta\", positive=True, real=True)\n", + "k = Symbol(\"k\", integer=True, nonnegative=True)\n", + "t = Symbol(\"t\", real=True)\n", + "H = Function(\"H\")\n", + "\n", + "H_na = Integral(t**(2*k) * exp(-eta * t**2), (t, -oo, oo))\n", + "display(Eq(H(2*k, eta), H_na))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# From factorial" + "Confirm that the equation displayed is the same as that in (defH).\n", + "\n", + "And then solve it:" ] }, { "cell_type": "code", - "execution_count": 47, + "execution_count": 21, "metadata": {}, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "1.0 1.0\n", - "1.0 1.0\n", - "3.0000000000000004 3.0\n", - "15.000000000000004 15.0\n", - "135135.00000000003 135135.0\n", - "316234143225.00006 316234143225.0\n" - ] + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ8AAAAXCAYAAAAV3cp8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAAJdklEQVR4nO2cf7BVVRXHP4A/8Ecp4miTkwUDAk3Ek4cJKRQqmmXyQ8mmRLCZxukXYWJFZV+XJqGTAs1kDEPKm9BsRhAj0JwUUhxRJEgQKZ0i+kERQSQWkUB/rH3gvnPPOffcc+9979Hc78yd/e7Ze6+93jnfvfbaa69zux06dIgmmjCzx4GNkm5u4BijgOlAK/B24HpJCxPazQP2SZpmZjOACcAA4D/AGmCGpE2N0rOJfOje2Qo00WXQAqxv8BgnA5uALwL/TmpgZt2AK4Gl4dIHgXuB9wMXAW8CPzez0xqsaxMVcEy9BZpZG3A50EfSGzXKuhG4B/ikpAfroV8BHVqBF4FPS1rQGTo0Gmb2NuBMYEP4fhLwA+AcYIKkrfUYR9IKYEUYY2FKs/OA44HVoc9lMV0nAXuAC4BlteiTxdWuwL2ugCz+lxkPM3sCGAOMlrQqReA84AbgBknzS66fB0wCpkcPw8x6A+OBjwCDgbOA/cBG4H7gfkkHU3RvDeWLuf7TBkDSOjNbCtxuZg9J2ttZujQQLbgn8GszGwAswb2QCyQleggNxDhguaQ3U+rfgnvMu0svmlme/XcvSf8I7cu4GkOncy+CmV0NfAB/TkPwe/CApGtz9D0T+DPuvd1KlXMxi/9JnkcrcABYm6HT+aGM39g7gH8C3y+5NjF83w6sBLbhq9wEYAFwuZlNlJT08IcGea9m6NIR+DbwPDAVmNnJujQCLTiBxgHzgdslzekkXcYCt2TUz8U9pOdS6i2j776Sv5O4Woquwj2Ab+BGYy/wR2BgFX3H4sZ2CcXnYiL/2xkPM+sLnAa8lLblMLMTgffgwauNJdfPAS4BFsRWq9/ge9jlpVbNzL4GvABcFZRfHBvnJDxI9kyKYekwSHrBzLYAN5jZrAxPqcvAzL4FfL1Cs8i7bAH6A/cBV0r6RZ3kVgUz6wf0BX6WUn8PcCFwoaQDSW0k3ZpjnDSuRvV1556ZTcFX9yL35kbcaLyGeyArq+g7Hvg78DQQxZOqmotp/I97HsNC+XyGMq2h33pJ/y25/qmg3I9jAz+VJETSX8L25w48KLY41qQFt5jrSi+aWS9gIX4T5gI3x/SoCDN7GegDnI4H7ybhpN0NLAK+mkDOh3C3bwwp5G4kCug8J1zPwrZQtuAr0yfwxSML1citFuOAJ5MWLjObDXwcn3y/LSg/QiJXS9BCg7hXBJIOGwuzLMeqPczsrXiQ+cHAjaJzERL4X8R4DA9lfMtyCb7dWZPRN47oxiftb4eG8vADNLPz8QfeC5go6eEqxopknICvKr/HLXhP/KauxifPdNzKz411fTaUHW48iugsaSewM4fsE3GvY1KQ/UMzGyXpl0nt88otiLFAW4KOc4FrcMOxpQ7jVOJqQ7jXCbgCOA54JEfbrLkICfxPMx7XmNnoFCFl8Y7g5rUAr+Q9YTGzY4DrwtfHE5pEAat1of2XgFnAFmCMpKJ70SFAD3zVng/cFbmmZrYKeAC31nHjEcWARmUJN7NpwKlV6LNB0tIG6ZwH7wUOAZskrTWzgcAyM3ufpD8VkJcKMzsZ6Be+dgfONrMWYBcesB0OXB3r8z3csI0DdoeTIYC9RYLXObnaKO51NMYDbwBPZDXKMRchgf+HjUc4X48s7pgcipV6Hmfh5N6eo1+EWXjsZIWkpJV8KPA6sNPMHsVdxTbgMzWeAJwbykck3Rmri1bbsvwXSXvMbB9wdgX504B3VqFPG0dyGtJQSOecaAFeLbmn38S9nJ+Y2UhJ/yooNwnDaL9ft/Bpw/fkayX9Ndbns6F8MnbdcDe6WuThaqO412Ews57Ah4DHJO2r0LzSXEzkf6nn0R84BVgjaUSKQr04skq8XFLVO5S7yzoly5kK3IRb8kkJ9T2BQcAOfHKcQf3yLCIDOT+h7oxQpu2pd+HR6VRIelcxtTJRi86ZkDQPmFfy/RDwsSKycoy1Co81lCFM0qUJfRLb14BMrtaDe2a2lfQFZGVC3KJN0pS88nNiDJ6Ul7llqTQXY2jH/1LjEW1ZEve6ARGJN8SCc5E17llhcMzs87h7vRm4WNKuhGZDgm7H42nMi+qYoHUuvt9NOlEYEsq0TMsTSMmMbDBq0flowbPAjzpgnEpcrQf35lC+dW3hSExna6xuQ5Xy82ACnsOxPK1BzrlYinb8L2o84sHSHaHsTQZCPGA2nqJ8saQdKU2jcW7C98DXmtm6WnMPzOxY3D17JcX9bLfXjfXtjhPidxXGmEYdYx616Hw0QdJdHTRUJa7WzL2ktuGodiywsMgxdjUwsx7AR4GnJO1JaTONfHMxal/G/2qNR7T3jhuP7cDf8L1y2uBfwfdWG/CgU1bEvjS772F8xb3bzLZJWpLRrxLeja8oaROtFbesmxPqBuAu94YKY0yjvjGPWnRuohyVuNoo7nUkRuHGMXHLUuVcjFDG/2OCsO64YdhP+1hGHImeh6RDZvY0cJWZ9ZP0WkzZW4Db8AlwaQ73aCg+IbZIOmBmV+DHaovM7CJJZUds4V2JyaS8qRnTv8xAhiPLQXjQLikBKTqizkzQaUDMoxadm4ihElcpwL0uiPHAQeDReEWBuRihjP+R5zEQD66sl7Q/qWc4ZuuPp8gmnbUvxjPULsMz4aJ+k4OyB4BngKkJAaOt0YQ3s+NwN319NCFCEsuH8X3xMjMbkfDQo9OGtHNqOOI5Ja3i0XFo2gp/afgfyh5Ig1GLzk0kI42rRbnXUJjZOPyoGiA6qh5R8nLhTknTQ9tuoe1z8ZOraudiDGX8jyZcni3LkNB+fUp69mJ8P3ld7HqfUPbAXXolfKaUtB8MHEssAChpMx4EOgV4zMxOj40zGD9eSw0Q4RPxIMlbj7LEoAhmdgr+QH4q6Q8Z8huBQjo3kYk0rhblXqPRgnvVk3GDB57zE10rzY0ZBrwDzxiOo9q5CKTzv1s9fwzI/IdbZgJDJXVY9N/MTsXz9++W9OWCMu4DrgeGSHopVvcF4LvASEmra1S3bsjSuYlsdBZXGw0zmwnMAPpKygzuVyEzkf/1/jGg2fh7DbfVWW4ljMTTa++pQUYr/tZlu8BjSA2fASzuSoYjIFHnJnKhs7jaaIwHflVHw5HK/7p6HmGwUcBo4Dt5U9U7GyEx6HVgnaThsbpB+HsVC1WnH8WpB7J0biIfjkaudjSy+F9343E0Irz0tAa4V9LnOlufPDgadW7i/wtN49FEE00UQvMHkJtooolC+B8c22f8qZDUrAAAAABJRU5ErkJggg==", + "text/latex": [ + "$\\displaystyle H{\\left(2 k,\\eta \\right)} = \\eta^{- k - \\frac{1}{2}} \\Gamma\\left(k + \\frac{1}{2}\\right)$" + ], + "text/plain": [ + " -k - 1/2 \n", + "H(2⋅k, η) = η ⋅Γ(k + 1/2)" + ] + }, + "metadata": {}, + "output_type": "display_data" } ], "source": [ - "import scipy\n", - "p = Symbol('p', integer=True)\n", - "\n", - "# n!! = special.gamma(n/2+1)*2**((m+1)/2)/sqrt(pi) n odd\n", - "# = 2**(n/2) * (n/2)! n even\n", - "\n", - "def fac2(n):\n", - " assert n % 2 == 1\n", - " return gamma(n/2+1)*2**((n+1)/2)/sqrt(pi) # n odd\n", - "\n", - "\n", - "for n in (-1,1,3,5,13,23):\n", - " v1 = scipy.special.factorial2(n)\n", - " v2 = float(fac2(n).evalf())\n", - " print(v1, v2)\n", - " assert np.isclose(v1, v2)\n" + "display(Eq(H(2*k, eta), H_na.doit().simplify()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And, recalling that we are looking for the normalization $\\normp(l,m,n,\\alpha)$, which is given by\n", + "$$\n", + "\\normp(l,m,n,\\alpha) = \\left(\\int x^{2l} y^{2m} z^{2n} \\exp(-2\\alpha\\|\\br\\|^2) d\\br\\right)^{-\\frac12} = \\biggl(H(2l,2\\alpha)~H(2m,2\\alpha)~H(2n,2\\alpha)\\biggr)^{-\\frac12}\n", + "$$" ] }, { "cell_type": "code", - "execution_count": 48, + "execution_count": 22, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcAAAAAVCAYAAADRoT5bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABJ0AAASdAHeZh94AAALQ0lEQVR4nO2df7BVVRXHP09UIDBTHLWUHxIpKvYuvGxsGlEDLNLSZzo4YyQ2lUORv8IxmWyx1MxGDLQpqcERnJE0gn6pUDAqaigGSagpmQiYP1PMQFMR6Y+1D5537jn3nnvPvvfcK+c7wxzePvvstc/6sfdea699bseOHTvYVaGqS4BHROTivPvSLlDV0cBUoAv4CHCOiMyt8sx0QCLFL4rIgY3oY4ECBQqkwW55dyBnlICHq1VS1bluEG8qmkW3Rjr9gUeB84H/1UBmHfDh0L+jauljgQIFCvjG7nGFqjoPGA8cIiKvN7dLfqCqXcAq4OsiMifm/oHAAcAajzQPBp4BzgTOBY4BngQmYBPHj4FPYBPI6SKyqQG0TwUmA8cCLwBfE5G7fdERkTuBOx3NuTU8+o6IvOCrH2FU0ldVvRDj+1kiMr8R9FsF1XQ+Q7ttPx7kiUIu+aIS/8smQFU9GpgITA0zVVVPB47DvKZOYC/gFhH5cooOHAA8B/wMmA50AydhXsBBwNvAI8BNwE0i8m6tLxmFiKxW1d8CV6jqrSKyNVKlhHkw67LSCqHTXSdjIb/NwAJgLvAGMA3YAvwGuAi4oAG0LwKuAKYA12KD/0iPdOrFUFV9DngLWAlME5H10UqqmiYmv4+I/MfVj9XXELrc9a919boO1GsrWe0khc43gr9tg7zGsEIu+c4flfgf5wH+APgvcEOk/Huu41uBfwHDq3U8hFOwcOsi4AzX9vPA3cAmzBM7DZgDjFfVM0TEx+bkD7HB9jzgqsi9Erb/t90DnXCbrwETRORFAFVdinmEw0XkFVe2HAsD+kQJk9uEwNNS1V9jPMgbK4FJwBPA/pgurVDVIwOexEArtPdm6P9J+hpgFKaz/6ilwxlRr634sJNKOh+GL/62E/Icw3Z1ueQ9f8Tyv8cEqKqHAmOBOSIS3d+50HX8n9hMXktYrRt4BbgX6AC+CNwRnqlVdRrwEPAl9zILa2g/FiLykKo+AZyrqldHVgYlEsKfri/TQkW9gR2qOjVUNl5E7os82gncHkx+DoOARZGBfhDwF490A9p3RMKMwzB5+aRTM0RkcYT+g8B64GzMQ417Znq1dqvoK6raDzgMWFFLVEFVJ2GryRNE5J60z4VQr61ktpMqOh+uN71aZ6rxt9nIUS6QUTbvB7lk5H+u80cS/6Me4FcdgduiPQjvI6lWWqT0hKp+EPgMMN95W3fF1RORF1R1NrayOT76Aq6tPYBvYd7EYRhTFgCXuHd5GlgmImeFHrsVc5vHAX8MlZeA6xK6PRv4VejvHwHPAteHyp6Nea4UqQMWfvx+pKwT+IVHugHtn8TQXuOZTmaIyFZVfQz4WMamEvXVoYStHJsW/oT6bMWnnZCs87WiIn9V9QTXz2uB+ZiejwZ6ufLJrr9HYB7AWGyxdR/wTZ974GnQAmNYs+QyDvgTcDVwC7bYHQP0A9YCF4rIygz0a0YL8B5i+B+dAMcC24EHU/ewOk4G9sT2vaphm7u+E72hqvsCS4CjgduxFzgZ20d7FngX2JfydPs/u+vOl1bVD2CDb2wGqIhsxvbvAtpbgM0iUuZNher0Az4ablNVBwADI2UDgQFxtOuhm0TbYSQWNvBCxxdUtQ8WAsmanFNNX0e56+qMdJoBL3biUKbzdSItfw8F7gcWAzcCn8OSsXqr6g3AL7GBax42OJ0E3Oz+3w7wJZtmySXY8z8UizQtxXh/JPB54PeqOkxEtmToQ7PQULvYOQG6QbQEPO55U7UbeB1bkSRCVXcHvuL+XBJT5TZs8jtfRK53z1yDudXjgSOAuTGDeBBqHB0q+7i7rk35DmkQtLkmVFbCkj4ei5RtJSY06ZO2m3wPxmOWq2u3PxZaBfOuBqlqCZtAN7k6U4ApIjLc/T0D+AMWr98fuAxbjc7L0I80+tr0BJgM8GUnEK/zNSElf4MJ8JPAMSKy1j17OSbrE7HBeJyIPODu7Ynp/mhV7SMib5Y323LwJZtmy+VY4NMislP/VXUhFiIsYZ54q6OhdhE+B3gQFrp4vq5uxneqD7YaXJxC0a8GRgB3ikiP1ZGqjsVWPfcRCvOJyMvABsxF3ge4PNqoiLyGbR4PChWXgCdF5I3a3qgiOl2bYaUcCTwqIu9E6v3NR6ZrDO1wdtlIbEX0d490wI5xPOz+9cU27R+mJ+/3w0LUAQ7GvIB1mEf6FjZgbszQjzT6OgrL9H08A52Gw5edBEjQ+VqRlr8Ak4LJz9HfgtllL+DiYPJz997G9KADWwS1NHzKpolyCTzAc8KTn0NgC30y9KEpaIZdhEOgA9z11bp7XI5x2Pm3iu6rqp4HfAfLEpwYUyUomxWT3RMw5uci8kwCic1YphAAIjIb2wdLBRGZlKJOWZsiMgOYESm7ErjSF90KtJdhoQNvdFzde7DBq1Kd6VisPfj7zLTt14CK+uqM53BgVaVMX1XdAAxOuH13zH7FvFr4lRK+7CSMHjpfB6rxtx8WYlsvInEr7sGuD3H7VIOBLRUygFtFLuBfNo2WS39se2cj7rxuBEPd9alKRFqE/w23i/AEGGQT+VwZnIad0bgjqYILl12HeSpj3P5UFMdh3kySa/sGlVOL+1LbV0sKtD6q6Wsnpt/Vwp+zgA9FykpY6vU8zJMJY0267tUEX3YSRladT8Pf3YBl0RuqOgSLyCwSkW2Re0EI/f4q9GeRv1zAv2waLZcStkBdmnCUbBR2VOvpKnRmkT//G24X4QnwJXcdgAeoai/gC8BdzvWMq3MBMBP7MsoYEXkppk5fzGV9KhqyVNWhWDLFisjRg3Cd3TBBVhN4gfZCNX1NlQAjIrOiZS7d+xRsT/me+rqXDr7sJFLfh85n4W9XhXsjsQG64sIkb7k4el5lk7dcVHUvzGtfnjA57kTe/G+WXYT3AJ8H/k3PvZssGI0JKdZ9VdVLsM6vwc6VJHW+L2YwcXtmM7G06qSsH7D36aBxK8QC+aCavrZLAowvOwnDh85X428w0K6KuddV4V6wP9XqcgH/smmGXAL+JvG+g/bIim6KXeycAN2K4F5gP1UdVvZo7ejGJq3fRW+o6mXYpuVqbOZ+uUI7r2JZk8NUNch2RFUnYwciodxVD+MYd/X2PcwC+SOFvo6iPAO3FeHLTsLIrPMp+fs2tvqOopIHGEyc7TAB+pZNM+USl+HeLotCaJJdRM8BLsRO0n+WSJq+qp6Kne0BCH7G5lP63geRXxaRqa5uh6v7QDQ0qapnYxmD27GszvNiNlQ3iPuJHRHZ4WhMAZap6gJHvxtjzt7A8e4Q5I0iEv3CyomOVhkjC7Q9YvXVpdqPwLJttyU82zCktRWfdhKBL51P4m9v7NjRWpfVGUUXsDEhySW3zNy8xrAQWkEukIMH2AK8hxj+x02AL2HnKX4auVfCPl0VxlDeyyraiP1OHFiq/EBsIzWKQ9y1F8kfg16OfUA6wMXYan4C8A3MK5wJfBdz62/Gfn2hx+l/Vd0bY+TtFTJEC7QvkvT1KGAP8lvplkhnK77txLfOJ/F3BMbfuH2mwVjoannMvd5YZu7qSpm5DUSJnMawJspld+LDn2ATYLO/ixugRH7zRyL/O6I/iKuql2IZlaNEJPZLKdWgqlcBlwJDRSS35BNV/Tb2ea9jRaRa1lmBNoQPfc0LjbAT3zrfzvzNAt+yKeSSHs20i7gfxJ2JfcWh7FB5DejGwk95Tn59MSYuLCa/9zV86Gte8GonDdL5duZvFniTTSGXmtE0uyibAN2J+4nAKnfYtWaIyOEiUqrnWY8Ygn1wemqVegXaGD70NS80wE6G4Fnn25m/WeBZNkMo5JIazbSLshBogQIFChQosCvg/4NessAgH0GVAAAAAElFTkSuQmCC", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAAmCAYAAADgMR6zAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAASBUlEQVR4nO2debRdVX3HP4EwJxSBMpVQRgkUeffxCqgQIC0iRa2U0gIqk1JkKBpiEQHhl68ESYhggOKyuloSEOpitExFRIbIZJgeM2kYInMUEDAEAoH0j9++vPNO7rv3nDsP+7PWW+e+ffbe53fP/d5z9/Dbvz1q2bJlRHoDSaPNbGmd65wCHAQ8ABwNfBI4E3gHmGxm99TzepHWUS/9NEKHod4pRC1GAo3SWYVrTiFqsCdZodUGRBqPpFGSDFi3zvXuCWwADOAPj6uA6cCXwt9ZkkbV85qRlnK6pJWqLdwoHYa6oxYjaWrSa16iBnub2JjqDaYDd5vZK1kyhx+9FTNk7Qdmm9kiM5sBrAKcZ2ZPmNkC4Eka8MMZaRmXAT+uoXwuHULUYqQmatUrEDUYycboVhsQaSySvgSMMbObMuRdFTDgKGAlST8CTjKzD0YoMg/YG7hb0u74w2OypBuAt4DxwGt1eBs9j6RxwMXAesBS4HQzu7yZNpjZg5IWSDrOzM7PUzaPDkP+qMUOptP1ClGDnUarNRdHproYSesB3we+m7HIT4HvAGsBawAnACeWyX8tsL6k54AfAPsD04C5wKPAD83sw6qMj6RZCkwys22BvYCZktZogR1nAcdI+njWAlXoEKIWO52O1WuCqMHOoqWaGxUd0LsXSRcCL5nZKRnyrge8AowCvhFeXwa8YGbjGmpoA5A0DRgws8+02pZGIOkh4PNm9nwLrn0k8E9Z720eHYb8XaXFSGfpNZSJGuxwmq25jpjmk7Qb8G+4Y99GwOFmNqulRrU5kjbHnR63yFikgD84AC4B3gRuAP4gaXUzW1x3I6tA0knAfsDWwBLgHnzo/dFU1gIw2Ka2ZalrBrC9mX22xLkBYMVW/DAFZgNTJf2tmf26XMYqdAjdp8W2JOq1LAWiButON2uuIxpTwBh8mPSi8BepzLeAW83shYz5i46R75rZ6+H15+pvVs3sAfwIuBd/2H0PuFnStgm7wR+GzdZKVtuysBMwJ50oaW38ff1LbaZWj5ktkfRzfNqu0o9TXh1C92mxXdmDqNeRiBpsDHvQpZrriMaUmd2A9wqQNKu11rQ/klbGRwNOzlFsbDi+VX+L6ke6FyLpYLzXuAvus4CkDYD1afLIVBbbQvoA7svxaeBF4Kv4kupvAhOBRcBKwG6Svgs8YWbbSloF+AUwzczuavgbKs9/A8dJGm9mT5bKUKUOobu0uDHwPHAg8HU87tB84AC8k3gO8Nd4Z3F/M3uunewP6T2h1xTdqMF98dhXE/BpyyPM7NZ2szekd5zmOqIxFcnNPrjT5A05yowJx0V1t6axjMUXUvwxkVbAg+TNa4VBCZazTdKOwO24Q/ZR+GohAWsDx+NOlJ8C7gN2Bp4DloT4NLOAW8zs4noYF1a/PAdcbWb75Sw+F38IHgyM5AtVjQ6hu7TYF45H45/168Dl+Ge5GG9o/gm4GpgMTGqOqSXpdb0m6UYNTgZOB/4VOBtvyPc31brl6RrNxcZUd7IPsNDMfpejTLEn9qcG2JMLScfjX/Qvm9mlFbKfi49A3Z1IKwCPlFnC3CxK2XY2cK2ZTQWQdCneI5tjZreEtA3xz+FeM1sW0nbFRzMelrRvqOtgM3ukBvsGwvGBYkLWe29mH0i6A+/tlmtM5dUhdJ8W3wQOMLOFoc5f4SNV483stZB2O7BhXY3PT6/rNUm3afAtXIOvhDqvwCOzt5qu0VxsTHUnu5F/iqudemJFAd9fLpOkc4BdgV1TDacCdZjikzSVyg/eiWZ2WxbbwvTjBHyYush7eM/s1ERaP/BQ8SEBYGZ3UP9QJjuEY/I+Z7r3gUHgc5LWNbNXS5yvRofQXVrsA64rNqQCmwBXFRtSibR7azE06rUig5TXa5Ju0+D1NjxY7pbAU7UaFjU3RGxMdRmS1sJXStyYs2jb9MRwsS4C/m+kDJJ+iPfuJ5rZM6nTBbzHUyszgZ9VyLOcj0sZ27YJx/sSaVsD88KDoEgBeDCvsVWwXK+LDPc+QXEFzqdI+DtATTqE7tPieam0fuC0VFof8JOaLI16rcSIei1Bt2kwHbS0n/r4lM4kag6IjaluZMtwzLN6CtqgJxZiQyWD4n0oqfj6kOI8uKRz8eHciWlnUkmrA1tR5osm6R+BI4AdgTWB3wEXAtOTPbrQe63Ug03XPaJtuP/QMqDYCxuL9+rS26v0Af+b57pl7BkNHAt8Dfg4sBBfTXMW3ut6ycwWZr33KYoPyVIBEavVIXSPFtfAQ0I8mEhbBxiXShsHrMMImo16bYpe03StBgP9+N6BI107ai6n5jqiMSVpDEMP5xWATSQVgNebufqlQyjG83k5Z7l26Ik9gMeEORS4C/hV4txtAJIuwJ1I9wX+GIaFARaZ2SJg+/D/w+nK5ftr/QzvET2FOwIvAf4OOAPvAR1arfEZbBvElwOfJOkSYAb+OW0paSszmx/yjwbGS9oIWGxmb1Rpz8rAdcBnwrXPx3+0p+Dfpw0Y6qFXvPcleCkcS8WQqlaH0H1aHEyULeCaeyyVtojUtEvUa1P1mqZrNRga9BtTYmQqaq56zXVEYwpfOpxcwqnwNxs4rBUG1QNJk/CWeFYGzewXFfIU46Pk7VG1vCdmZpeF6aFD8Q1DS017HBOO6Xgxwr8ABWC+lQ6qdy7+kJgGnGpmSwEknYA/nA6RNN3MHq/yLZS1zcyelXQKviLl28CVwJ74l/ku4M9D/lPwTYEnA/+BrwSrhgvwh8RpwNSEo+YsfLUMhOHrjPc+zRvhWGrz1mp1CN2jxT5ci28nzvUDjxa1F+jD/T/SW41EvdI0vabpNg0m30c/8D5QSjdRc1SnuY5oTAXntVGV8nUgk4C/zJF/Nh5DoxzFvYjezWlLO/TEYMjh74FSJ82srA7M7MeU2Cle0s74F/l/zOykVJn3Jc3GY5rsTOmHTEUq2RbyfB9f8ptkl1SeS/Coy1UjaSd8mP46Mzs9Vf8cSU/g/glJJ8qy974ERY2V2v+qWh1CF2vRzH6A79uWTJsKTE2mRb0Oq78Zek3TzRq8GVg5nTdqblj9uTXXEY2pbsXMNm1AtauE43s5y7W8JxbYAe811bKctRTH4Q3yxZKmlDi/XTh2y+bfx4XjGSOcL64kSz4U8t77osZWKXGuWh1C92sxC1Gvw2m0XtP0ogaj5oaTS3OxMdV9vBOOWR4YSVreEwuOgJ8AHjezJXWufq9wPKhCvrwxkdqVvfCHwW9HOL858HszexGqvvdFjb1T4ly1OoTu12IWol6H02i9pulFDUbNDSeX5mJjqoU0yGeq6Cu0ak5z2qEnti1ud9Zh+0xIWhWfa59jZrvXs+52JLzf9YAHk3FYEud3wDcMT66GqeberxaOpTRTrQ6hi7WYhajX5c43Q69pekqDUXPLnc+tudGh4HKVRepHmXniSdTfZ6q4hHSdHPXC0MOjlT4ChXCsd/yQ4v3P4njaDXwQ/tYb4XwxyF7yoVAIxzz3fs1wXFjiXLU6hO7WYhaiXofTDL2m6TUNRs0NJ7fmRkM2p7AiISbKxcGIpcDpZnZ51vKRIRrkM/V0OG6UtYCk1YAVw7+t7IkVf3jrurmomb0j6WFge0n7mdly8VXkWxHcba3fgqZmgrPofHzp8BfMLLmB6IlAcY+ppGNlNff+L8Lx6RLncusQul+LWYh6bYleP6IXNRg1V7vmqpnmWwpMMrPBECPifkk3pJb/RlrHU3iLe9McZcYmXreyJ1YU7hmStgPeBh6rU2P9BOB64EpJN+NxqFbAH7ADwEpmtkkdrtMunImPZF4p6ef4SNEe+Jz/83jgyGSvq5p7X/xxKhWZuRodQm9oMQtRr83Va5Je1WDUXA2ay+2Vb2Yvm9lgeP0KHv107WQeSR+TtFDSFom0afINPrsSSZdL+lar7QiN2odwQWRlTOJ1K+OqzMFXWLwdjsbQ0Gqtdd+EL+u9Er8338CDx20D3AwcXo/rtAtmdhE+jfwC7lB6KP6A2AWPMPyaJTYgrvLe/xXwIcM3KS3WV40OoQe0mPH6Ua9N1GuKntRg1Fxtmhu1bNmQu5SkXwN/A3zNzP4rkX4QHkp+rJm9n0gfwINXbUcCSTOAdc3s8ETajcAjZnZCTXegTZH0CTzI12Zm9maLbTkHOAr4s+TnVSZ/H0PRcDc3s2cbaF6kS5B0DbCxme0wwvlcOgxlohYjDaGSXhP5ogYjuUmPTO2Ah9zfP5U+gDeEkg2ptYGLgCOTGeV7ox0B/GeqjgKtceZsCmb2CPAM8JVW24L3LFbDN/TMQrsMa0c6BEkrAJ8ErimTLa8OIWox0gAy6rVI1GAkNx/5TIUpubWAk4GZktayof1uBkg4YklaBV9VNs3M7krVuQ8+RHZnIv8GwPpk2KVa0sb4UNuBwNfxL8B8fGPEMcA5+PYyjwL7W4a9+RJ17ouHnZ+Az48eYWa3limarGMA3/jw08CLwFfxfXu+aWbF6KzX4MOFF2Sps1GY2Z2SngE+C8xJnw+RX6cAC8zsGHxuGNwf7vVm2RnpaHbEl1KPGIW4kg4hajHSNEbUa9RgpB4kR6YGcPFciDcWvgggaRQ+YnV/4v9ZwC1WemfuCcD9qdgNBTxQ2rwMNvWF49H4fj0748HWZuGNmZPxnu66+L48WSjWORnfyqEPb4ydk6WwpB2B3+D7A24P3BNsOwU4NZF1LrBTWA3Sas5n5OBrr+IbVx4mqZ+h0bRBW35/sEikFPvhMWnmV8hXTocQtRhpDuX0GjUYqZnkar4B3DP9XUlX4lN9s/Gdk9dkyKt9F3yU6GFJ+4a0g8M0F3jcpJcYTgGfJsyypLIAvAkcYGYLAYLj+oHAeDN7LaTdDmyY7W1SwJczHhCc5pF0Be7Nn4WzgWvN99BC0qX4TtJzzOyWRL6XgJXw5eBll982gZ/gO29PTI++mdkzkq4DPs/w1QrTGmFI2JrgoHCto/HRxjPxBvZkM7unEdeNNIYQCfhQsm0yPqIOoblajDrsTSrpNT4PI/Ug2Zj6aPQJ93W4XdKaeCPrPcJeNGZ2B+VXAa7G8kHRCmSY4gv04RsPJuvYBLiq2JBKpN2bo87riw2pwJb48u2yhCnKCcDERPJ7+D04NZW9uE1By0emzGyxpFOBE/ERtTQHA/8OfAF4FjjNzLL4E+RC0p74dOgA/uC4Co/XcSB+vy6StHupKLSRtuWfgafN7MZKGTPoEJqgxajDniaLXuPzMFIT6cbU1eH1b4HfA3+PT2s9amZZNyx9FfhYKq0AnJuxfAE4L5XWD5yWSuvDe71Z6zy/RJ2DGcpuE473JdK2BuaFhmWSYoiIP2S0q9H8FDhE0gQz+03yRPCHy+0sL2kqQ9FhR2Kimd0WXvfjKz4XATPCaOZ5ZvZEqO9JfMq2Xe5ZpAySVgS+A3w5R7ERdQjVaTHqMJKFrHqNz8NIrawAIGkzvCFwP0BoFV+FT/UNcz7PwIP4HjaEulcHtiLDSj5JawBbJPNKWgd3CEymjcNb81XVGcjamFoLd6j/INQ3Fv/yLC6RdzvgxdSoWssIn+NXAIWh7nowE29glvubm8g/D9gbN2J33P9tsqQNw2cznqHduSPtzyTg0sS0fkWiDiMtZBI59ZqTmUQdRhgamSo6nz+UOHcFcBM+pXVZjjp/CUyXtE6Ylts+pD+coWwx72AirQAsAR5LpS0iwzRdqTpDA21jsjWmBvF9i06SdAkwA3gZ2FLSVimHxgn4+28bzGyBpOn4yF56dK+a+l7FRx+zci2wj6Tn8Onf/fH7NBfX3OTo5NkZBOfcrc3syIqZU0QdRppNLXrNStRhpEiyMfW4mb2bOHcn8AY+v5t512oze0TSXHwO+AK84TPfzD4ayZF0GL5qcDMzW5Ao3hfyJrem6cenGZem8j2UFF2GOpORbPuB94HHK5U3s2clnQIcD3wb9yfbE7gOuAtfblvchfof8GXgbYWZ/VLSGpI2Td2bZlx7GR64MckCfH/HSGexN+7nURVRh5EmU5NeG0HUYfcyLAJ6vZC0N+4jtW2pFXyShLfI+1KNpFquWVOddSh/LPBFM9srb9lIJBKJRCKdS+69+bIQVk1cgE+llWIf4Nh6NaTqVGet5d/H9+uJRCKRSCTSQzRkZCoSiUQikUikV2jIyFQkEolEIpFIr/D/ZcTX+h+I8P4AAAAASUVORK5CYII=", "text/latex": [ - "$\\displaystyle \\left(\\frac{1}{2 \\alpha}\\right)^{l + m + n + 1.5} \\Gamma\\left(l + \\frac{1}{2}\\right) \\Gamma\\left(m + \\frac{1}{2}\\right) \\Gamma\\left(n + \\frac{1}{2}\\right)$" + "$\\displaystyle \\frac{1}{N^{2}{\\left(l,m,n,\\alpha \\right)}} = \\left(\\int\\limits_{-\\infty}^{\\infty} t^{2 l} e^{- 2 \\alpha t^{2}}\\, dt\\right) \\left(\\int\\limits_{-\\infty}^{\\infty} t^{2 m} e^{- 2 \\alpha t^{2}}\\, dt\\right) \\int\\limits_{-\\infty}^{\\infty} t^{2 n} e^{- 2 \\alpha t^{2}}\\, dt$" ], "text/plain": [ - " l + m + n + 1.5 \n", - "⎛ 1 ⎞ \n", - "⎜───⎟ ⋅Γ(l + 1/2)⋅Γ(m + 1/2)⋅Γ(n + 1/2)\n", - "⎝2⋅α⎠ " + " ∞ ∞ ∞ \n", + " ⌠ ⌠ ⌠ \n", + " ⎮ 2 ⎮ 2 ⎮ 2 \n", + " 1 ⎮ 2⋅l -2⋅α⋅t ⎮ 2⋅m -2⋅α⋅t ⎮ 2⋅n -2⋅α⋅t \n", + "────────────── = ⎮ t ⋅ℯ dt⋅⎮ t ⋅ℯ dt⋅⎮ t ⋅ℯ dt\n", + " 2 ⌡ ⌡ ⌡ \n", + "N (l, m, n, α) -∞ -∞ -∞ " ] }, "metadata": {}, @@ -275,15 +458,16 @@ }, { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAAVCAYAAAC3+8IMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABJ0AAASdAHeZh94AAAKdUlEQVR4nO2de7BXVRXHPxdQYITwUaMzTvGIKA3ieq82mIqYoJEPuKZjk5FUU46PSCYcw7LFMjXKSa42hjU4gpPkYzAqEUseIgIpUISa5PgAjUjikYGPQKA/1j73nnvuOed3zu93fr9zf3C+M3fOvfu5zlpr77X3Xmuf23DgwAEKFChQoECBOPTIm4A8oaqLgc0i8pUqtX81cAUwwCW9ANwsIgtCyt4NvAe8CVwEfBz4H/AnYKqIPF8NGgsUKFAgCbrlTUDOaALWpqmgqrNVdVrC4v8Arnf9nAwsAear6qcCbTYAFwLzgVHAz4HPAJ8F3gcWqerRaegsUKBAgSwRurNQ1TnAWGCgiLxdW5Kygao2A2uAb4jIrJD8jwJHktJYpIGI/DaQ9D1VvRI4FVjvSz8F6Ak8LSLnBuicALwFnAb8vlKaomSrqpOB24HLRGRupf10dZTSjwrarfuxkycKueSLOP53MhaqegowAZgSmEwuBs4EGoHhQF/gfhH5cgICjgX+ia2YpwEtwHnAMOB4YA/wHHAvcK+I7E/7kkGIyFpVnQ/8UFUfEJHdgSLNwH5gXaV9JYGqdgcuAfoAKwPZ44EFIvJ+SNW+2A5wZ6C9JM6mo0TkP746obJ1aHbPPydoNzOUq1eV6lQC/UjN4xL8rSsUcskPec61cfwP21ncAvwXmBlI/74jfDd2vPKJUoT7MA6b8B7BJsyZwBZgKfA6cCx2Tj8LGKuql4hIFp73HwHPAJOAWwN5zcBLYcqYJVR1GLAK6IXxrkVEngsUGwfcGNHEHZhBWxXVRUz37wX+jpIt2FHZbuClmPaqgXL1KguditMPP5LyOI6/9YZCLvkh77k2lP8djIWqDgFGA7NE5N1AA5Md4S9jVm9pihdoAbYDTwHe+fwCv1VT1RuAZ4EvuJeZl6L9UIjIs6q6AbhCVacHrGgif4Wj6wZfUk/ggKpO8aWNFZHlEU38HVsh9AMuBuao6ijPYa2qg4FBwB9C+r4dOB04XUT2RbzjtFLv4NqKlK2qHoE51Fem3dWp6kRslXKWiDyZpq5DuXpVsU6V0A9/uWmliCkxdmqOQi5tfdVcLjnyHqrI/+DO4muugweDFIhIG8GqcQa9I1T1A5ijdq6b8JaElRORf7mIoFswJ28nY6GqhwFXAxOxyW078DDmRO4BvAYsEpHLfNUewLZjY+g4ITcBNyd4hbuBh3x//xjYDNzpS9scVVlE9mBCB1jrtsSTga+7tPHA4uD2WFVnAF/EFO7VBHSWQqRsMWPWjRofQUF5epWlThGtH2kRx19U9SxH50+BucAPgJFAd5d+paP3RGxlORpbmCwHrhKR1yugLTUOIbmMAf4ITAfuxxaGZwNHYH7FySLyTAX9p0ZXmGsJ4X/QWIwG9mHhmlnhfOBw4DcJyu51z05n9y4a6HHMGfwo9gLnA9dik/V+4GhAAlVXuGfbS6vqQFe25OQoIjuAHT46dgE7ROTl6Fqx6IZNAh7GAXP8BVT1DuBSzFBsKLOfIOJk2+SeVXP2Z4xMdMqhk36UiVJjx+PxEOBpYCFwD/A5bMHQU1VnAr/GBvkcbCCfB9znfu/qqEe5nOSeQ4DVwBMY7z8JfB74naoOFpFdFdBQK1SV/23Gwh1FNAIvZuwEagHexqx3JFS1B+Ddd3g8pMiDmKH4tojc6erchm3XxgInArNDJvHV7jnSl+Y5c/eq6lBf+j4ReTH+dZJDVacDC4A3MCfVl2ifAFDVDwEjsOMpr85dmDNuPLBTVY9zWbvL9a8kkG0uzu0KkJVOQbh+pELCseMZi08DI0Rkvat7E3aWfA42cY0RkVUu73BsVzpSVXuJSNAH1dVQz3I5AzhNRNrGgKrOw45pGrEdXldHVfnvv2dxPLYl3lIWmeFE9cJWTgsTKPp0YCjwmIh0WEmo6mhshbAc+JmXLiLbgI3Y1uso4KZgoyLyFubs+ogv2ZscV2CRAd7PQ2SL44BfYX6LxZixGysiC13+BcBqEXnTV+cqzLAsxmTh/fh9JGlRSrZNwLtAZoayWshKpzxE6EdaJBk73qQ00TMUrv9dmA53B67zDIXL24PpTgN2LNJlUcdy8XYWX/UbCgdvPPSqgIaaoBb89x9DHeOeHUI0K8QYLFQ0dlukqpOA7wAbsFV1EF5aa4jn3mPML0TkjYgudmBRAACIyFRgajzp4RCRiRmWHYddxPPXaUhNVGlEytYp2QnAmignuq/sRqB/RPbSkPPVOWn4lRBZ6ZQfHfSjDMSOHbfCHQK8KiJhK7n+joawc/X+wC4R2R7VeSGXSJSSSx/gY8Am4LGQIoPc85WoDroI76EG/PcbCy9SIEsrehEW19vp8xYeVPUaLDz0b8DZzkcQxJnYGVvUlukd4kPsetP+fl0JK7Az6mojTrbDMT1IcgTVil1k9KORdr/LxkDeumTkpUJWOuVHpfpRauwMx3bxi4IZqjoA2xU/IiJ7A3l9gMGYjyMOrRRyCUMpuTRiu7YnIkL1m7ALsa/F9NFK/ryHGvDfbyy2uucxZAB3Ce0CYInb0oSVuRaYATyPEb81pExvbCv0ioi8E8gbhMUgrwwc5fjLdMOEGSfwXCAiP6lRV3GyTezcFpHWYJoLExyH+YueLI+8ZMhKpwLls9CPUmMnjsfNMXknYZNZrCEv5BKJsuWiqn2x3eCyuDtfefPe9VcT/vt9FluAf2MhqVlgJCak0G2Rql6PEb8Oi/qJIr43NmDCYq1nYJFFUR59sPdpoEY3tbso4mRbT87trHTKjyz0o9TY8SalNSF5zTF53nl6V5dNvcrF428U7xuojwjBmvC/zVg46/kU8EF3UaxStGATfPD7SKjqjZiTZS1m5bbFtLMTu8k4WH0f4FP7xtKF7s8jY+qPcM80F1sOKpSQbRP2ddsXak5YemSlU35UrB8Jxk4TdkQQ9uXguJ2FZ2S6urGod7msD8mrp0VUTfgfvGcxD7vVdy7tF8m8Tsdj4ZxgUT4Ap6rqbPf7NhGZ4so2uLKrgsdDqno5FrW0D4tumhTiBNooIrPBBO76uAb7+urDrv8WjDn9gFHuksk9IrI60NY5rq9OjDzE0Em2LjRzKPDX4Hl5rZBUr7LUqQCy0o/QsaOqPbGw7vUuuimIZmBThAM7tyi1Qi6xRryqyGuuDaAT/8OMxVYsBveuQF4jcHkgbRDtEQObaA/vPBn4MOb8CWKge3bHLtSFYRkw2/f3ddjq91Lgm9huYwbwXWy7eB/2fyM63ERU1X4YIx+NiZQ6VBAm22HAYeS7emokmV5lrVNZ60fU2BmK8TjsXLw/dnywLCSvJxaltrZUlFqV0MjBL5cehB9BgRmLPL6VBvnOtZH8bwj+pzxVnYpFFjWJyF8iOoiFqt6KhaYOEpHcHMuq+i3ssxxniEipiJKDHlnINi9UQ6ey1o965m+5KOSSL2rJ/7B/fjQDu1Ha6YJbCrRgRxt5GoreGBPnFYaiDVnINi9kqlNV0o965m+5KOSSL2rG/07Gwt3+mwCscZeJUkNEThCRxnLqZogBwC+p7ObzQYUsZJsXqqBTA8hYP+qZv+WikEu+qCX//w8lQz8AqrJJcgAAAABJRU5ErkJggg==", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjEAAAAhCAYAAAAoJvpGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAAQ2klEQVR4nO2dfbRVxXXAfyAJIOIXSc0HCiLRaNF34Zkuo4FIBRPMB2BwmTQ1wcQmJSYGrC4rxm53NCn5QDBdqbY19ZFUm1SxJBFjKmrwA9RI81D8KtVqUiNJEEWIIaB9/WPP4c0779x7z7137j3vvsxvLdbROXNm5s3suWefPXvPDOnp6SESiUQikUik3RhWdAPyoKrTgAuATuAtwNki0lVooyIDAlVdAnSKyMyi2xKJRCKR1tIWSgywH7AJ+Lb7F4kklIDugtvQVqjqucCngfEu6VHgChFZnZH3GmAX8CvgdOAo4PfA/cDFIrKpFW2ORCKRLIYW3YA8iMitIrJYRG4C/q/o9kQGFCXgZ9UyqWqXql7W9NYUWHcN9fwvcBEwBTgeuBNYparHpcobAnwQWAWcDPw9cCLwp8CrwBpVPThQ8yORSKRm2sUSE6kRVV0BzAIOF5HfFt2eelDVTuAh4C9E5NqM+28CDiGgJUZVxwK/AD6MWStOADYDZ2IWwSuxF/8mYJ6I/LwJdc8BFgBTgS3AOSJyV6h6ROT7qaRLVHUB8E7gYS/9HcBw4F4ReU+qrWcB24GTgB822qZK8qqqi7B+/6iI3NBoXQOZajLfQLlt/3tQJHFciqVS/0clZhCiqu8AzgIu8CeGqs4D3o1ZLzqA0cD1IvLnOco8BPgl9jV+GTAXeB9wLPBWYDfwCHAdcJ2INGwxE5ENqroKuFxVvysiO1NZSsDvgCcbrcujw10XAAJsA24EuoBXgMXADuDfgfOBhU2o+3zgcuCzwFLsBT45YD17UdV9gDMwBW1d6vYcYLWIvJrx6GjMkvtiRpl5ogUOEpGXXP5MefXodNf/zFFuEOqdK43Okxwy34z+bRuK+g2L41Ls+6NS/0clZnDyJeBl4OpU+hcw4duJLSm8vYYyZ2MvrZuxl97VwPPAXcDPMYvI6cC1wCxVPUNEQoS+/S3wAHAe8OXUvRLwiIi8FqAev8ztwJki8isAVb0ds8y8XURecGlrgTcHrDep+2VX9xZXz01YHwRFVY8F1gMjMHmYKyKPpLLNBi4tU8RVmAVsfaVqKtzb5f13OXlNmOLa+F8VygtNvXMlxDypJPM+ofq3nSjyN+wPfVyKfn9k9n9UYgYZqnokMAO4VkR+l7q9CBO+/8Y06lqWKOYCLwB3A4mvxGpfY1bVxcCDwIcwgVxZ55+xFxF5UFWfAD6tqktSGnqJMktJri2LvaThQI+qXuClzRKRe1KPdgC3JAqM4zDg5kSB8dJ+2oS6VycKjGMiNl4h6wGzXpWAA4B5wApVPTlx1FXVicAE4McZ9V4JvAt4VyUFUkQuK3fPK6uSvKKqozBn4nW1WPdUdT72VTddRH6S9zmPeudKw/Okisz7+S6r1phq/dtqChwXaHBsBsO4NNj/hb4/yvV/VGIGH5/AhOR76Ru+X4VqpY+Fvqjq/pgz5w3upXVnVj4R2eKiWb6EOYL2U2JU9XXAucB87OX0ArZccxEmj/8DrBGRj3qPfRczQc6k70u1hFkEsrgG+Dfv/78CPAd8w0t7LuO5UioP2FLO36TSOoB/bELdf5dRd3fgehCR3fQqRxucaXsR8EmXNge4I8M/ZRlmlZouIk+XaVctlJVXRwn7gmvZUhLUN1dCzhPKy3ytVOxfVZ3u2rkUuAGT82nAPi59gWvvMdiX+AxMWb4H+ExIn7A8DIDfsFaNy0zgP4AlwPXYx8opwCjMb22RiDzQQP01MwD6HjL6vy2UGFXdD/siBftBO0xVS8C2Vk+iNmAG8BoWAhuK9wOvx/xAqrHHXfv5UbhIltswh9FbMCF8P+ZX8hwWeXYw5ovic5+77hVcVd0XeBtlIpNEZBvmz5LUvQOTl0yrhsszCjjCL1NVxwCHptIOBcY0u27HZMwEG6SeCgzFXk4Js4EVqTZehTk4TxeRJ+qoI4tq8jrFXTcEqq+ZBJknjn4yXyd5+/dI4F7gR8C3gPdiiuxwVb0a+Ffs5bMCe8G8D9vu4uQG2tZKQo1Nq8Yl8YE7ErP43o71/R8DpwE/UNWJIrKjgTa0iqbOi7ZQYrBoEN90pe7fCuyLvi1R1YXAgTU80i0iqyqUNwr7cn08sKPYXOC32JdBWVR1GPAx97+3ZWT5HqbAfF5EvuGe+RpmopwFHAN0ZbyEk2WbaV5aEg78MOFIyuz20krYviiPptJ2UmaZJ1TdToEaS+B9cNwGgauxSKjRwJ/R+2JCVd+IRWXN8575JuaEOAd40UWGAezMcnLM2Y488tpyp94GCDVPIFvmayJn/yZKzJ8AJ4jIw+7ZL2K+CqdiL9SZIrLe3Xs9JvvTVHWEiOzqX+yAI9TYtHpcpgIniche+VfVldhySwmziA10mjov2kKJcWt3Q4puRxNYCIyrIf8KbM+OcrwVMwM/X3+T+qKqI7Cvsh/l+LFaAkwCbhWRPl8pqjoD+/q4B2/JRES2quozmLnx98AX04WKyHZV3YX5oSSUgM0i8kqtf1MFOlyZ/g/LZGBTKkKnA9gYIgIro25fIZiMfZk8FrAegDcB/+Ku2zFFcJY3Zh8AfpryC/qMu96RKksx82495JHXKVgE2uN11tESQs2ThDIyXyt5+xdgfqLAuPp3uHlZAi5MFBh3b7eqPolZKEfR11l1wBFybFo4Lokl5mxfgXEkc2FEA21oCa2YF22hxAxWRGR84CLHuGu/sNcGmImF31Y0BarqecBfAU9gX+xpkrTlGV7niXD/g4j8okwV2zAPdgBE5BrMJyQXIjI/R55+ZYrI14Gvp9KuAK5oQd1rMDNssHpy5ptNSlkWkWZ8RFSUV/cDeDTwUCUHYveyLfcxcFfG+v2KvH1VA6HmiU8fma+Dav07CluueFpEsr58x7k2ZPltjAN2pJzd0+U/Q/HjAuHHptnjsh+2VP4scGtGlgnu+lSlSgZI/zd9XkQlZnCReLmH1NBPx2L4+21Jn6Cqn8UcbB8DTnG+GmnejVkVypkJX6Fy2OJIev++SHO5D/OBaDbV5LUD+42qtpS0nP7LsiV6/XqeSd3rzte8mgg1T3walfk8/TsUWJO+oarjgYOwqLw9qXuJj+K9VepfTvHjAuHHptnjUsJWHm7P+OADs55tx4IgKrGc4vu/6fMiKjEFEtonBvi1u46pkCc3ahuhfQC4U0S2l8mzEFiG7WB7ioj8OiPPSMz891R6+UdVJ2D7DaxLLV/4eYZi/VRt0kYCICJfbVFV1eQ1l1OviCxPp7lQ0tmYj9VP6mtePkLNk1T+EDLfSP92Vrg3GXvJVlQuix4XV1/QsSl6XFR1NGY9W1tGwdlL0f3fqnkxzN0IsSlZpAwVTPELCesT8zzwGyx0OQTTsImWaQpU1YuwdcxuzPFva5lyRmI/elk+JMuwqJhy3uhgf88Q4kGPg41q8touTr2h5olPCJmv1r/Jy/KhjHudFe4l/hoDfVwg/Ni0YlyS/i3X90Noj2i9lsyLYVDbercLL/0O8EfYi+dyEbkx7/ORXkL7xIhIj6reDXzIhd81Gj0zF1M80mftoKqXYk64G4BTq5gAX8SieSaq6nFeBMQCbNMjqGyROsFdg50fFCmeHPI6hf6RYQORUPPEp2GZz9m/u7Gv4DSVLDGJ8tMOSkzosWnluGRFXraLYg8tmhf1LCe9CiwUkW4XZrlBVW8NHNIbqZ+V2I6H7yEVAqyqc7AQWbDIFIB3qmqX+++tInKByzvE5V2fXuZR1Y9jAvgaFm10XoaT2DMi0gV7J20XdhbQGlW90dU/FxPwA4CT3UZH3xKR9E64p7q6+k2GSNuTKa8ujHcSFgW2p8yzTSPvXAk5T1KEkvly/Tsc29LgYbGND9N0As+WcdwtLGKsqN8wj4EwLlCAJWYA9D1k9H/NSoyIPI8LDRPbYW8rtkGZf9DgQZiX8Yki8pRLWwJ0isjMWutsB9yL+X4RWVpwU1Zia64fA76ZulcCPp5Km0Cvt/uzQLJl/fFYCOXyjDoOd9d9KH8A4lrs0MSEC7Gv6jOBT2HWmWXAX2Mm0m9jp0b32aVRVQ/AJsMtFSKXIu1LOXk9FngdxX1xlsg3V0LPk9AyX65/J2H9m+V3MQ5bBlibcW84FjG2oVLEWBMpUdBvWAvHZRjZS0lgSkyrzxFLKFHc+6Ns/w/p6enxM92B7dfxSRH5Zy/9I9h5C6P9ryK147FXiMikVGVfA94gImd7abdhh/VdWKbRbY3agXprsSPVM52YWtiWi7FInykikrmrbI4yvgxcDEwQkcIcalX1c9g2+lNFpFo0RKQNCSGvRdGMeRJa5tu5fxsh9NjEcclPK+fF0FS+Kdhx2fNS6Z2YAuIrMAdjX8+fSlW0L3AOtnW1T4ky27QPBsROAH4aqHoseQtYhu222W/juBqYi5nyi1RgRmITYWVUYAY1IeS1KILOkybJfDv3byMEG5s4LjXTsnkxzMt0BOZcuRhYrqoHishL7nYnntnRmRRXAUtEZF2qvtOAHnrPOMD5zhxCDo9uVR2LbYf+YWx54QRgM7YMsR9wJWaq2gTMkxxnJ3llzgEWYFs5bwHOEe9QqypldAJfBU7Ezvn5BLYu+HkROcll+wHwEfov47QUEdmlqmcB01V1VD3+SiJydBOaVivjsUMWu4ptRqSZhJDXomjCPBlPYJlv5/5thMBjM544Lrlp5bzwfWI6Mafd6zD/hdnACuegMwU74TRx2OnCYr+/k1HZVGy91A/bLmGOYE/maGyHuy7ADgLchp1y3IVtiLYY2IGFbZ1P+TW1rDLPBy7HHEyXYgrR5HIPJaid8LsWM/39pWuXYr5Ai7ysDwJfUNWRUvCx9yJyN3bsedsiIo9T/5b2kTZiMMhrCJol87F/GyOOS7FU6n9/OakTeFTsfIOV9C4pTQT2p9fB7iTMKjJHVbvdv2O9csZhS1I+JWw5Ko8jWAnbjfBMEVnrlmlux074PUNE7hWRjZhS8eYc5SVlvuzKXONC2m4C3pjz+aXAD0XkChHZjCl007ATg/1jxX+JOcu9JWe5kUgkEolE6sS3xEyhd8loJbBWVffHlJvdwCMAbj0q7UvjMxJI77xaIv/mQB2Y97FfxmHY9tcvpNLSobiVylwtIlu8tInkOIXYLYVNBaZ7ybuxPrg0lT2xvozM2a5IJBKJRCJ14isjvhLzABb+9UGXvqlMzHoWW7EzN3xK5FdiSsD9qbTJwPpUWgf5HYVLGc9PztmmZG3PD3k7Cngyw8HrYHf9Tc52RSKRSCQSqZOhAKp6OPYC3gC2ORlwM7ak1MepNwc/wzbrwZW9L3YiZ1WFQ+1U1SP8vKo6Bos399MOxfYxqKtMR14l5kDMUfk1V95o4BLMPyfNJOC59OY+kUgkEolEwpNYYhKn3o3evZuw3QSPpzYl5sfA0U75ADjOXbO2UE6T5O320kr033q8hG34k2db/X5luraNJZ8S042d1XCxqh4FXI9t9neEqr4tlXcq9vdHIpFIJBJpMr4S85hz6k24D3iJvk69VXGOuA9iIdJgCsdm8U4vVtX5qtqjdty7T4fL64eaTcaWs15N5dsoInsPFMxR5s5UmXuwo78rPu/i3C/BoqU2YpFRM7AQ73Xe8yOw2Ph/SvdJJBKJRCKR8PTZsTcUqvpe4CrgmKyIJLWDEuYBHSnlpJE6GyozwPPnArNF5NRan41EIpFIJFI7laKM6kZEbsM2fBtbJstpwLmhFJhAZTb6/B7gc3U+G4lEIpFIpEaaYomJRCKRSCQSaTZNscREIpFIJBKJNJv/B0hP3uB3BqUPAAAAAElFTkSuQmCC", "text/latex": [ - "$\\displaystyle \\left(\\frac{1}{2 \\alpha}\\right)^{L + \\frac{3}{2}} \\Gamma\\left(l + \\frac{1}{2}\\right) \\Gamma\\left(m + \\frac{1}{2}\\right) \\Gamma\\left(n + \\frac{1}{2}\\right)$" + "$\\displaystyle \\frac{1}{N^{2}{\\left(l,m,n,\\alpha \\right)}} = \\left(\\frac{1}{2 \\alpha}\\right)^{l + m + n + \\frac{3}{2}} \\Gamma\\left(l + \\frac{1}{2}\\right) \\Gamma\\left(m + \\frac{1}{2}\\right) \\Gamma\\left(n + \\frac{1}{2}\\right)$" ], "text/plain": [ - " L + 3/2 \n", - "⎛ 1 ⎞ \n", - "⎜───⎟ ⋅Γ(l + 1/2)⋅Γ(m + 1/2)⋅Γ(n + 1/2)\n", - "⎝2⋅α⎠ " + " l + m + n + 3/2 \n", + " 1 ⎛ 1 ⎞ \n", + "────────────── = ⎜───⎟ ⋅Γ(l + 1/2)⋅Γ(m + 1/2)⋅Γ(n + 1/2)\n", + " 2 ⎝2⋅α⎠ \n", + "N (l, m, n, α) " ] }, "metadata": {}, @@ -291,87 +475,409 @@ } ], "source": [ + "N = Function(\"N\")\n", + "a = Symbol(\"alpha\", positive=True, real=True)\n", "l,m,n = symbols(\"l m n\", integer=True, nonnegative=True)\n", - "alpha = Symbol(\"alpha\", positive=True, real=True)\n", "\n", - "def current_squared(l,m,n):\n", - " L = l+m+n\n", - " N = pi**1.5 / (2 * alpha) ** (L + 1.5)\n", - " N *= fac2(2*l-1) * fac2(2*m-1) * fac2(2*n-1) / 2**L\n", - " return N\n", + "H_l = H_na.subs({k:l, eta:2*a})\n", + "H_m = H_na.subs({k:m, eta:2*a})\n", + "H_n = H_na.subs({k:n, eta:2*a})\n", "\n", - "def gammabased(l,m,n):\n", - " def fac2pm1(p):\n", - " return fac2(2*p-1)\n", - " \n", - " L = l+m+n\n", - " N = pi**1.5 / (2 * alpha) ** (L + 1.5)\n", - " N *= fac2pm1(l) * fac2pm1(m) * fac2pm1(n) / 2**L\n", - " return N\n", - "\n", - "\n", - "display(current_squared(l,m,n).simplify(), Z2_expanded)" + "N_val = (H_l * H_m * H_n)\n", + "display(Eq(N(l,m,n,a)**(-2), N_val))\n", + "display(Eq(N(l,m,n,a)**(-2), N_val.doit().simplify()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Next steps" + "We might finally ask - what is $H$ for odd arguments? And we can see: the integrand is odd, so it is zero.\n", + "And just in case, let's check: " + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAMgAAAAkCAYAAADM3nVnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAAI/0lEQVR4nO2debAcVRWHvxACCImkQAKiQQuEhBSSF6IBDVsKiaCoVEQFJCzlUgoFhgDFpvzyE5AELEhCgSwuCYhlQYVFEKsQFKMCiiRhkbAqCaBsCVsWKZbnH+d20m94M2/mpefNvNBf1avbc/tO95m+ffrcc869/QZ0dnZSUlLSPRsWfUDb04DDgAXAd4E9gPOA1cBUSfcUfc6SkmaxQZEHs/0ZYBtgLKEg1wMzgMPT3/m2BxR5zpKSZlK3gtgeYHtgD83GAHMlrZB0AbAxMFvSYklPAY8AH+i1tCUlfUyPQyzbmwACvgMMsn0pcLqkt7tp/ihwAHC37X0IBZlq+1bgNWAksKwo4UuKxfZw4GpgGPAWcLak61orVWupx4JcCZwGDAU2A04BTq3S9mZga9tLgR8DhwDTgb8DDwEXSXpnHWUuaR5vAVMkjQImAjNtb9ZimVrKgFpRLNvDgOeAAcAJafta4BlJw/tEwoKwPR0YK2n/VsvSX7B9P3CQpKdbLUur6GmI1UEoB8A1wKvArcCLtjeVtKqJstWF7dOBScAI4A3gHmII+FBF0w5gUZNl2Rs4mQhSbAscI2lOM8+5Lti+ANhV0me72TcWGPheVg7oWUEyh/p/kpan7c83UZ7esC9wKXAvocw/BG63PSonM4SCXNXowW3PAZ6SNK2O5oOJoeRVvTlXCxgHzK+stL0FIf+3+lyiNqMnBRmSyteaLUhvqXz62Z5MWLrxhE+E7W2ArUkWJI2rfwbsBExKEbYiZLmVsLCZYrUM24uAm4AtiLzUO8DFks62vRGwAhgE7G37+8BiSaNsbwzcCEyXdFdLhG8jelKQwalc0WxBCmQIEXx4OVfXQSQqH7U9gsjPLATGS1rd5xI2GduDgJ2JazENmAVMjl2+EngB+BTwD2B3YCnwRspRzQH+IOnqgmQZno5/g6RJRRyzL6nXgrzebEFqYftE4ELg65J+1UPzWYSluDtX1wE8CBwMXEGEL2cWLWcbMQrYCDhJ0o0Atn8KnAUMlfSc7Q8S/XqvpM7UZk/ga8ADtg9Ox5os6cF1kGVsKhdkFQ32Z1Ow/WFiOH4AsCXwX8JyWtKah2t/sSDZRb6vViPbFwJ7AntW5Gk6gB2BnwNflPSnGsc4AzgjV7Ux0Gn75FzdgZL+XL/464btc4Aze2g2QdKdaXs0sJw0xExslcrnUjkGuD9TDgBJf6Hg2RXAbqnM911d/dksbO8A3EXke24iEtjjgO8BB9geL2kZ9BMLQlzQFcBj1RrYvgg4lLhR/lWxu4MYVh1OjMlrcRkRys6YATwLzM7VPVuX1MUxE/hlD22W5rY7gPu6eUgskfRK7vPCQqSrzbssCHX0Z5O5lFCOEyRdnFWmB+yJwLlEYry9LUjKXeSTku/YzraPzMbJtmcRQ4MJkh6pOMamhPWYDPwRuNr23pLyHbaGFPlanvv+68BySU8U86u6yPZl4JvAJ4H3A0uAXwAz8je3pJeAlxo49GjCv8gzhq5h7tHA7xqX+t3Y3hA4DvgGEfh4nrgJzycsyH8kPV9vfzaTZD0mAk8Bl1TsFvBtYLLtkyStbHcLsgCYCxxFmMTf5/bdCWD7EuLmPxh4OUWsAFZIWgHsCnQCD0m61/ZI4Gbb4yQVaglsDwY+lj5uAGxnu4NQsKW5dgMJi3Ao8ARwHZHDOZB4eo1Iv7m3jCZmQOQZQ9frtyEw0va2wKqcZWmIFBG7BdifUMCLiTH9NOJabMPaoV6P/dkHTEjlbZWzOiS9bvuvhALtAdzR1hZE0rW2hxIXdK6kK7ppdmwq76ioN9FJHcDjuWjVWcQN+BvbexWc7PwEYaXyMpi4KY7O1c8ilGM68ANJbwHYPoW4UY60PUPSw40KkJzPLclZixSdGk080TPOJIaPU4HLiaUJveESQjnOAs7JOfxzgMzXWwB192fl75lCTHOql0VZYKIKI1JZbXj3OKEgO1GHgrTagsBaJ6/akKjm9HlJlxF+Rfa5E/hqvSeXdHQDbe9k7cyDbrG9O6HUN0k6veL7b9qeC3yaCL82rCCSnulGhh2Jh93CXLtriNkRvcb2OGKIeIuksyvkmG97MRFuzjvjNfuzG6YAH2lArLlENKoam6fy1Sr7s/qh0OY+SGI34E0iTLs+cDxxA69Ki8sq2SWVRUaTxgDLmjBt5PhUnltlfzZzO68MDfWnpI/2SrKCaGsLkpy/jwMPS3qjFTI0gYmpPKyHdksKPGelg14UEwkl+FuV/dsDL2S+Xpv0Z2YhNq+yP6t/BdrfgowCNqF+c9zWpLU1WwHzJe3TV+eVdFrRx0y/ZRiwMJ9Lye3fjZiwmY+UNdyfTfBBHk3lTlX275jKx6B+BWmVD9KRyr6I1/cFmW+wPqyqfDv9DauyP0ts5pWhI5WN9OcUivVBsiDKRNsb5CNZtocQc/hWEbPCqyuI7fcB2RLbVlmQLVPZtpMlG0HSatsPALvaniTp+so2abrH3VVWbLYNKaDwOBEq/oKkNVl726cSSxCgq4PecH8W7YNIetL2bcTw8DgiLJ1hYlHg5ZJWQm0LMiS33SoLkl3cc23vAqwE/tnPl4GeAvwWmGf7duABwiH/EJFhHiRpuxbK1wjnEU/sebZ/TUxj2ZfwM54GhtPVgrRLfx5L5GFm294PWExEDScQQ6s103pqRUoG57ZblQeZT0RKVqZSrDXT/RJJtxFh3HnEjXQCkejcGbgdOKZ10jWGpKuIIdAzRNDhKEIxxhPJ2WWSluTat0V/SnqSyFnNIRTjJGAHIj+1RzYPC2osubU9mrWRj+0l/bt5IpeUtCe1LEg7DLFKSlrKGh8kZUWnEctLjyXGjxBvulj+7q+WlKz/5J30l4jJcqvTqrMjUv2i8lU9Je9Vuvggtm8GDqpoc4ikeUWetHx/b0l/odIHmUxMYHsNuB/4UhOUo3x/b0m/oUseJK0JOKL7pt3Ti+Wga97fC1yQ1j7PlrQ4HS97f++LjchRUtIMivj3BzNpbDlo+f7ekn5DzVePNoM0fPoJ8DliaeZXgL2AHxERs6mSbuhToUpKqtDnClJS0p8o+hUvJSXrFf8HDd9X9cHn/0IAAAAASUVORK5CYII=", + "text/latex": [ + "$\\displaystyle \\int\\limits_{-\\infty}^{\\infty} t^{2 k + 1} e^{- \\eta t^{2}}\\, dt = 0$" + ], + "text/plain": [ + "∞ \n", + "⌠ \n", + "⎮ 2 \n", + "⎮ 2⋅k + 1 -η⋅t \n", + "⎮ t ⋅ℯ dt = 0\n", + "⌡ \n", + "-∞ " + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "H_odd = Integral(t**(2*k + 1) * exp(-eta * t**2), (t, -oo, oo))\n", + "display(Eq(H_odd, H_odd.doit()))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can use SymPy to generate a python function that uses the [gamma](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.gamma.html) function from scipy" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Help on function _lambdifygenerated:\n", + "\n", + "_lambdifygenerated(l, m, n, alpha)\n", + " Created with lambdify. Signature:\n", + " \n", + " func(l, m, n, alpha)\n", + " \n", + " Expression:\n", + " \n", + " (1/(2*alpha))**(l + m + n + 3/2)*gamma(l + 1/2)*gamma(m + 1/2)*gamma(n + 1/2)\n", + " \n", + " Source code:\n", + " \n", + " def _lambdifygenerated(l, m, n, alpha):\n", + " return ((1/2)/alpha)**(l + m + n + 3/2)*gamma(l + 1/2)*gamma(m + 1/2)*gamma(n + 1/2)\n", + " \n", + " \n", + " Imported modules:\n", + "\n" + ] + } + ], + "source": [ + "f = lambdify((l,m,n, a), N_val.doit().simplify(), modules=\"scipy\")\n", + "help(f)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "$$\n", + "\\def\\br{\\mathbf{r}}\n", + "\\def\\bA{\\mathbf{A}}\n", + "$$\n", + "\n", + "## Integrals over sets of basis functions\n", + "\n", + "The above visualisation uses the formulas defined earlier to numerically evaluate\n", + "each atomic orbital from a linear combintation of primitive Gaussians. Introducing\n", + "the basis set allows us to replace the problem of solving a system partial differential equations (e.g. the Kohn-Sham equations) with an algebraic system of equations that can be solved using standard matrix eigenvalue techniques.\n", + "The individual elements of the matrices in this system are single integrals over pairs of\n", + "primitives:\n", + "$$\n", + "\\tag{5}\n", + "M_{AB} = \\int p(\\br - \\bA; \\nu_A) ~\\op p(\\br - \\bB; \\nu_B) d\\br.\n", + "$$\n", + "The operator $\\op$ in its most general definition is\n", + "a function that transforms the $\\psi(\\br)$ function to represent a physical observable.\n", + "That sounds fancy but is just a complicated way of saying we are breaking\n", + "down something we can measure (e.g. the energy, dipole moment, etc) into components\n", + "that we can evaluate through these integrals. \n", + "\n", + "## The overlap integral\n", + "\n", + "The simplest operator $\\op$ is simply an \n", + "identity mapping, $\\op p(\\br) := p(\\br)$, used in conjunction with the trivial $g(\\br_1, \\br_2) = 1$, that we name the _overlap integral_:\n", + "$$\n", + "\\tag{6}\n", + "S_{\\mu\\nu} = \\int p(\\br - \\br_\\mu; l_\\mu, m_\\mu, n_\\mu, \\alpha_\\mu) ~ p(\\br - \\br_\\nu; l_\\nu, m_\\nu, n_\\nu, \\alpha_\\nu) ~ d\\br.\n", + "$$\n", + "\n", + "Gaussian basis functions of the form above have a few convenient properties that help\n", + "make it easier to evaluate integrals of the form in equation (5). To show this\n", + "we start by expanding the overlap integral into its full three-dimensional form:\n", + "$$\n", + "\\tag{9}\n", + "\\tilde{S}_{\\mu \\nu} = \\iiint p_\\mu(\\br) p_\\nu(\\br) dx dy dz,\n", + "$$\n", + "and observing that a primitive $p_\\mu(\\br) = p(\\br - \\br_\\mu; l_\\mu, m_\\mu, n_\\mu, \\alpha_\\mu)$ can be written as a product of Cartesian components:\n", + "$$\n", + "\\tag{10}\n", + "p_\\mu(\\br) = N_\\mu\n", + "h(x - x_\\mu; l_\\mu, \\alpha_\\mu)\n", + "h(y - y_\\mu; m_\\mu, \\alpha_\\mu)\n", + "h(z - z_\\mu; n_\\mu, \\alpha_\\mu)\n", + "$$\n", + "where we are using $h$ from (defh), but not assuming that $l,m,n$ are even.\n", + "We directly see that we can separate the three-dimensional integral in equation (9) as a\n", + "product of three one-dimensional integrals:\n", + "$$\n", + "\\tag{11}\n", + "\\tilde{S}_{\\mu \\nu} = N_\\mu N_\\nu \\tilde{S}_{\\mu \\nu}^{(x)} \\tilde{S}_{\\mu \\nu}^{(y)} \\tilde{S}_{\\mu \\nu}^{(z)}\n", + "$$\n", + "where we introduce the one-dimensional overlap as\n", + "$$\n", + "\\tag{12}\n", + "\\def\\twelve{\n", + "\\tilde{S}_{\\mu \\nu}^{(x)} = \\int_{-\\infty}^\\infty \n", + "\\left((t - x_\\mu)^{l_\\mu} e^{-\\alpha_\\mu (t - x_\\mu)^2} \\right)\n", + "\\left((t - x_\\nu)^{l_\\nu} e^{-\\alpha_\\nu (t - x_\\nu)^2} \\right) dt\n", + "}\n", + "\\twelve,\n", + "$$\n", + "and use the same definition for the y and z components. We note that if $\\br_\\mu=\\br_\\nu$, we have the normalization integral from (defH), because we can gather the $t - x_\\mu$ terms\n", + "$$\n", + "\\tilde{S}_{\\mu \\nu}^{(x)} = \\int_{-\\infty}^\\infty \n", + "(t - x_\\mu)^{l_\\mu + l\\nu} e^{-(\\alpha_\\mu + \\alpha_\\nu) (t - x_\\mu)^2} dt\n", + "$$\n", + "and substitite $u = t - x_\\mu$, $du = dt$ to get\n", + "$$\n", + "\\tilde{S}_{\\mu \\nu}^{(x)} = \\int_{-\\infty}^\\infty \n", + "u^{l_\\mu + l\\nu} e^{-(\\alpha_\\mu + \\alpha_\\nu) u^2} du = H(l_\\mu + l_\\nu, \\alpha_\\mu + \\alpha_\\nu)\n", + "$$\n", + "and in particular\n", + "$$\n", + "\\tilde{S}_{\\mu \\mu}^{(x)} = \\int_{-\\infty}^\\infty \n", + "u^{2l_\\mu} e^{-2\\alpha_\\mu u^2} du = H(2l_\\mu, 2\\alpha_\\mu)\n", + "$$\n", + "as before." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "And now, evaluate over the range $(0,1)$ rather than $(-\\infty,\\infty)$:" + "### The hard case: non-coincident centres\n", + "Let us return our attention back to the general overlap elements in equation (12), \n", + "repeated here: \n", + "$$\n", + "(12) = \\twelve\n", + "$$\n", + "but let's replace $\\nu,\\mu$ subscripts with more visually distinguishable symbols:\n", + "$$\n", + "\\tilde{S}((A,\\alpha,i), (B,\\beta,j)) = \\int_{-\\infty}^\\infty \n", + "\\left((t - A)^i e^{-\\alpha (t - A)^2} \\right)\n", + "\\left((t - B)^j e^{-\\beta (t - B)^2} \\right) dt\n", + "$$\n", + "and grouping the polynomial terms and the Gaussian terms we get\n", + "$$\n", + "\\tilde{S}((A,\\alpha,i), (B,\\beta,j)) = \\int_{-\\infty}^\\infty \n", + "\\biggl((t - A)^i (t - B)^j \\biggr)\n", + "\\biggl( e^{-\\alpha (t - A)^2} e^{-\\beta (t - B)^2} \\biggr)\n", + "dt\\\\\n", + "$$\n", + "To make progress we start by rewriting the product of two\n", + "Gaussian functions as a single Gaussian:\n", + "$$\n", + "\\tag{16}\n", + "e^{-\\alpha (t - A)^2} e^{-\\beta (t - B)^2} = \n", + "G ~ e^{-\\gamma (t - C)^2}, \\\\\n", + "\\gamma = \\alpha + \\beta \\\\\n", + "C = \\frac{\\alpha A + \\beta B}{\\gamma}\\\\\n", + "G = \\exp\\left(-\\frac{\\alpha \\beta}{\\gamma} (A - B)^2\\right) \\\\\n", + "$$\n", + "Next we rewrite the product of the leading polynomial terms in terms of the Gaussian product center $C$:\n", + "$$\n", + "\\def\\CA{C_A}\n", + "\\def\\CB{C_B}\n", + "\\def\\tc{t_C}\n", + "\\tag{17}\n", + "(t - A)^i (t - B)^j = (t - C + C - A)^i (t - C + C - B)^j = (\\tc + \\CA)^i (\\tc + \\CB)^j,\n", + "$$\n", + "where we introduce\n", + "$$\n", + "\\tag{18}\n", + "\\tc = t - C\\\\\n", + "\\CA = C - A\\\\\n", + "\\CB = C - B.\n", + "$$\n", + "There isn't a universal term that covers the product of binomial expansions that we have arrived at in equation (17),\n", + "and using the [binomial theorem](https://en.wikipedia.org/wiki/Binomial_theorem) we have:\n", + "$$\n", + "\\tag{19}\n", + "(\\tc + \\CA)^i = \\sum_{k=0}^i \\binom{i}{k} \\tc^{i-k} \\CA^k \\\\ \n", + "(\\tc + \\CB)^j = \\sum_{l=0}^j \\binom{j}{l} \\tc^{j-l} \\CB^l \n", + "$$\n", + "\n", + "$$\n", + "\\tag{20}\n", + "(\\tc + \\CA)^i (\\tc + \\CB)^j = \\sum_{k=0}^i \\sum_{l=0}^j \\binom{i}{k} \\binom{j}{l} \\tc^{i-k} \\tc^{j-l} \\CA^k \\CB^l \n", + "$$\n", + "we can use SymPy to help us take the product of these expansions for some representative values of $i$ and $j$" ] }, { "cell_type": "code", - "execution_count": 49, + "execution_count": 28, "metadata": {}, "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIcAAAAoCAYAAADUrekxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABJ0AAASdAHeZh94AAAHqElEQVR4nO2cf5BVZRnHPytY/lgaxyGgGKoxlYGx2Ni0GSGECjRTY8gZwwkLpEYslBxIMfXrV1HBQpdxdEomJYlqMqwBlWLUSItSskC3JmKoCC0IxDJHS4Xtj+ddOHvZu/fu7rl7uXK/M2fec877vud57nuf8/x8721oa2sjL9geB8wFmoF3AtMlLcuNQB19iiNyfl4j0ApcDrya87Pr6GM05Kk5srD9MvCluuaoXeStOep4E6EuHHUURdnCYbvBdr9KMlPHoYX+pQbYPgoQcAlwpO27gPmS9laaucMRtocBy4FBwBvAjZLurwYv5WiOpcBVwHHAscA84MoK8nS44w1gjqSRwCSgxfax1WCky2jF9iBgB9AAXJbOfwA8J2lYJ+MbgRPT5XpgIbAK2CPpb/myfhDthUCzpImVpNPXsL0JOEfS9r6mXcqsNBGCAbAC+DfwMLDL9jGSXikY/0HgZ5lrp+PbwOd6yqTt+cAUYDjwP+DXhGlrLeB1Y09pHIqw3Qz0q4ZgQGnhGJja/0rak84/UWywpHUcEKY8MR64C9iQnn8D8IjtkRm+moD7KkC7KrB9PPF5Pl8tHkoJx4DUvlRpRrqCpDOz17anEVpsDLDa9hBgMElzJBv9LeBkYIqkv/YVr+ltvxU4HXgemAEMAS6XNMb2+YQWPlnStjRnCXAOcLqknbbfCvwYWChpfV/xXohSDmljal+uNCPdxACC9xfTdRORrt9sezjwFOHYjeljwTgVeIIwre8nzJ+BrwLXpmErgWeBa9KcucBU4KwkGA3AMuAxSctz4muY7TbbD3RnXrma4z89Y6tiWEJoiV+l6yZiwScDdxPhX0sV+FoMrJa0AMD2d4HVwOOSHgOQ1Gb7auAh21uBq4GPStqSnjEGuAB4xvbkdG+apGd7wVdzan/bnUmlhOOQ0xy2bwPGAmMzuZYm4CTgHuA8ST/PidYC4q3vChMkrUum7cPAhEzfa4SGuzY7QdJa2xuABcC5kjZk+n5B/pnr0al9ujuTakpz2L4d+DTxhfw509UEPABcCByfI8kW4DslxrSH6CNS+5tM33Bgc/rC98P2R4BRhHO9s/dslkTtaA7ba4GJwPmSVmbuNwD3Ap8FFkm6KtO3hFC3EyT9MXP/GEJrTCNs/XLb4yQVXQjbnwJmAqcCbwO2JbqLsplfSbuB3WV+rOOANmBvojGA0Do7CmiPAn4EzCYiv1uADg53T2C7P/BF4GLCEd9JRHi3Eprj75J2prFlrX8p9VUpzTEP2AfcWFCv+Xpi7O4CwbgTmE5ohhdtD0lHI+H4tQGtklYAtxMRzNBCorb72f4e8EMiWXc/sYD7gJsIs9RTbCQ0wfzkFK8A/gG81/ZJif67gTXAYkn3EGWJibbH94Iutt9C5J9aCOG8A3gUuJ7wwYbQ0aSUtf5ViVYkbSLqByOIN57kpF1BZGBnFUy5lBDUR4kFbz/mEiZli6T2zUXXAb8EViWtksUSwiwtBEZImiVpDnAKkdG9yPbIHn6mvxCaYhawiXihPkZsflqf8hY/IRzWG9KcVkJAb+kJzQzuJDTBdcBoSfMkzSA00sw0Zr8mLXf9S6XPNxBZz5sllXLMuoVUYPoToXYXE9L+U8KhfC1PWoneh4joZpWkyZ30fwH4JjBD0r15068UbJ8GPAk8KOncTvr/QAjBeZJWZ+6XXP+qRSuStttuIYp6dxBv7pRKCEbCbELtv2L7+k76T0ltre1xmZ3am4r0v5DaDj5YOetf7WhlV+b84k5qNXliUmqnlhi3rYI8VAKTCAF4skj/CcA/JT3fSV+X6181zWH7QsIB2kFKL3Owr5EXraOAtxPJqDMqQaMaSJ9rEPA7SQf5B7ZHE78CWNNJX8n1L9chzVVz2D6bSBG3EtHGZmBm8vIrgfZi4MAuR9Ue9qZjUJH+dj+xg0kpd/2LCofto4H2MCc3zWF7LBFKPgecKWkXUWfoDyzKi04WKZJ5Bhhpe0oxvmptG6Sk14EtwFDbHZxR21cS2xwgE8Z2Z/2LRitpo0979m5MHtVB203AOqJINlbS1kxfe2Q0TtITvaXVCe1JwEPEIjxCCMsRwFAig3ikpHflTbfSsH0RsV/mdeD7hJkYD7yP8CmGAe+RtK2769+VWWnMnPdac9g+kYjz2wiJ3VowZH5qv9ZbWp1B0lqijL6SWLjLiBh/BCEs0ytBt9KQdB8wh9AEU4kk1naigNcGvJAEo9vr35XmGMWBnVUnpCRPHYcRuopWBmTOu3RIbV9KpGTfAfye2CCbu2moo2+x36zYPs32w+mnBxC2CmLTzJ6DZh6YdwGRlr4Z+ACRTFlju+bsdx0dkdUcu4GPA6/aXgp8Jt3fKGlfF8+4AlgmaWm6nm37LCJmnl98Wh2HOvZrjrQ/4kHgaCIuPjt1LSw2OVUDm4G1BV3tzl8dNYzCaGUaUWp+iagsfjJb7+8EA4lcSOGGlZ1E1q2OGkYHh1TSvzhgTuo4zNHbCuRuIn07uOD+YAp2QNVRe+iVcKTy7tPERpMsJhJRSx01jJK/si8DtxH7Np8idmBdQlQCv5HDs+uoInL526eUBPsKkQRrBb4s6fFeP7iOqqJi/wlWR+2j1rbE1dGH+D9jv+5kjCDNVwAAAABJRU5ErkJggg==", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIMAAAAVCAYAAABlol04AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/OQEPoAAAACXBIWXMAABJ0AAASdAHeZh94AAAFWUlEQVR4nO2aXagVVRTHf9c0tc9bBmVppUWRX/fgReqhwkyL6KGvlxK0WxFipJRFIlT/lmnZS4mUiZR5laAUw6iHlFCxjLSr3fwoyyBN0MhKMrUbfpwe9h4d58zMOWfO3CtX/cNhz1mzZs36r71mrzX7nLpiscjpCjO7CDgkaf/J9uVkwczqgSOS/qk7nZPhDE5E1zihmTUDdwH9JB3oWJfygZk1Ai3A45LeiTnfB9gJ3CBpa+RcIn8zGw/MBh6TNK+9/M8DFcSgHtgLDJPU0iVGYRgwBpjRGRLBzJ42s6KZjQ7LJa0HlgIvm9l5MZc2AAeBnyL2yvFv9OP6Wn3PCzXEoAAcBjYDlCQDMB3YB7ydp8PtiLTJeRW4DJgYc64AbJJ0NCIvx78RaAO2VO1p+6GWGPwoqQ0iyWBm1wEjgUWS/s3N1QSYWZPP6OE1mGkE9hN5wgEkrQO2AuPMLJr4DUBrxJ9U/mbWHRgIbJR0uAafo3ZrjUPWGBQIxSDaMzwK1AEfJjg9GrgbGAb0xi0xPwOzJb2XgURmmNkMYHJIdNTMguOxkhb64w+Al4BRwLKQfgF4PWI2lT8wGOgGrDezAvACMBzoAawBJknaXDWZjMghBg3A+8GXaKaMBI4AX8fc+HxgAdAf+AJ4E/gI6AfMM7PJ0WvaGRuAZn/8FWChz6qQ3ho/jgoEZnYucA2RlYEU/h7BcnyVt1sE3vW+jAJW+Kaso1BLDLoBA4hbGXyACsAPCY1TEegj6bew0Myexy1PjwCvVc8nGyQt8oF/GGiWNDdB9Rs/3hqSDfHjxkBQAX+AoX68EbhZ0reh6xfgGs8ngFcqZ5IdNcZgAHA2oWQIrwxXAGcBuxNuvD+aCF6+G9gFXFwZhVwRTM6GJAVJf+MavitD4gZgm6SDIVkqf49gZXgunAgeb/lxcDmnc0bWGBSAXZL+CAThnqGXH/fGGfS7dU/ieobrgQs4MZmiwYlevx23vMZhZajWBWiW1JRmExeIQ8CmMnp/AZcGXyTNAeZEdMrx74ab6B24chlF8KD0SHOkHeKQKQa4vm9tWCGcDEH3XELGzIYAy72xdbiG5E9cA9kPGAt8V8aZmUB9RFYA7sHVve2Rc61pxsysK25yvpf0X5l79+Q4vyQk8vcYhFtWP0l4kwgmeEeZ+8wkpzhkiYGZ9cSViAeAF8MK4WT43Y+9KMVCHIHbJK2KODTVH7akeSJpZlRmZk24IMyP2q0AA3ATl7g8+nt0wfn+Sxl7afzheInYnnD+Pj8uT7tJznHIEoNn/GcxMD+sF06G3cAeXAkIG+qLa7iWxSRCPa5xhI7fkSv4MbU84fjUUWalIYF/CEEylPRGZtYbGAdso0wy5IyCHyuOgaRpwLQ4pWM1X1IRWA1cYmbXhnTa/Njf100AzKwX7n28D65ctFZMIR8ET/C+Mno3+XFlmlIK/wBBo/aQf/MAwG/zLgC6AxPz3IyqALnGILrptARXS+7EbSYhaY+ZrQBGAGvN7HPchtNd3vhRXM1qo2MRrETTzWwQcADYImlxRO8O3N7BxxXYLOEPx2rzENwTeA7QamZLcQlwP3A58JSkzzKzyYZcYxDddFqCq51jI/IHcX1DX2A8bsNmEu59ugtl+oX2gKTVwARcACYA4viyCYCZXQjcC3wqaWcFZpP4D8TV5hbgdlznPg5owv1GMULSrGxMsiPvGJT8n8HMpuAmeWjMu3SngplNAGYBt0j6ssJrThn+UF0M4n61fAP4FZgac67TwL9CTQGWVJoIHqcEf6g+BiXJ4Gv/GKAl3Ch1QlwNzAWereaiU4g/VBmDM397O4Nj+B9XmjyElRxHagAAAABJRU5ErkJggg==", "text/latex": [ - "$\\displaystyle \\int\\limits_{0}^{1} x^{2 k} e^{- \\alpha x^{2}}\\, dx$" + "$\\displaystyle \\left(a + t\\right)^{i} \\left(b + t\\right)^{j}$" ], "text/plain": [ - "1 \n", - "⌠ \n", - "⎮ 2 \n", - "⎮ 2⋅k -α⋅x \n", - "⎮ x ⋅ℯ dx\n", - "⌡ \n", - "0 " + " i j\n", + "(a + t) ⋅(b + t) " ] }, + "execution_count": 28, "metadata": {}, - "output_type": "display_data" - }, + "output_type": "execute_result" + } + ], + "source": [ + "x, a, b = symbols(\"t, a, b\", real=True)\n", + "i, j = symbols(\"i, j\", nonnegative=True, integers=True)\n", + "poly_terms = (x + a)**i * (x + b) **j\n", + "poly_terms" + ] + }, + { + "cell_type": "code", + "execution_count": 29, + "metadata": {}, + "outputs": [ { "data": { - "image/png": "iVBORw0KGgoAAAANSUhEUgAAANMAAAAXCAYAAACRfnp7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAABJ0AAASdAHeZh94AAAIpUlEQVR4nO2be7BXVRXHPzx8UJqSlUyWDQyBMiFXHiGJmCljTgRcErXgKtUUk5KPBAkrv61KYypeOiQq5sUIH4X5ArF8IIKWyKORwcuAgmiZyGMsfBVw+2PtH57OPb/XOb97L1x/35k75/722Xutvc/Za6+1vnufdo2NjbQkzGwJ8JykSc2oYwgwEegHfBz4uqT6PHXnAO8ArwGjgJ7Au8BfgCmS1jVXP6toW2jfCjprgDXNrOMIYB1wGfB2vkpm1g4YDtwLfB74NfA54AvAHuARM/twM/e1ijaCji2pzMy6AMcCa8PvDwK3Aj2AUZK2VEKPpMXA4qCjvkDVAcBhwHJJZ8f6Wge8AZwKPJC1T2Y2DzgH6Crpzdi9K4DpwBhJC7LqqqIpzKwf8CzwLUlzm0NHixoT7pXeBjaYWU/gHtxLnSoprwdpRowEFknak3DvSNxz78qqxMwGAHXAxLghBfQL12ez6soKMzsXOB1/V33w5/A7SWNLaHss8A/cw/8YqAW+BPQGjgP+AzwH3AbcJmlf5UeQDEmrzOxe4Kdmdqek3fnqph1HS4d5NaETI4GngFskjW0lQwIYgYd4SZiFe9CnK6DnWuBfwI157vcN9zdWQFdW/BCYgL+rv5fZdgQ+p+4BRgO3AAOBvwIzgYXAZ4C5wN0hzG5J/BzoAlxapF6qcWT2TGb2M+AHRaqdIWkp/oI+DfwGGC7piQrJLRtm1h3oBjyccG86MBgYLGlvGvkRWT2As4C5SYtGCHV7Ak9KyswGmdk4fMVM+2yuAF4BNuEe6vEy2tYCO4BlQC4fXRT1QGZ2NfAM8BWc8FmYoo+pIOkZM2sAxpvZ1AKeMdU4mhiTmR0CXAKMw1/yDuD3wGTc+DYDj0gaE5rMBOYXGcfWcK3Brf1rQLHEvhy5aTASeDQhf5kBXIBPxhdj92qBu4HVwHmSXooLNTMDrgHGS7oZ+Ab+Qu7K048afBVcFZPTGajHX+QsYJKk/5Y1whSQtN94fCilwcw+hBM3C8IC9Fge+f8MDOq1OOmT2phSzFWAO/HQbSjJC2nqcfyfMQXmagmemD8YlA0DLsdd/j7cCBQRuh3YXsLAP4B7pTp8tfutmQ2RtDpPZ0uSmwEjgHmxPs4CzscNqSGhzTbgzziRcBX+IqPtu4XylXgIAO6V9uJUexL6hut+YzKzgbjxdQZGS/pDyaNqPQwDDgX+WELd3KKQlKuWhDRzNWBFuCYaExnGEfdMd4XOXSbp+tDpX+Ju/xygF1AvaVMJiuI4CWgE1klaaWYnAA+Y2WcllRubF4SZHQF0Dz/bA8ebWQ2wU9JWM/socApwbqTNbNzQRwK7AvMIsDuXrEpaEbzTDqB/gurr8RdxsaR9IYSrAZ7PQzzAe+TDqtCP7wFTgQZgqKQDIY8qBbXAm8CfClUys47AheHnkgz60s7VleE6JI/c1ONoH7l5Fr6KPgnckCsPHmIL7vo6Az8ppKQAaoCNkbzhGnyVuD94rUqiP84SrgE6ARb+z/X9y8BKSa9F2lyMM1ePAq9G/iZGBUt6F1gPnBAtN7PhOONzs6QcK3cc0CHIyYe+wL+B7WZ2HzANWAAMPFgMycwOB74IPCTpnSLVp+LJ+2JJSZ6hFH2p56qkN/BN+uMrPY6oZ6oL15kJiXBO8E2SXi6iJBGS5gBzIr8bgfPSyCpB11I8T8mHJiyepHKYpQZggJl9QtIrZtYJz/G2A1dH6h0Tron0enh5J+Lh42rgY5SxD2JmW4BP5bn9eELOM0/SuFJkl4mh+EZ5wdDIzC4FrsSfX12hukWQda7uxPc748g0jqgxnY7HgPlc71vAdYWUHERYAdyRoX0un+qFhxVTgK7ANyVFDSfnhQ/PI6cP/g4Ow489zS9zQ3EmcHSsrIb38sEtsXtry5BdDkbhey+L8lUwswk4kbIeOFPSzgz6ss7VTiSfjMk0jo6hQifc7b0g6a1Y4254SPNULCw6aCHpFxlF7DcmM3sBJx2exinpKLaF6zEkI0c+XInnb2PNbJWkmaV0IqleoMZH4PnC0lLkZIGZdcDD5sdCCJVU53JgBn7E60xJ25Lqlagv01w1s/b4ArS50uPI5Uyd8LAoiXefga+cqZmXNoioZ5qFL0qXJIQcrwKv47RtEqInHy7ASYhpZjaqst1tVgzBF4vE0MjMJuNzaC3OkqY2pICsc7VnaL82Vp55HDlj2gXsBrqb2UkRAd/B9zmgaTjxfsYmnO4ejZMON0pqcng3GNcy4CNhkziOvni40RDYvmHAy8B8MzuluTpfYdTiE/u++A0z+xGeqK/CV/KCWx1mVm9mjcG75kPWuZp7rvHN6MzjaJf7BMPMbsCPkbyOb3x1CQruB47CN6ZuAm6VtDJJ2PsJZrYRp9+3AT0KhAZfxdm5CZJmR8oPxSfFGkkDI+W98JxuDzCo3G2IrCcgzGwkvj0APgfOBl7EmTOA7ZImhrrtgJeArZIGx+RchG8678UZt6TnsyX6aYyZ3Y4n9HWS8m7YZ5mrZnYHvgh2zREUlRpHlICYhH/Hcz7wbXwFmAF8HzgZuB0YTwse/zjA0YAb0+R8hhSwEDe4C4HZkfLewCHEPkeRtD6EeQ8DD5nZoGIreoVRA1wUK+sW/sAnXW67oD/wSZwIiaNruHbAN1KT8AQ+UXPojW8T5CUAAlLNVTM7Cl8oHowxfRUZR7uW/jiwrcDMlgODgCPjiXBC3Sk4u9Q3KRw8WGFm1+FMZjdJm4vVLyLraHwzfJqkqyrQvSQd38U31k+TtDxSXpFxtMbHgQc9QljQB9hQzJACZuDnCNNueB+oqAX+ltWQAk7D6e7pFZDVBIEFnAIsjBpSQEXGUfVMKRBOgm/AD0OOKVY/tBkCnAH8qsDRoiqaCWZ2Ih4W1qtCH6HG0dIfB7YVnByuJYdskpbhzF4VrQBJz+OnxZsNVc9URRUVQjVnqqKKCuF/hg/HHMc7Mw0AAAAASUVORK5CYII=", "text/latex": [ - "$\\displaystyle \\frac{\\alpha^{- k - \\frac{1}{2}} \\gamma\\left(k + \\frac{1}{2}, \\alpha\\right)}{2}$" + "\\begin{align*}(a + t)^0 (b + t)^0 &= 1\\\\(a + t)^0 (b + t)^1 &= b + t\\\\(a + t)^0 (b + t)^2 &= b^{2} + 2 b t + t^{2}\\\\(a + t)^1 (b + t)^0 &= a + t\\\\(a + t)^1 (b + t)^1 &= a b + t \\left(a + b\\right) + t^{2}\\\\(a + t)^1 (b + t)^2 &= a b^{2} + t \\left(2 a b + b^{2}\\right) + t^{2} \\left(a + 2 b\\right) + t^{3}\\\\(a + t)^2 (b + t)^0 &= a^{2} + 2 a t + t^{2}\\\\(a + t)^2 (b + t)^1 &= a^{2} b + t \\left(a^{2} + 2 a b\\right) + t^{2} \\cdot \\left(2 a + b\\right) + t^{3}\\\\(a + t)^2 (b + t)^2 &= a^{2} b^{2} + t \\left(2 a^{2} b + 2 a b^{2}\\right) + t^{2} \\left(a^{2} + 4 a b + b^{2}\\right) + t^{3} \\cdot \\left(2 a + 2 b\\right) + t^{4}\\\\\\end{align*}" ], "text/plain": [ - " -k - 1/2 \n", - "α ⋅γ(k + 1/2, α)\n", - "───────────────────────\n", - " 2 " + "" ] }, + "execution_count": 29, "metadata": {}, - "output_type": "display_data" + "output_type": "execute_result" } ], "source": [ - "ZZ01 = Integral(x ** (2 * k) * exp(-a * x**2), (x, 0, 1))\n", - "ZZ01_expanded = simplify(ZZ01.doit())\n", - "display(ZZ01, ZZ01_expanded)" + "from itertools import product\n", + "\n", + "\n", + "# output monomial terms with ascending orders x^n, n=0, 1, 2, ...\n", + "output = r\"\\begin{align*}\"\n", + "for ival, jval in product(range(3), range(3)):\n", + " p = simplify(expand(poly_terms.subs(i, ival).subs(j, jval)))\n", + " output += f\"(a + t)^{ival} (b + t)^{jval} &= \"\n", + " output += \" + \".join(\n", + " [latex(p.coeff(x, n) * x**n) for n in range(ival + jval + 1)]\n", + " )\n", + " output += r\"\\\\\"\n", + "\n", + "output += r\"\\end{align*}\"\n", + "IPython.display.Latex(output)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By using the symmetry properties of the binomial expansion we can rewrite the terms in the series of Equation (20) as:\n", + "$$\n", + "\\tag{21}\n", + "\\binom{i}{k} \\binom{j}{l} \\tc^{i-k} \\tc^{j-l} \\CA^k \\CB^l = \\binom{i}{k} \\binom{j}{l} \\tc^k \\tc^l \\CA^{i-k} \\CB^{j-l}.\n", + "$$\n", + "Recognising that we want to collect terms of same power $\\tc^s$ we introduce the constraint that $k+l=s$ allows us to\n", + "write the terms as:\n", + "$$\n", + "\\tag{22}\n", + "\\binom{i}{s-l} \\binom{j}{l} \\tc^s \\CA^{i-(s-l)} \\CB^{j-l}.\n", + "$$\n", + "Replacing the summation variable $l=0, 1,\\cdots, j$ with $t$ and incorporating the constraint we arrive at:\n", + "$$\n", + "\\tag{23}\n", + "(\\tc + \\CA)^i (\\tc + \\CB)^j = \n", + "\\sum_{s = 0}^{i+j} \\tc^s \n", + "\\sum_{\\substack{t=0 \\\\ s - i \\le t \\le j}}^s \n", + "\\binom{i}{s-t} \\binom{j}{t} \\CA^{i-(s-t)} \\CB^{j-t},\n", + "$$\n", + "Finally we can write this series as:\n", + "$$\n", + "\\tag{24}\n", + "(\\tc + \\CA)^i (\\tc + \\CB)^j = \\sum_{s = 0}^{i+j} B(i, j, \\CA, \\CB, s) \\tc^s,\n", + "$$\n", + "where we introduce $B(i, j, \\CA, \\CB, s)$ as the coefficient of $x^s$ in the expansion:\n", + "$$\n", + "\\tag{25}\n", + "B(i, j, \\CA, \\CB, s) = \\sum_{\\substack{t=0 \\\\ s - i \\le t \\le j}}^s \\binom{i}{s-t} \\binom{j}{t} \\CA^{i - (s - t)} \\CB^{j-t}\n", + "$$\n", + "An evaluation strategy for this variable-sized loop is explored in an accompanying [notebook](./binom_factor_table.ipynb).\n", + "\n", + "Combining this result with the one we derived earlier for the product of two\n", + "Gaussian functions Equation 16 gives us:\n", + "$$\n", + "\\tag{26}\n", + "\\tilde{S}_{\\mu \\nu}^{(x)} = \\int_{-\\infty}^\\infty \n", + "\\left((t - X_\\mu)^{l_\\mu} e^{-\\alpha_\\mu (t - X_\\mu)^2} \\right)\n", + "\\left((t - X_\\nu)^{l_\\nu} e^{-\\alpha_\\nu (t - X_\\nu)^2} \\right) dt = \\\\\n", + "\\int_{-\\infty}^\\infty \n", + " e^{-\\alpha_\\mu \\alpha_\\nu (X_\\mu - X_\\nu)^2 / \\gamma} e^{-\\gamma \\tc^2}\n", + "\\sum_{s=0}^{l_\\mu+ l_\\nu} B(l_\\mu, l_\\nu, C_\\mu, C_\\nu, s) \\tc^s dt\n", + "$$\n", + "Swapping the order of integration and the summation we get\n", + "$$\n", + " e^{-\\alpha_\\mu \\alpha_\\nu (X_\\mu - X_\\nu)^2 / \\gamma} \n", + "\\sum_{s=0}^{l_\\mu+ l_\\nu} B(l_\\mu, l_\\nu, C_\\mu, C_\\nu, s) \n", + "\\int_{-\\infty}^\\infty \\tc^s e^{-\\gamma \\tc^2} dt\n", + "$$\n", + "and see that we need to evaluate\n", + "integrals of the form that we already defined above at (Hdef)\n", + "$$\n", + "\\tag{27}\n", + "H(s,\\eta) = \\int_{-\\infty}^\\infty t^s e^{-\\eta t^2} dt\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "And we recall that for $s$ odd, this is zero, so we evaluate it only for $s$ even." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This result can also be written in terms of the [double factorial](https://en.wikipedia.org/wiki/Double_factorial#Additional_identities) which gives us two possible computation strategies for this integral:\n", + "$$\n", + "\\tag{28}\n", + "H(2s, \\eta) = \\int_{-\\infty}^{\\infty} t^{2s} e^{-\\eta t^2} dt \\\\\n", + "= \\eta^{-s - \\frac{1}{2}} \\Gamma\\left(s + \\frac{1}{2}\\right) \\\\\n", + "= \\frac{(2s-1)!!}{(2\\eta)^s} \\sqrt{\\frac{\\pi}{\\eta}}.\n", + "$$\n", + "The last form agrees with Equation (3.15) derived by [Fermann and Valeev](http://arxiv.org/abs/2007.12057).\n", + "\n", + "Using the function $H(2s, \\eta$ allows us to write the one-dimensional overlap integral as:\n", + "$$\n", + "\\tag{29}\n", + "\\tilde{S}_{\\mu \\nu}^{(x)} = \n", + "e^{-\\alpha_\\mu \\alpha_\\nu (X_A - X_B)^2 / \\gamma}\n", + "\\sum_{s=0}^{\\lfloor(l_\\mu + l_\\nu)/2 \\rfloor} B(l_\\mu, l_\\nu, CA_x, CB_x, 2s)\\;H(2s, \\gamma)\n", + "$$\n", + "substituting this back into Equation (11) gives us the overlap of two primitive Gaussians:\n", + "$$\n", + "\\tag{30}\n", + "\\tilde{S}_{\\mu \\nu} = \\iiint p_\\mu(\\br) p_\\nu(\\br) dx dy dz \\\\\n", + "= N_\\mu N_\\nu \\exp\\left(-\\frac{\\alpha_\\mu \\alpha_\\nu |\\mathbf{A}-\\mathbf{B}|^2}{\\alpha_\\mu + \\alpha_\\nu}\\right) \\times\\\\\n", + "\\sum_{s=0}^{\\lfloor(l_\\mu + l_\\nu)/2 \\rfloor} B(l_\\mu, l_\\nu, CA_x, CB_x, 2s)\\;G(\\gamma, 2s) \\times \\\\\n", + "\\sum_{s=0}^{\\lfloor(m_\\mu + m_\\nu)/2 \\rfloor} B(m_\\mu, m_\\nu, CA_y, CB_y, 2s)\\;G(\\gamma, 2s) \\times \\\\\n", + "\\sum_{s=0}^{\\lfloor(n_\\mu + n_\\nu)/2 \\rfloor} B(n_\\mu, n_\\nu, CA_z, CB_z, 2s)\\;G(\\gamma, 2s)\n", + "$$" ] } ],