Docstring equations in evals and integrals modules

docs/notes.tex

@@ -855,19 +855,19 @@ \subsubsection{Laplacian of density}
 \end{itemize}
 \subsubsection{Hessian of density}
 \begin{equation}
-  H[\rho(\mathbf{r}_n)]
-  =
-  \begin{bmatrix}
-    \frac{\partial^2}{\partial x^2} \rho(\mathbf{r}_n) &
-    \frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}_n) &
-    \frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n)\\\\
-    \frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}_n) &
-    \frac{\partial^2}{\partial y^2} \rho(\mathbf{r}_n)&
-    \frac{\partial^2}{\partial y \partial z} \rho(\mathbf{r}_n)\\\\
-    \frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n) &
-    \frac{\partial^2}{\partial z^2} \rho(\mathbf{r}_n)&
-    \frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n)\\
-  \end{bmatrix}\\
+	H[\rho(\mathbf{r}_n)]
+	=
+	\begin{bmatrix}
+		\frac{\partial^2}{\partial x^2} \rho(\mathbf{r}_n) &
+		\frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}_n) &
+		\frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r}_n)\\\\
+		\frac{\partial^2}{\partial y \partial x} \rho(\mathbf{r}_n) &
+		\frac{\partial^2}{\partial y^2} \rho(\mathbf{r}_n)&
+		\frac{\partial^2}{\partial y \partial z} \rho(\mathbf{r}_n)\\\\
+		\frac{\partial^2}{\partial z \partial x} \rho(\mathbf{r}_n) &
+		\frac{\partial^2}{\partial z \partial y} \rho(\mathbf{r}_n)&
+		\frac{\partial^2}{\partial z^2} \rho(\mathbf{r}_n)\\
+	\end{bmatrix}\\
 \end{equation}
 
 Example:

gbasis/contractions.py

@@ -17,27 +17,27 @@ class GeneralizedContractionShell:
 
     .. math::
 
-        g_i (\mathbf{r} | \vec{a}, \vec{R}_A) =
-        N_g (\alpha, \vec{a})
+        g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a}) =
+        N_g (\alpha, \mathbf{a})
         (x - X_A)^{a_x} (y - Y_A)^{a_y} (z - Z_A)^{a_z}
-        \exp{-\alpha_i |\vec{r} - \vec{R}_A|^2}
+        \exp{-\alpha_i |\mathbf{r} - \mathbf{R}_A|^2}
 
-    where :math:`\vec{a} = (a_x, a_y, az)` is the angular momentum components in Cartesian
-    coordinates, :math:`\vec{R}_A` is the coordinate of the center :math:`A`, :math:`N_g` is the
+    where :math:`\mathbf{a} = (a_x, a_y, az)` is the angular momentum components in Cartesian
+    coordinates, :math:`\mathbf{R}_A` is the coordinate of the center :math:`A`, :math:`N_g` is the
     normalization constant of the primitive, and :math:`\alpha_i` is the exponent of the primitive.
 
     A **Cartesian contraction** is a linear combination of Cartesian primitives:
 
     .. math::
 
-        \phi_{\vec{a}, \vec{R}_A} (\mathbf{r}) =
-        N_{\phi} (\vec{a}, \vec{R}_A) \sum_i d_i g_i (\vec{r} | \vec{a}, \vec{R}_A)
+        \phi (\mathbf{r} | \mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) =
+        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_i g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
 
     where :math:`d_i` is the contraction coefficient of the primitive and :math:`N_{\phi}` is the
     normalization constant of the contraction.
 
     Note that the Cartesian contraction depends on the angular momentum components in the Cartesian
-    coordinate, :math:`\vec{a}`. At a given angular momentum, :math:`\ell`, we have
+    coordinate, :math:`\mathbf{a}`. At a given angular momentum, :math:`\ell`, we have
     :math:`\frac{(\ell + 1) (\ell + 2)}{2}` different components. We can linearly transform these
     contractions to obtain the :math:`2 \ell + 1` **spherical contractions**. Since we can convert
     the contractions in one coordinate system to another, we use the term **contraction** to
@@ -62,8 +62,8 @@ class GeneralizedContractionShell:
 
     .. math::
 
-        \phi_{j, \vec{a}, \vec{R}_A} (\mathbf{r}) =
-        N_{\phi} (\vec{a}, \vec{R}_A) \sum_i d_{ij} g_i (\vec{r} | \vec{a}, \vec{R}_A)
+        \phi_j (\mathbf{r} | \mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) =
+        N_{\phi} (\mathbf{R}_A, \mathbf{a}, \mathbf{d}, \boldsymbol{\alpha}) \sum_i d_{ij} g_i (\mathbf{r} | \mathbf{R}_A, \mathbf{a})
 
     Attributes
     ----------
