{ "cells": [ { "cell_type": "markdown", "id": "8890a039", "metadata": {}, "source": [ "# Mcmurchie Davidson Scheme\n" ] }, { "cell_type": "markdown", "id": "f43933a5", "metadata": {}, "source": [ "## Overlap\n", "\n", "Let us frist consider the simple overlap integral\n", "\\begin{equation}\n", "S_{ab} = \\langle G_a|G_b \\rangle\n", "\\end{equation}\n", "\n", "it can be factorized in the three Cartesian directons\n", "\\begin{equation}\n", "S_{ab} = S_{ij}S_{kl}S_{mn}\n", "\\end{equation}\n", "\n", "The Obara-Saika recurrence relations for the Cartesian overlap integral over one direction is:\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "S_{i+1,j} &= X_{PA}S_{ij} + \\frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})\\\\\n", "S_{i,j+1} &= X_{PB}S_{ij} + \\frac{1}{2p}(iS_{i-1,j}+jS_{i,j-1})\\\\\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "With the boundary condtion\n", "\n", "\\begin{equation}\n", "S_{0,0} = \\sqrt(\\frac{\\pi}{p}) exp(-\\mu X_{AB}^2)\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 1, "id": "4a449dfb", "metadata": {}, "outputs": [], "source": [ "class Overlap(object):\n", " \"\"\"The Obara-Saika scheme for three-dimensional overlap integral over\n", " primitive Gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " a : float \n", " Gaussian exponent facotr.\n", "\n", " b : float \n", " Gaussian exponent facotr.\n", "\n", " i : int\n", " Angular momentum quantum number.\n", "\n", " j : int\n", " Angular momentum quantum number.\n", "\n", " A : float\n", " Coordinate in on direction.\n", "\n", " B : float\n", " Coordinate in on direction.\n", "\n", " Returns\n", " -------\n", " result : float\n", " The non-normalizecd overlap interal.\n", " \"\"\"\n", " def __init__(self):\n", " self.p = 0\n", " self.mu = 0\n", " self.P = ()\n", " \n", " def __call__(self, pga, pgb):\n", " \"\"\"Evaluates non-normalizecd overlap integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " pga: PrimitiveGaussian\n", " The first primitive gaussian orbital.\n", "\n", " pgb: PrimitiveGaussian\n", " The second primitive gaussian orbital.\n", " \n", " C: List[float,float,float]\n", " Coordinate of nuclei.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " result = 1\n", " for r in range(3):\n", " result *= self.S1d(r, pga, pgb)\n", " return result\n", " \n", " def S1d(self, r, pga, pgb):\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a + b\n", " mu = (a*b)/(a+b)\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " XAB = A-B\n", "\n", " if pga.shell[r] > 0:\n", " return self.recursive(r, *self.gaussian_factory(r, pga, pgb))\n", " elif pgb.shell[r] > 0:\n", " return self.recursive(r, *self.gaussian_factory(r, pgb, pga))\n", " else:\n", " # Starting from the spherical Gaussians.\n", " S00 = np.power(np.pi/p,0.5)*np.exp(-mu*XAB[r]**2)\n", " return S00\n", " \n", " def recursive(self, r, pga, pgb, pga_1, pga_2, pgb_1):\n", " term1 = term2 = term3 = 0\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a + b\n", " mu = (a*b)/(a+b)\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " P = (a*A+b*B)/p\n", " XPA = P-A\n", "\n", " if XPA[r] != 0:\n", " term1 = XPA[r] * self.S1d(r, pga_1, pgb)\n", " if pga_1.shell[r] >= 0:\n", " term2 = pga_1.shell[r] * (1 / (2 * p)) * self.S1d(r, pga_2, pgb)\n", " if pgb.shell[r] >= 0:\n", " term3 = pgb.shell[r] * (1 / (2 * p)) * self.S1d(r, pga_1, pgb_1)\n", " return term1 + term2 + term3\n", " \n", " def gaussian_factory(self, r, pga, pgb):\n", " ca = pga.coefficient\n", " cb = pgb.coefficient\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", "\n", " A = pga.origin\n", " B = pgb.origin\n", "\n", " i,k,m = pga.shell\n", " j,l,n = pgb.shell\n", "\n", " if r == 0:\n", " pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)\n", " pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)\n", " pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)\n", " return pga, pgb, pga_i_1, pga_i_2, pgb_j_1\n", " elif r == 1:\n", " pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)\n", " pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)\n", " pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)\n", " return pga, pgb, pga_k_1, pga_k_2, pgb_l_1\n", " elif r == 2:\n", " pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)\n", " pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)\n", " pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)\n", " return pga, pgb, pga_m_1, pga_m_2, pgb_n_1" ] }, { "cell_type": "code", "execution_count": 4, "id": "f75f87b3", "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-8.880194157791243e-05\n" ] } ], "source": [ "# The overlap integral between two primitive Gaussian type orbital\n", "S_integral = Overlap()\n", "print(S_integral(pga,pgb))" ] }, { "cell_type": "markdown", "id": "9692f9d3", "metadata": {}, "source": [ "## Kinetic\n", "\n", "Consider the integral of electrons kinetic energy\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "T_{a,b} &= \\langle G_a \\vert -\\sum^N_{i=1}\\frac{\\hbar^2}{2m_i}\\boldsymbol{\\nabla}_i^2 \\vert G_b \\rangle\\\\\n", "T_{a,b} &= T_{ij}S_{kl}S_{mn} + S_{ij}T_{kl}S_{mn} + S_{ij}S_{kl}T_{mn}\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "\n", "The Obara-Saika recurrence relations for the Cartesian kinetic integrals over one direction is:\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "T_{i+1,j} &= X_{PA}T_{ij} + \\frac{1}{2p}(iT_{i-1,j}+jT_{i,j-1}) + \\frac{b}{p}(2aS_{i+1,j}-iS_{i-1,j})\\\\\n", "T_{i,j+1} &= X_{PB}T_{ij} + \\frac{1}{2p}(iT_{i-1,j}+jT_{i,j-1})+ \\frac{a}{p}(2bS_{i,j+1}-jS_{i,j-1})\\\\\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "With the boundary condtion\n", "\n", "\\begin{equation}\n", "T_{0,0} = \\left[a-2a^2(X_{PA}^2 + \\frac{1}{2p})\\right] S_{00}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 5, "id": "d5f71a86", "metadata": {}, "outputs": [], "source": [ "overlap = Overlap()\n", "class Kinetic(object):\n", " \"\"\"The Obara-Saika scheme for three-dimensional kinetic energy integral over\n", " primitive Gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " a : float \n", " Gaussian exponent facotr.\n", "\n", " b : float \n", " Gaussian exponent facotr.\n", "\n", " i : int\n", " Angular momentum quantum number.\n", "\n", " j : int\n", " Angular momentum quantum number.\n", "\n", " A : float\n", " Coordinate in on direction.\n", "\n", " B : float\n", " Coordinate in on direction.\n", "\n", " Returns\n", " -------\n", " result : float\n", " The non-normalizecd kinetic interals in one dimension.\n", " \"\"\"\n", " def __init__(self):\n", " self.p = 0\n", " self.mu = 0\n", " self.P = ()\n", " \n", " def __call__(self, pga, pgb):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " pga: PrimitiveGaussian\n", " The first primitive gaussian orbital.\n", "\n", " pgb: PrimitiveGaussian\n", " The second primitive gaussian orbital.\n", " \n", " C: List[float,float,float]\n", " Coordinate of nuclei.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " Sij = overlap.S1d(0,pga,pgb)\n", " Skl = overlap.S1d(1,pga,pgb)\n", " Smn = overlap.S1d(2,pga,pgb)\n", "\n", " Tij = self.T1d(0,pga,pgb)\n", " Tkl = self.T1d(1,pga,pgb)\n", " Tmn = self.T1d(2,pga,pgb)\n", "\n", " Tab = Tij*Skl*Smn+Sij*Tkl*Smn+Sij*Skl*Tmn\n", " return Tab\n", "\n", " def T1d(self, r, pga, pgb):\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a + b\n", " mu = (a*b)/(a+b)\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " P = (a*A+b*B)/p\n", " XAB = A-B\n", " XPA = P-A\n", "\n", " if pga.shell[r] > 0:\n", " return self.recursive(r, *self.gaussian_factory(r, pga, pgb))\n", " elif pgb.shell[r] > 0:\n", " return self.recursive(r, *self.gaussian_factory(r, pgb, pga))\n", " else:\n", " # Starting from the spherical Gaussians.\n", " S00 = np.power(np.pi/p,0.5)*np.exp(-mu*XAB[r]**2)\n", " T00 = (a-2*a**2*(XPA[r]**2+1./(2*p)))*S00\n", " return T00\n", " \n", " def recursive(self, r, pga, pgb, pga_1, pga_2, pgb_1):\n", " term1 = term2 = term3 = term4 = term5 = 0\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a + b\n", " mu = (a*b)/(a+b)\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " P = (a*A+b*B)/p\n", " XPA = P-A\n", "\n", " if XPA[r] != 0:\n", " term1 = XPA[r] * self.T1d(r, pga_1, pgb)\n", " if pga_1.shell[r] >= 0:\n", " term2 = pga_1.shell[r] * (1 / (2 * p)) * self.T1d(r, pga_2, pgb)\n", " if pgb.shell[r] >= 0:\n", " term3 = pgb.shell[r] * (1 / (2 * p)) * self.T1d(r, pga_1, pgb_1)\n", " term4 = (2*a*b) / p * overlap.S1d(r, pga, pgb)\n", " if pga_1.shell[r] >= 0:\n", " term5 = pgb.shell[r] * (b / p) * overlap.S1d(r, pga_2, pgb)\n", " return term1 + term2 + term3 + term4 - term5\n", "\n", " def gaussian_factory(self, r, pga, pgb):\n", " ca = pga.coefficient\n", " cb = pgb.coefficient\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", "\n", " A = pga.origin\n", " B = pgb.origin\n", "\n", " i,k,m = pga.shell\n", " j,l,n = pgb.shell\n", "\n", " if r == 0:\n", " pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)\n", " pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)\n", " pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)\n", " return pga, pgb, pga_i_1, pga_i_2, pgb_j_1\n", " elif r == 1:\n", " pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)\n", " pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)\n", " pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)\n", " return pga, pgb, pga_k_1, pga_k_2, pgb_l_1\n", " elif r == 2:\n", " pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)\n", " pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)\n", " pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)\n", " return pga, pgb, pga_m_1, pga_m_2, pgb_n_1" ] }, { "cell_type": "code", "execution_count": 6, "id": "22845b36", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0016734344532709933\n" ] } ], "source": [ "# The kinetic integral between two primitive Gaussian type orbital\n", "T_integral = Kinetic()\n", "print(T_integral(pga,pgb))" ] }, { "cell_type": "markdown", "id": "cb7a7840", "metadata": {}, "source": [ "## Nuclear Attraction\n", "\n", "\n", "\\begin{equation} \n", "\\begin{aligned}\n", "V_{a,b} &=\\langle G_a \\vert -\\sum^N_{i=1}\\sum^M_{\\alpha=1} \\frac{Z_\\alpha e^2}{\\textbf{r}_{i\\alpha}} \\vert G_b \\rangle \\\\\n", "V_{a,b} &= V_{ijklmn}^{000} = \\Theta_{ijklmn}^{0}\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "\n", "The Obara-Saika recurrence relations for the Cartesian nuclear attraction integrals:\n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\Theta_{i+1,j,k,l,m,n}^{N} &= X_{PA}\\Theta_{ijkln}^{N}\n", "+\\frac{1}{2p}(i\\Theta_{i-1,j,k,l,m,n}^{N} + j\\Theta_{i,j-1,k,l,m,n}^{N})\\\\\n", "&-X_{PC}\\Theta_{ijklmn}^{N+1}\n", "-\\frac{1}{2p}(i\\Theta_{i-1,j,k,l,m,n}^{N+1} + j\\Theta_{i,j-1,k,l,m,n}^{N+1})\\\\\n", "\\Theta_{i,j+1,k,l,m,n}^{N} &= X_{PB}\\Theta_{ijkln}^{N}\n", "+\\frac{1}{2p}(i\\Theta_{i-1,j,k,l,m,n}^{N} + j\\Theta_{i,j-1,k,l,m,n}^{N})\\\\\n", "&-X_{PC}\\Theta_{ijklmn}^{N+1}\n", "-\\frac{1}{2p}(i\\Theta_{i-1,j,k,l,m,n}^{N+1} + j\\Theta_{i,j-1,k,l,m,n}^{N+1})\\\\\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "With the boundary condtion\n", "\n", "\\begin{equation}\n", "\\Theta_{000000}^{N} = \\frac{2\\pi}{p}(-2p)^{-N}K_{ab}^{xyz}R_{000}^{N} \n", "= \\frac{2\\pi}{p} K_{ab}^{xyz}F_N(pR_{PC}^2)\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 7, "id": "44395d35", "metadata": {}, "outputs": [], "source": [ "class NuclearAttraction(object):\n", " \"\"\"The Obara-Saika scheme for three-dimensional nuclear attraction integral over\n", " primitive Gaussian orbitals.\n", "\n", " Attributes\n", " ----------\n", " p : float \n", " The total exponent.\n", "\n", " mu : float\n", " The reduced exponent.\n", "\n", " P : List[float,float,float] \n", " The centre of charge coordinate.\n", "\n", " C : List[float,float,float] \n", " The coordinate of given nuclei.\n", "\n", " Kab : float\n", " The pre-exponential factor.\n", "\n", " Methods\n", " -------\n", " __init__(self)\n", " Initialize the instance.\n", " \"\"\"\n", " def __init__(self):\n", " \"\"\"Initialize the instance.\n", " \"\"\"\n", " self.p = 0\n", " self.mu = 0\n", " self.Kab = 0\n", " self.P = []\n", " self.C = []\n", " self.boys_dict = {}\n", "\n", " def __call__(self, pga, pgb, C):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " pga: PrimitiveGaussian\n", " The first primitive gaussian orbital.\n", "\n", " pgb: PrimitiveGaussian\n", " The second primitive gaussian orbital.\n", " \n", " C: List[float,float,float]\n", " Coordinate of nuclei.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " l_total = sum(pga.shell) + sum(pgb.shell)\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a + b\n", " mu = (a*b)/(a+b)\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " P = (a*A+b*B)/p\n", "\n", " RAB = np.linalg.norm(A-B)\n", " RPA = np.linalg.norm(P-A)\n", " RPB = np.linalg.norm(P-B)\n", " RPC = np.linalg.norm(P-C)\n", "\n", " Kab = exp(-mu*RAB**2)\n", " \n", " self.p = p\n", " self.mu = mu\n", " self.P = P\n", " self.C = C\n", " self.Kab = Kab\n", "\n", " # Build boys function F_{N}(x)\n", " N = l_total\n", " x = p*RPC**2\n", " boys_pre_factor = (2*np.pi)/p*Kab\n", " boys_function = boys(l_total, x)\n", " Theta_N_000000 = boys_pre_factor * boys_function\n", " self.boys_dict = {l_total: Theta_N_000000}\n", "\n", " while N >= 1:\n", " boys_function = boys_recursion(N, x, boys_function)\n", " N -= 1\n", " Theta_N_000000 = boys_pre_factor * boys_function\n", " self.boys_dict[N] = Theta_N_000000\n", "\n", " result = self.V(0, pga, pgb)\n", " return result\n", "\n", " def V(self, N, pga, pgb):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga : PrimitiveGaussian\n", " The first primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The second primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " vlaue : float\n", " Integral value.\n", " \"\"\"\n", " if pga.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pga, pgb))\n", " elif pga.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pga, pgb))\n", " elif pga.shell[2] > 0:\n", " return self.recursive(2, N, *self.gaussian_factory(2, pga, pgb))\n", " elif pgb.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pgb, pga))\n", " elif pgb.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pgb, pga))\n", " elif pgb.shell[2] > 0:\n", " result = self.recursive(2, N, *self.gaussian_factory(2, pgb, pga))\n", " return result\n", " else:\n", " return self.boys_dict[N]\n", "\n", " def recursive(self, r, N, pga, pgb, pga_1, pga_2, pgb_1):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " r : int\n", " Cartesian index 0, 1, 2. \n", "\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga_1 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pga_2 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb_1 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " term1 = term2 = term3 = term4 = term5 = term6 = 0\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " p = a+b\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " P = (a*A+b*B)/p\n", " C = self.C\n", "\n", " XPA = np.array(P) - np.array(A)\n", " XPC = np.array(P) - np.array(C)\n", "\n", " if np.array_equal(P,A) is False:\n", " term1 = XPA[r] * self.V(N, pga_1, pgb)\n", " if pga_1.shell[r] > 0:\n", " term2 = pga_1.shell[r] * (1 / (2 * p)) * self.V(N, pga_2, pgb)\n", " if pgb.shell[r] > 0:\n", " term3 = pgb.shell[r] * (1 / (2 * p)) * self.V(N, pga_1, pgb_1)\n", " if np.array_equal(P,C) is False:\n", " term4 = XPC[r] * self.V(N+1, pga_1, pgb)\n", " if pga_1.shell[r] > 0:\n", " term5 = pga_1.shell[r] * (1 / (2 * p)) * self.V(N+1, pga_2, pgb)\n", " if pgb.shell[r] > 0:\n", " term6 = pgb.shell[r] * (1 / (2 * p)) * self.V(N+1, pga_1, pgb_1)\n", "\n", " result = term1+term2+term3-term4-term5-term6\n", " return result\n", "\n", " def gaussian_factory(self, r, pga, pgb):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " r : int\n", " Cartesian index 0, 1, 2. \n", "\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " result : Tuple(pg,pg,pg,pg)\n", " Tuple of 4 PrimitiveGaussian orbital instance. \n", " \"\"\"\n", " ca = pga.coefficient\n", " cb = pgb.coefficient\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", "\n", " A = pga.origin\n", " B = pgb.origin\n", "\n", " i,k,m = pga.shell\n", " j,l,n = pgb.shell\n", "\n", " if r == 0:\n", " pga_i_1 = PrimitiveGaussian(ca, A, (i - 1, k, m), a)\n", " pga_i_2 = PrimitiveGaussian(ca, A, (i - 2, k, m), a)\n", " pgb_j_1 = PrimitiveGaussian(cb, B, (j - 1, l, n), b)\n", " return pga, pgb, pga_i_1, pga_i_2, pgb_j_1\n", " elif r == 1:\n", " pga_k_1 = PrimitiveGaussian(ca, A, (i, k - 1, m), a)\n", " pga_k_2 = PrimitiveGaussian(ca, A, (i, k - 2, m), a)\n", " pgb_l_1 = PrimitiveGaussian(cb, B, (j, l - 1, n), b)\n", " return pga, pgb, pga_k_1, pga_k_2, pgb_l_1\n", " elif r == 2:\n", " pga_m_1 = PrimitiveGaussian(ca, A, (i, k, m - 1), a)\n", " pga_m_2 = PrimitiveGaussian(ca, A, (i, k, m - 2), a)\n", " pgb_n_1 = PrimitiveGaussian(cb, B, (j, l, n - 1), b)\n", " return pga, pgb, pga_m_1, pga_m_2, pgb_n_1 " ] }, { "cell_type": "code", "execution_count": 8, "id": "80cc7338", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "-8.543259485351415e-05\n" ] } ], "source": [ "# The nuclear attraction integral between two primitive Gaussian type orbital\n", "V_integral = NuclearAttraction()\n", "print(V_integral(pga,pgb,FCenter[0]))" ] }, { "cell_type": "markdown", "id": "be27bf0f", "metadata": {}, "source": [ "## Electron Repulsion\n", "\n", "\n", "\\begin{equation} \n", "V_{abcd} = \\langle G_a G_b\\vert \\sum^N_{i=1}\\sum^N_{j>i} \\frac{e^2}{\\textbf{r}_{ij}}\\vert G_c G_d \\rangle \n", "\\end{equation}\n", "\n", "The source and target integrals:\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\Theta_{0000;0000;0000}^{N} &= \\frac{2\\pi^{2.5}}{pq\\sqrt{p+q}}\n", "K_{ab}^{xyz}K_{cd}^{xyz}F_N(\\alpha R_{PQ}^2)\\\\ \n", "\\Theta_{i_xj_xk_xl_x;i_yj_yk_yl_y;i_zj_zk_zl_z}^{0} &= g_{i_xj_xk_xl_x;i_yj_yk_yl_y;i_zj_zk_zl_z}\n", "\\end{aligned}\n", "\\end{equation}\n", "\n", "The Obara-Saika two electron recurrence relation\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\Theta_{i+1,j,k,l}^{N} &= \n", "X_{PA}\\Theta_{ijkl}^{N}\n", "-\\frac{\\alpha}{p}X_{PQ}\\Theta_{i,j,k,l}^{N+1} \n", "+\\frac{i}{2p}\\left(\\Theta_{i-1,j,k,l}^{N}-\\frac{\\alpha}{p}\\Theta_{i-1,j,k,l}^{N+1}\\right)\\\\\n", "&+\\frac{j}{2p}\\left(\\Theta_{i,j-1,k,l}^{N}-\\frac{\\alpha}{p}\\Theta_{i,j-1,k,l}^{N+1}\\right)\n", "-\\frac{k}{2(p+q)}\\Theta_{i,j,k-1,l}^{N+1} \n", "-\\frac{l}{2(p+q)}\\Theta_{i,j,k,l-1}^{N+1} \n", "\\end{aligned}\n", "\\end{equation}\n", "Using the horizontal recurrence relation, a similar relation may be written down for increments in j, replaceing $X_{PA}$ with $X_{PB}$. \n", "\n", "\\begin{equation}\n", "\\begin{aligned}\n", "\\Theta_{i,j,k+1,l}^{N} &= \n", "X_{QC}\\Theta_{ijkl}^{N}\n", "-\\frac{\\alpha}{q}X_{PQ}\\Theta_{i,j,k,l}^{N+1} \n", "+\\frac{k}{2q}\\left(\\Theta_{i,j,k-1,l}^{N}-\\frac{\\alpha}{q}\\Theta_{i,j,k-1,l}^{N+1}\\right)\\\\\n", "&+\\frac{l}{2q}\\left(\\Theta_{i,j,k,l-1}^{N}-\\frac{\\alpha}{q}\\Theta_{i,j,k,l-1}^{N+1}\\right)\n", "-\\frac{i}{2(p+q)}\\Theta_{i-1,j,k,l}^{N+1} \n", "-\\frac{j}{2(p+q)}\\Theta_{i,j-1,k,l}^{N+1} \n", "\\end{aligned}\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": 9, "id": "5e59d476", "metadata": {}, "outputs": [], "source": [ "class ElectronRepulsion:\n", " \"\"\"The Obara-Saika scheme for three-dimensional nuclear attraction integral over\n", " primitive Gaussian orbitals.\n", "\n", " Attributes\n", " ----------\n", " p : float \n", " The total exponent.\n", "\n", " mu : float\n", " The reduced exponent.\n", "\n", " P : List[float,float,float] \n", " The centre of charge coordinate.\n", "\n", " C : List[float,float,float] \n", " The coordinate of given nuclei.\n", "\n", " Kab : float\n", " The pre-exponential factor.\n", "\n", " Methods\n", " -------\n", " __init__(self)\n", " Initialize the instance.\n", " \"\"\"\n", " def __init__(self):\n", " \"\"\"Initialize the instance.\n", " \"\"\"\n", " self.alpha = 0\n", " self.R = []\n", " self.boys_dict = {}\n", " \n", " def __call__(self, pga, pgb, pgc, pgd):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " pga: PrimitiveGaussian\n", " The first primitive gaussian orbital.\n", "\n", " pgb: PrimitiveGaussian\n", " The second primitive gaussian orbital.\n", " \n", " C: List[float,float,float]\n", " Coordinate of nuclei.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " l_total = sum(pga.shell) + sum(pgb.shell) + sum(pgc.shell) + sum(pgd.shell)\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " c = pgc.exponent\n", " d = pgd.exponent\n", " p = a + b\n", " q = c + d\n", " mu = (a*b)/(a+b)\n", " nu = (c*d)/(c+d)\n", " alpha = (p * q) / (p + q)\n", " self.alpha = alpha\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " C = np.array(pgc.origin)\n", " D = np.array(pgd.origin)\n", " P = (a*A+b*B)/(a+b)\n", " Q = (c*C+d*D)/(c+d)\n", " R = (p*P+q*Q)/(p+q)\n", " self.R = R\n", "\n", " RAB = np.linalg.norm(A-B)\n", " RCD = np.linalg.norm(C-D)\n", " RPQ = np.linalg.norm(P-Q)\n", "\n", " # Build boys function F_{N}(x)\n", " Kab = np.exp(-mu*RAB**2)\n", " Kcd = np.exp(-nu*RCD**2)\n", " boys_pre_factor = (2*np.pi**(5/2))/(p*q*np.sqrt(p+q))*Kab*Kcd\n", " N = l_total\n", " x = alpha*RPQ**2\n", " boys_function = boys(l_total, x)\n", " Theta_N_0000_0000_0000 = boys_pre_factor * boys_function\n", " self.boys_dict = {l_total: Theta_N_0000_0000_0000}\n", "\n", " while N >= 1:\n", " boys_function = boys_recursion(N, x, boys_function)\n", " N -= 1\n", " Theta_N_0000_0000_0000 = boys_pre_factor * boys_function\n", " self.boys_dict[N] = Theta_N_0000_0000_0000\n", " result = self.Eri(0, pga, pgb, pgc, pgd)\n", " return result\n", " \n", " def Eri(self, N, pga, pgb, pgc, pgd):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgc : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgd : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " vlaue : float\n", " Integral value.\n", " \"\"\"\n", " if pga.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pga, pgb, pgc, pgd))\n", " elif pga.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pga, pgb, pgc, pgd))\n", " elif pga.shell[2] > 0:\n", " return self.recursive(2, N, *self.gaussian_factory(2, pga, pgb, pgc, pgd))\n", " elif pgb.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pgb, pga, pgd, pgc))\n", " elif pgb.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pgb, pga, pgd, pgc))\n", " elif pgb.shell[2] > 0:\n", " return self.recursive(2, N, *self.gaussian_factory(2, pgb, pga, pgd, pgc))\n", " elif pgc.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pgc, pgd, pga, pgb))\n", " elif pgc.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pgc, pgd, pga, pgb))\n", " elif pgc.shell[2] > 0:\n", " return self.recursive(2, N, *self.gaussian_factory(2, pgc, pgd, pga, pgb))\n", " elif pgd.shell[0] > 0:\n", " return self.recursive(0, N, *self.gaussian_factory(0, pgd, pgc, pgb, pga))\n", " elif pgd.shell[1] > 0:\n", " return self.recursive(1, N, *self.gaussian_factory(1, pgd, pgc, pgb, pga))\n", " elif pgd.shell[2] > 0:\n", " return self.recursive(2, N, *self.gaussian_factory(2, pgd, pgc, pgb, pga))\n", " else:\n", " return self.boys_dict[N]\n", "\n", " def recursive(self, r, N, pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " r : int\n", " Cartesian index 0, 1, 2. \n", "\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga_1 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pga_2 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb_1 : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " result : float\n", " Integral value.\n", " \"\"\"\n", " term1 = term2 = term3 = term4 = term5 = term6 = term7 = term8 = 0\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " c = pgc.exponent\n", " d = pgd.exponent\n", " p = a + b\n", " q = c + d\n", " alpha = (p*q)/(p+q)\n", " #alpha = self.alpha\n", "\n", " A = np.array(pga.origin)\n", " B = np.array(pgb.origin)\n", " C = np.array(pgc.origin)\n", " D = np.array(pgd.origin)\n", " P = (a*A+b*B)/(a+b)\n", " Q = (c*C+d*D)/(c+d)\n", "\n", " XPA = P - A\n", " XPQ = P - Q\n", "\n", " if XPA[r] != 0:\n", " term1 = XPA[r] * self.Eri(N, pga_1, pgb, pgc, pgd)\n", " if XPQ[r] != 0:\n", " term2 = alpha/p*XPQ[r] * self.Eri(N+1, pga_1, pgb, pgc, pgd)\n", " if pga_1.shell[r] > 0:\n", " term3 = pga_1.shell[r] * (1 / (2 * p)) * self.Eri(N, pga_2, pgb, pgc, pgd)\n", " term4 = pga_1.shell[r] * (alpha / (2 * p ** 2)) * self.Eri(N+1, pga_2, pgb, pgc, pgd)\n", " if pgb.shell[r] > 0:\n", " term5 = pgb.shell[r] * (1 / (2 * p)) * self.Eri(N, pga_1, pgb_1, pgc, pgd)\n", " term6 = pgb.shell[r] * (alpha / (2 * p ** 2)) * self.Eri(N+1, pga_1, pgb_1, pgc, pgd)\n", " if pgc.shell[r] > 0:\n", " term7 = pgc.shell[r] * (1 / (2 * (p + q))) * self.Eri(N+1, pga_1, pgb, pgc_1, pgd)\n", " if pgd.shell[r] > 0:\n", " term8 = pgd.shell[r] * (1 / (2 * (p + q))) * self.Eri(N+1, pga_1, pgb, pgc, pgd_1)\n", "\n", " return term1 - term2 + term3 - term4 + term5 - term6 + term7 + term8\n", "\n", " def gaussian_factory(self, r, pga, pgb, pgc, pgd):\n", " \"\"\"Evaluates nuclear attraction integral over two primitive gaussian orbitals.\n", "\n", " Parameters\n", " ----------\n", " r : int\n", " Cartesian index 0, 1, 2. \n", "\n", " N : int\n", " Order of the boys function F_{N}(x).\n", "\n", " pga : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " pgb : PrimitiveGaussian\n", " The primitive gaussian orbital.\n", "\n", " Return\n", " ------\n", " result : Tuple(pg,pg,pg,pg)\n", " Tuple of 4 PrimitiveGaussian orbital instance. \n", " \"\"\"\n", " ca = pga.coefficient\n", " cb = pgb.coefficient\n", " cc = pgc.coefficient\n", " cd = pgd.coefficient\n", "\n", " a = pga.exponent\n", " b = pgb.exponent\n", " c = pgc.exponent\n", " d = pgd.exponent\n", "\n", " A = pga.origin\n", " B = pgb.origin\n", " C = pgc.origin\n", " D = pgd.origin\n", "\n", " ix,iy,iz = pga.shell\n", " jx,jy,jz = pgb.shell\n", " kx,ky,kz = pgc.shell\n", " lx,ly,lz = pgd.shell\n", "\n", " if r == 0:\n", " pga_1 = PrimitiveGaussian(ca, A, (ix - 1, iy, iz), a)\n", " pga_2 = PrimitiveGaussian(ca, A, (ix - 2, iy, iz), a)\n", " pgb_1 = PrimitiveGaussian(cb, B, (jx - 1, jy, jz), b)\n", " pgc_1 = PrimitiveGaussian(cc, C, (kx - 1, ky, kz), c)\n", " pgd_1 = PrimitiveGaussian(cd, D, (lx - 1, ly, lz), d)\n", " return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1\n", " elif r == 1:\n", " pga_1 = PrimitiveGaussian(ca, A, (ix, iy-1, iz), a)\n", " pga_2 = PrimitiveGaussian(ca, A, (ix, iy-2, iz), a)\n", " pgb_1 = PrimitiveGaussian(cb, B, (jx, jy-1, jz), b)\n", " pgc_1 = PrimitiveGaussian(cc, C, (kx, ky-1, kz), c)\n", " pgd_1 = PrimitiveGaussian(cd, D, (lx, ly-1, lz), d)\n", " return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1\n", " elif r == 2:\n", " pga_1 = PrimitiveGaussian(ca, A, (ix, iy, iz-1), a)\n", " pga_2 = PrimitiveGaussian(ca, A, (ix, iy, iz-2), a)\n", " pgb_1 = PrimitiveGaussian(cb, B, (jx, jy, jz-1), b)\n", " pgc_1 = PrimitiveGaussian(cc, C, (kx, ky, kz-1), c)\n", " pgd_1 = PrimitiveGaussian(cd, D, (lx, ly, lz-1), d)\n", " return pga, pgb, pgc, pgd, pga_1, pga_2, pgb_1, pgc_1, pgd_1" ] }, { "cell_type": "code", "execution_count": 10, "id": "481695d8", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.9060888184873294e-08\n" ] } ], "source": [ "# The overlap integral between two primitive Gaussian type orbital\n", "Eri_integral = ElectronRepulsion()\n", "print(Eri_integral(pga,pgb,pga,pgb))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.6" }, "vscode": { "interpreter": { "hash": "3b4dd1ce6c8ed31661adb8e55d1f9c30ec54ed472785ff10576acd13e5f2c743" } } }, "nbformat": 4, "nbformat_minor": 5 }