@@ -73,11 +73,11 @@ class GeneralizedContractionShell:
         Angular momentum of the contractions.
 
         .. math::
-            \ell = \sum_i (\vec{a})_i = a_x + a_y + a_z
+            \ell = \sum_i (\mathbf{a})_i = a_x + a_y + a_z
 
     angmom_components_cart : np.ndarray(L, 3)
         The x, y, and z components of the angular momentum vectors
-        (:math:`\vec{a} = (a_x, a_y, a_z)` where :math:`a_x + a_y + a_z = \ell`).
+        (:math:`\mathbf{a} = (a_x, a_y, a_z)` where :math:`a_x + a_y + a_z = \ell`).
         `L` is the number of Cartesian contracted Gaussian functions for the given angular
         momentum, i.e. :math:`(\ell + 1) * (\ell + 2) / 2`.
         Property of `GeneralizedContractionShell`.
@@ -123,7 +123,7 @@ def __init__(self, angmom, coord, coeffs, exps, coord_type, icenter=None):
         angmom : int
             Angular momentum of the set of contractions.
 
-            :math:`\sum_i \vec{a} = a_x + a_y + a_z`
+            :math:`\sum_i \mathbf{a} = a_x + a_y + a_z`
 
         coord : np.ndarray(3,)
             Coordinate of the center of the contractions.
@@ -229,7 +229,7 @@ def angmom(self):
         angmom : int
             Angular momentum of the set of contractions.
 
-            :math:`\sum_i (\vec{a})_i = a_x + a_y + a_z`
+            :math:`\sum_i (\mathbf{a})_i = a_x + a_y + a_z`
 
         """
         return self._angmom
@@ -243,7 +243,7 @@ def angmom(self, angmom):
         angmom : int
             Angular momentum of the set of contractions.
 
-            :math:`\sum_i (\vec{a})_i = a_x + a_y + a_z`
+            :math:`\sum_i (\mathbf{a})_i = a_x + a_y + a_z`
 
         Raises
         ------
@@ -365,7 +365,7 @@ def angmom_components_cart(self):
         -------
         angmom_components_cart : np.ndarray(L, 3)
             The x, y, and z components of the angular momentum vectors
-            (:math:`\vec{a} = (a_x, a_y, a_z)` where :math:`a_x + a_y + a_z = \ell`).
+            (:math:`\mathbf{a} = (a_x, a_y, a_z)` where :math:`a_x + a_y + a_z = \ell`).
             `L` is the number of Cartesian contracted Gaussian functions for the given
             angular momentum, i.e. :math:`(\ell + 1) * (\ell + 2) / 2`
 
@@ -392,10 +392,9 @@ def angmom_components_sph(self):
         We define spherical components as regular solid harmonics :math:`R_{\ell, m}` such that
 
         .. math::
-            \begin{align}
-                R_{\ell, m} &= C_{\ell, m}, \quad &m = 0 ... \ell
-                R_{\ell, -m} &= S_{\ell, m}, \quad &m = 1 ... \ell
-            \end{align}
+
+                R_{\ell, m} = C_{\ell, m}, \quad m = 0 ... \ell\\
+                R_{\ell, -m} = S_{\ell, m}, \quad m = 1 ... \ell
 
         where :math:`C_{\ell, m}` and :math:`S_{\ell, m}` are Cosine- and Sine-like real
         regular solid harmonics (respectively).
@@ -439,7 +438,7 @@ def norm_prim_cart(self):
         For a Cartesian primitive with exponent :math:`\alpha_i`, the normalization constant is:
 
         .. math::
-           N(\alpha_i, \vec{a}) = \sqrt {
+           N(\alpha_i, \mathbf{a}) = \sqrt {
            \left(\frac{2\alpha_i}{\pi}\right)^\frac{3}{2}
            \frac{(4\alpha_i)^{a_x + a_y + a_z}}{(2a_x - 1)!! (2a_y - 1)!! (2a_z - 1)!!}}
 

gbasis/evals/density.py

@@ -65,6 +65,13 @@ def evaluate_density_using_evaluated_orbs(one_density_matrix, orb_eval):
 def evaluate_density(one_density_matrix, basis, points, transform=None, threshold=1.0e-8):
     r"""Return the density of the given basis set at the given points.
 
+    .. math::
+
+        \rho(\mathbf{r}) = \sum_{ij} \gamma_{ij} \phi_i(\mathbf{r}) \phi_j(\mathbf{r})
+
+    where :math:`\mathbf{r}` is the point at which the density is evaluated, :math:`\gamma_{ij}`
+    is the one-electron density matrix, and :math:`\phi_i` is the :math:`i`-th basis function.
+
     Parameters
     ----------
     one_density_matrix : np.ndarray(K_orbs, K_orbs)
@@ -116,23 +123,24 @@ def evaluate_deriv_reduced_density_matrix(
 
     .. math::
 
-        \left.
+        &\left.
         \frac{\partial^{p_x + p_y + p_z}}{\partial x_1^{p_x} \partial y_1^{p_y} \partial z_1^{p_z}}
         \frac{\partial^{q_x + q_y + q_z}}{\partial x_2^{q_x} \partial y_2^{q_y} \partial z_2^{q_z}}
         \gamma(\mathbf{r}_1, \mathbf{r}_2)
-        \right|_{\mathbf{r}_1 = \mathbf{r}_2 = \mathbf{r}_n} =
+        \right|_{\mathbf{r}_1 = \mathbf{r}_2 = \mathbf{r}}\\\\
+        &=
         \sum_{ij} \gamma_{ij}
         \left.
         \frac{\partial^{p_x + p_y + p_z}}{\partial x_1^{p_x} \partial y_1^{p_y} \partial z_1^{p_z}}
         \phi_i(\mathbf{r}_1)
-        \right|_{\mathbf{r}_1 = \mathbf{r}_n}
+        \right|_{\mathbf{r}_1 = \mathbf{r}}
         \left.
         \frac{\partial^{q_x + q_y + q_z}}{\partial x_2^{q_x} \partial y_2^{q_y} \partial z_2^{q_z}}
         \phi_j(\mathbf{r}_2)
-        \right|_{\mathbf{r}_2 = \mathbf{r}_n}
+        \right|_{\mathbf{r}_2 = \mathbf{r}}
 
     where :math:`\mathbf{r}_1` is the first point, :math:`\mathbf{r}_2` is the second point, and
-    :math:`\mathbf{r}_n` is the point at which the derivative is evaluated.
+    :math:`\mathbf{r}` is the point at which the derivative is evaluated.
 
     Parameters
     ----------
@@ -202,6 +210,24 @@ def evaluate_deriv_density(
 ):
     r"""Return the derivative of density of the given transformed basis set at the given points.
 
+    .. math::
+
+            &\frac{\partial^{L_x + L_y + L_z}}{\partial x^{L_x} \partial y^{L_y} \partial z^{L_z}}
+            \rho(\mathbf{r})\\\\
+            &=
+            \sum_{l_x=0}^{L_x} \sum_{l_y=0}^{L_y} \sum_{l_z=0}^{L_z}
+            \binom{L_x}{l_x} \binom{L_y}{l_y} \binom{L_z}{l_z}
+            \sum_{ij} \gamma_{ij}
+            \frac{\partial^{l_x + l_y + l_z} \rho(\mathbf{r})}{\partial x^{l_x} \partial y^{l_y} \partial z^{l_z}}
+            \frac{
+                \partial^{L_x + L_y + L_z - l_x - l_y - l_z} \rho(\mathbf{r})
+            }{
+                \partial x^{L_x - l_x} \partial y^{L_y - l_y} \partial z^{L_z - l_z}
+            }
+
+    where :math:`L_x, L_y, L_z` are the orders of the derivative relative to the :math:`x, y, \text{and} z` components,
+    respectively.
+
     Parameters
     ----------
     orders : np.ndarray(3,)
@@ -286,6 +312,16 @@ def evaluate_density_gradient(
 ):
     r"""Return the gradient of the density evaluated at the given points.
 
+    .. math::
+
+            \nabla \rho(\mathbf{r})
+            =
+            \begin{bmatrix}
+            \frac{\partial}{\partial x} \rho(\mathbf{r})\\\\
+            \frac{\partial}{\partial y} \rho(\mathbf{r})\\\\
+            \frac{\partial}{\partial z} \rho(\mathbf{r})
+            \end{bmatrix}      
+
     Parameters
     ----------
     one_density_matrix : np.ndarray(K_orb, K_orb)
@@ -353,6 +389,14 @@ def evaluate_density_laplacian(
 ):
     r"""Return the Laplacian of the density evaluated at the given points.
 
+    .. math::
+
+            \nabla^2 \rho(\mathbf{r})
+            =
+            \frac{\partial^2}{\partial x^2} \rho(\mathbf{r})
+            + \frac{\partial^2}{\partial y^2} \rho(\mathbf{r})
+            + \frac{\partial^2}{\partial z^2} \rho(\mathbf{r})
+
     Parameters
     ----------
     one_density_matrix : np.ndarray(K_orb, K_orb)
@@ -444,6 +488,22 @@ def evaluate_density_hessian(
 ):
     r"""Return the Hessian of the density evaluated at the given points.
 
+    .. math::
+
+        H[\rho(\mathbf{r})]
+        =
+        \begin{bmatrix}
+            \frac{\partial^2}{\partial x^2} \rho(\mathbf{r}) &
+            \frac{\partial^2}{\partial x \partial y} \rho(\mathbf{r}) &
+            \frac{\partial^2}{\partial x \partial z} \rho(\mathbf{r})\\\\
+            \frac{\partial^2}{\partial y \partial x} \rho(\mathbf{r}) &
+            \frac{\partial^2}{\partial y^2} \rho(\mathbf{r})&
+            \frac{\partial^2}{\partial y \partial z} \rho(\mathbf{r})\\\\
+            \frac{\partial^2}{\partial z \partial x} \rho(\mathbf{r}) &
+            \frac{\partial^2}{\partial z \partial y} \rho(\mathbf{r})&
+            \frac{\partial^2}{\partial z^2} \rho(\mathbf{r})\\
+        \end{bmatrix}
+
     Parameters
     ----------
     one_density_matrix : np.ndarray(K_orb, K_orb)
@@ -566,15 +626,17 @@ def evaluate_posdef_kinetic_energy_density(
 
     .. math::
 
-        t_+ (\mathbf{r}_n)
+        \begin{split}
+        t_+ (\mathbf{r})
         &= \frac{1}{2} \left.
-          \nabla_{\mathbf{r}} \cdot \nabla_{\mathbf{r}'} \gamma(\mathbf{r}, \mathbf{r}')
-        \right|_{\mathbf{r} = \mathbf{r}' = \mathbf{r}_n}\\
+          \nabla_{\mathbf{r}_1} \cdot \nabla_{\mathbf{r}_2} \gamma(\mathbf{r}_1, \mathbf{r}_2)
+        \right|_{\mathbf{r}_1 = \mathbf{r}_2 = \mathbf{r}}\\
         &= \frac{1}{2} \left(
-          \frac{\partial^2}{\partial x \partial x'} \gamma(\mathbf{r}, \mathbf{r}')
-          + \frac{\partial^2}{\partial y \partial y'} \gamma(\mathbf{r}, \mathbf{r}')
-          + \frac{\partial^2}{\partial z \partial z'} \gamma(\mathbf{r}, \mathbf{r}')
-        \right)_{\mathbf{r} = \mathbf{r}' = \mathbf{r}_n}\\
+          \frac{\partial^2}{\partial x_1 \partial x_2} \gamma(\mathbf{r}_1, \mathbf{r}_2)
+          + \frac{\partial^2}{\partial y_1 \partial y_2} \gamma(\mathbf{r}_1, \mathbf{r}_2)
+          + \frac{\partial^2}{\partial z_1 \partial z_2} \gamma(\mathbf{r}_1, \mathbf{r}_2)
+        \right)_{\mathbf{r}_1 = \mathbf{r}_2 = \mathbf{r}}\\
+        \end{split}
 
     Parameters
     ----------
@@ -641,7 +703,7 @@ def evaluate_general_kinetic_energy_density(
 
     .. math::
 
-        t_{\alpha} (\mathbf{r}_n) = t_+(\mathbf{r}_n) + \alpha \nabla^2 \rho(\mathbf{r}_n)
+        t_{\alpha} (\mathbf{r}) = t_+(\mathbf{r}) + \alpha \nabla^2 \rho(\mathbf{r})
 
     where :math:`t_+` is the positive definite kinetic energy density.
 

gbasis/evals/electrostatic_potential.py

@@ -14,6 +14,17 @@ def electrostatic_potential(
 ):
     r"""Return the electrostatic potentials of the basis set in the Cartesian form.
 
+    .. math::
+
+        - \left(
+            - \sum_A \frac{Z_A}{|\mathbf{R}_C - \mathbf{R}_A|}
+            + \sum_{ab} \gamma_{ab} \int \phi_a(\mathbf{r}) \frac{-1}{|\mathbf{r} - \mathbf{R}_C|} \phi_b(\mathbf{r}) d\mathbf{r}
+        \right)
+
+    where :math:`\mathbf{R}_C` is the coordinate of a unitary point charge, :math:`\mathbf{R}_A` is the
+    coordinate of the nucleus :math:`A`, :math:`Z_A` its charge, and :math:`\gamma_{ab}` is the
+    one-electron density matrix.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell

gbasis/evals/eval_deriv.py

@@ -147,6 +147,17 @@ def evaluate_deriv_basis(
 ):
     r"""Evaluate the derivative of the basis set in the given coordinate system at the given points.
 
+    The derivative (to arbitrary orders) of a basis function is given by:
+
+    .. math::
+
+        \frac{\partial^{m_x + m_y + m_z}}{\partial x^{m_x} \partial y^{m_y} \partial z^{m_z}}
+        \phi (\mathbf{r})
+
+    where :math:`m_x, m_y, m_z` are the orders of the derivative with respect to x, y, and z, 
+    :math:`\phi` is the basis function (a generalized contraction shell), and :math:`\mathbf{r}_n` are
+    the coordinate of the points in space where the basis function is evaluated.
+
     Parameters
     ----------
     basis : list/tuple of GeneralizedContractionShell