From b1884d105ac5a0b8be494eb58ef31ae0c2545351 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Wed, 18 May 2022 12:50:34 -0400
Subject: [PATCH 01/11] Add second derivative and laplacian to atomic dens

---
 bfit/density.py | 110 +++++++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 104 insertions(+), 6 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index ad92da5..8a4b5ac 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -249,7 +249,7 @@ def radial_slater_orbital(exponent, number, points):
         slater = norm.T * pref * np.exp(-exponent * points).T
         return slater
 
-    def phi_matrix(self, points, deriv=False):
+    def phi_matrix(self, points, deriv=None):
         r"""
         Compute the linear combination of Slater-type orbitals on the given points.
 
@@ -271,8 +271,9 @@ def phi_matrix(self, points, deriv=False):
         ----------
         points : ndarray, (N,)
             The radial grid points.
-        deriv : bool
-            If true, use the derivative of the slater-orbitals.
+        deriv : int, optional
+            If it is one, return the derivative of the slater-orbitals.
+            If it is two, return the second derivative of the slater-orbitals.
 
         Returns
         -------
@@ -287,12 +288,21 @@ def phi_matrix(self, points, deriv=False):
           zero instead. See "derivative_radial_slater_type_orbital".
 
         """
+        if deriv is not None:
+            if deriv not in [1, 2]:
+                raise ValueError(
+                    f"The deriv parameter {deriv} should be either one or two. Higher-order"
+                    f"is not supported."
+                )
+
         # compute orbital composed of a linear combination of Slater
         phi_matrix = np.zeros((len(points), len(self.orbitals)))
         for index, orbital in enumerate(self.orbitals):
             exps, number = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]]
-            if deriv:
+            if deriv == 1:
                 slater = self.derivative_radial_slater_type_orbital(exps, number, points)
+            elif deriv == 2:
+                slater = self.second_derivative_radial_slater_type_orbital(exps, number, points)
             else:
                 slater = self.radial_slater_orbital(exps, number, points)
             phi_matrix[:, index] = np.dot(slater, self.orbitals_coeff[orbital]).ravel()
@@ -387,7 +397,7 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
         Returns
         -------
         slater : ndarray, (N, M)
-            The Slater-type orbitals evaluated on the :math:`N` grid points.
+            First derivative of Slater-type orbitals evaluated on the :math:`N` grid points.
 
         Notes
         -----
@@ -407,6 +417,63 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
         deriv = deriv_pref * slater
         return deriv
 
+    @staticmethod
+    def second_derivative_radial_slater_type_orbital(exponent, number, points):
+        r"""
+        Compute the derivative of the radial component of Slater-type orbital on the given points.
+
+        The derivative of the Slater-type orbital is defined as:
+
+        .. math::
+            \frac{d^2 R(r)}{dr^2} = \bigg[
+                \bigg(\frac{n-1}{r} - \alpha \bigg)^2 - \frac{n-1}{r^2} \bigg]
+                 N r^{n-1} e^{- \alpha r},
+
+        where
+        :math:`n` is the principal quantum number of that orbital,
+        :math:`N` is the normalizing constant,
+        :math:`r` is the radial point, distance to the origin, and
+        :math:`\alpha` is the zeta exponent of that orbital.
+
+        Parameters
+        ----------
+        exponent : ndarray, (M, 1)
+            The zeta exponents of Slater orbitals.
+        number : ndarray, (M, 1)
+            The principle quantum numbers of Slater orbitals.
+        points : ndarray, (N,)
+            The radial grid points. If points contain zero, then it is undefined at those
+            points and set to zero.
+
+        Returns
+        -------
+        slater : ndarray, (N, M)
+            Second derivative of Slater-type orbitals evaluated on the :math:`N` grid points.
+
+        Notes
+        -----
+        - At r = 0, the derivative is undefined and this function returns zero instead.
+
+        References
+        ----------
+        See wikipedia page on "Slater-Type orbitals".
+
+        """
+        slater = SlaterAtoms.radial_slater_orbital(exponent, number, points)
+        # Consider the case when dividing by zero.
+        with np.errstate(divide='ignore'):
+            # derivative
+            # [n - 1] / r - alpha
+            deriv_pref = (number.T - 1.) / np.reshape(points, (points.shape[0], 1)) - exponent.T
+            deriv_pref[np.abs(points) < 1e-10, :] = 0.0
+
+            # [n - 1] / r^2
+            deriv_pref2 = (number.T - 1.) / np.reshape(points**2.0, (points.shape[0], 1))
+            deriv_pref2[np.abs(points) < 1e-10, :] = 0.0
+
+        deriv = (deriv_pref**2.0 - deriv_pref2) * slater
+        return deriv
+
     def positive_definite_kinetic_energy(self, points):
         r"""
         Positive definite or Lagrangian kinetic energy density.
@@ -463,6 +530,37 @@ def derivative_density(self, points):
             The derivative of atomic density on the grid points.
 
         """
-        factor = self.phi_matrix(points) * self.phi_matrix(points, deriv=True)
+        factor = self.phi_matrix(points) * self.phi_matrix(points, deriv=1)
         derivative = np.dot(2. * factor, self.orbitals_occupation).ravel() / (4 * np.pi)
         return derivative
+
+    def laplacian_of_atomic_density(self, points):
+        r"""
+        Return the Laplacian of the atomic density on a set of points.
+
+        Parameters
+        ----------
+        points : ndarray(N,)
+            The :math:`N` radial grid points.
+
+        Returns
+        -------
+        laplacian : ndarray(N,)
+            The Laplacian of the atomic density on the grid points.
+
+        """
+        molecular_orb = self.phi_matrix(points)
+        molecular_orb_deriv = self.phi_matrix(points, deriv=1)
+        molecular_orb_sec_deriv = self.phi_matrix(points, deriv=2)
+
+        # phi_i * phi_i^\prime / r
+        with np.errstate(divide='ignore'):
+            factor = 2.0 * molecular_orb * molecular_orb_deriv / \
+                     np.reshape(points, (points.shape[0], 1))
+            factor[np.abs(points) < 1e-10] = 0.0
+
+        factor = (
+            molecular_orb_deriv**2.0 + factor + molecular_orb * molecular_orb_sec_deriv
+        )
+        laplacian = np.dot(2. * factor, self.orbitals_occupation).ravel() / (4 * np.pi)
+        return laplacian

From bf96c9c8cef776b8bf6780efed3810c8b4f489c6 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Wed, 18 May 2022 12:51:02 -0400
Subject: [PATCH 02/11] Add tests for second deriv of molecular orbital

---
 bfit/test/test_density.py | 28 +++++++++++++++++++++++++++-
 1 file changed, 27 insertions(+), 1 deletion(-)

diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index 1b76a7d..d0d7eea 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -83,6 +83,21 @@ def test_derivative_slater_type_orbital_be():
     assert_almost_equal(orbitals, expected, decimal=6)
 
 
+def test_second_derivative_of_slater_type_orbital_against_finite_difference():
+    r"""Test second derivative of radial Slater-type orbital against finite-difference."""
+    # load Be atomic wave function
+    be = SlaterAtoms("Be")
+
+    # check values of a single orbital at various points
+    for r in np.arange(0.1, 10., 0.1):
+        exponent, number, pt = np.array([[12.683501]]), np.array([[1]]), np.array([r])
+        actual = be.second_derivative_radial_slater_type_orbital(exponent, number, pt)
+
+        first_deriv = lambda x: slater(exponent[0, 0], number[0, 0], x, derivative=True)
+        desired = (first_deriv(r + 1e-8) - first_deriv(r)) / 1e-8
+        assert_almost_equal(actual, desired, decimal=3)
+
+
 def test_positive_definite_kinetic_energy_he():
     r"""Test integral of kinetic energy density of helium against actual value."""
     # load he atomic wave function
@@ -143,7 +158,7 @@ def test_phi_derivative_lcao_b():
     # load Be atomic wave function
     b = SlaterAtoms("b")
     # check the values of the phi_matrix at point 1.0
-    phi_matrix = b.phi_matrix(np.array([1]), deriv=True)
+    phi_matrix = b.phi_matrix(np.array([1]), deriv=1)
 
     def _slater_deriv(r):
         # compute expected value of 1S
@@ -173,6 +188,17 @@ def _slater_deriv(r):
     assert_almost_equal(phi_matrix[1, :], _slater_deriv(2.), decimal=4)
     assert_almost_equal(phi_matrix[2, :], _slater_deriv(3.), decimal=4)
 
+    # Test second derivative of phi at point 1.0, 2.0, 3.0
+    pts = np.arange(0.5, 5., 0.1)
+    actual = b.phi_matrix(pts, deriv=2)
+    for i, r in enumerate(pts):
+        desired = (np.array(_slater_deriv(r + 1e-8)) - np.array(_slater_deriv(r))) / 1e-8
+        assert_almost_equal(actual[i], desired, decimal=3)
+
+    # Test raises error if value is not correct.
+    assert_raises(ValueError, b.phi_matrix, points=pts, deriv=3)
+    assert_raises(ValueError, b.phi_matrix, points=pts, deriv=b)
+
 
 def test_coeff_matrix_be():
     r"""Test the coefficients of Beryllium is correctly parsed."""

From 2b1eeb78bcbaf23afca385a4a7938654768d1213 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Wed, 18 May 2022 12:52:04 -0400
Subject: [PATCH 03/11] Add hydrogen test for laplacian atomic dens

---
 bfit/test/test_density.py | 14 ++++++++++++++
 1 file changed, 14 insertions(+)

diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index d0d7eea..4ddf39c 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -442,6 +442,20 @@ def test_kinetic_energy_heavy_element_ce():
     assert_almost_equal(integral, c.kinetic_energy, decimal=2)
 
 
+def test_laplacian_of_hydrogen_slater():
+    r"""Test laplacian of hydrogen Slater."""
+    # Slater is N e^(-r)
+    h = SlaterAtoms("h")
+
+    pts = np.arange(0.01, 10., 0.1)
+    lap = h.laplacian_of_atomic_density(pts)
+
+    # Analytic derive it
+    solution = np.exp(-2.0 * pts) * (pts - 1) * 16.0 / pts
+    solution /= (4.0 * np.pi)
+    assert_almost_equal(lap, solution, decimal=4)
+
+
 def test_raises():
     r"""Test raises error of SlaterAtoms."""
     assert_raises(TypeError, SlaterAtoms, 25)

From deaa88444b4d924cf090e8903e476f36e26ce0f0 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:24:46 -0400
Subject: [PATCH 04/11] Add option for normalization to Slater orbital

---
 bfit/density.py | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index 8a4b5ac..2e01a41 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -204,7 +204,7 @@ def orbitals_cusp(self):
         return self._orbitals_cusp
 
     @staticmethod
-    def radial_slater_orbital(exponent, number, points):
+    def radial_slater_orbital(exponent, number, points, normalized=True):
         r"""
         Compute the radial component of Slater-type orbitals on the given points.
 
@@ -227,6 +227,8 @@ def radial_slater_orbital(exponent, number, points):
             The principal quantum numbers :math:`n` of :math:`M` Slater orbitals.
         points : ndarray, (N,)
             The radial :math:`r` grid points.
+        normalized : bool
+            If true, returns the normalization constant :math:`N`.
 
         Returns
         -------
@@ -242,11 +244,16 @@ def radial_slater_orbital(exponent, number, points):
         """
         if points.ndim != 1:
             raise ValueError("The argument point should be a 1D array.")
-        # compute norm & pre-factor
-        norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number))
-        pref = np.power(points, number - 1).T
+        # compute pre-factor
+        with np.errstate(divide='ignore'):
+            pref = np.power(points, number - 1, dtype=np.float64).T
         # compute slater orbitals
-        slater = norm.T * pref * np.exp(-exponent * points).T
+        slater = pref * np.exp(-exponent * points).T
+        # compute normalization
+        if normalized:
+            norm = np.power(2. * exponent, number) * \
+                   np.sqrt((2. * exponent) / factorial(2. * number))
+            slater *= norm.T
         return slater
 
     def phi_matrix(self, points, deriv=None):

From 8d90e7b23f8e0e298f54e5a3148680befdd75fb8 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:27:31 -0400
Subject: [PATCH 05/11] Fix taking derivative of slater orbital

When n=1, the first component (n-1)/r does not exist
---
 bfit/density.py           | 24 ++++++++++++++----------
 bfit/test/test_density.py | 13 +++++++++++--
 2 files changed, 25 insertions(+), 12 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index 2e01a41..1c2761f 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -406,22 +406,26 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
         slater : ndarray, (N, M)
             First derivative of Slater-type orbitals evaluated on the :math:`N` grid points.
 
-        Notes
-        -----
-        - At r = 0, the derivative is undefined and this function returns zero instead.
-
         References
         ----------
         See wikipedia page on "Slater-Type orbitals".
 
         """
-        slater = SlaterAtoms.radial_slater_orbital(exponent, number, points)
+        norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number))
+        slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True)
+        # Calculate the un-normalized Slater with number less than one
+        slater_minus = SlaterAtoms.radial_slater_orbital(
+            exponent, number - 1, points, normalized=False
+        )
         # Consider the case when dividing by zero.
-        with np.errstate(divide='ignore'):
-            # derivative
-            deriv_pref = (number.T - 1.) / np.reshape(points, (points.shape[0], 1)) - exponent.T
-            deriv_pref[np.abs(points) < 1e-10, :] = 0.0
-        deriv = deriv_pref * slater
+        deriv = np.zeros((len(points), number.shape[0]))
+        # Compute -\alpha * slater
+        deriv -= exponent.T * slater
+        # Compute (n-1)*slater/r, note that this only occurs when n!=1.
+        deriv_pref = norm.T * slater_minus * (number.T - 1.0)
+        i_numb_one = np.where(number == 1)[0]
+        deriv_pref[:, i_numb_one] = 0.0
+        deriv += deriv_pref
         return deriv
 
     @staticmethod
diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index 4ddf39c..938a1eb 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -25,7 +25,7 @@
 from bfit.density import SlaterAtoms
 from bfit.grid import ClenshawRadialGrid
 import numpy as np
-from numpy.testing import assert_almost_equal, assert_equal, assert_raises
+from numpy.testing import assert_allclose, assert_almost_equal, assert_equal, assert_raises
 import scipy
 
 
@@ -35,6 +35,8 @@ def slater(e, n, r, derivative=False):
     slater = norm * np.power(r, n - 1) * np.exp(-e * r)
 
     if derivative:
+        if n == 1:
+            return - e * slater
         return (((n - 1) / r) - e) * slater
     return slater
 
@@ -49,6 +51,12 @@ def test_slater_type_orbital_be():
     # check values of a single orbital at r=2.0
     orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[2]]), np.array([2]))
     assert_almost_equal(orbital, slater(0.821620, 2, 2.0), decimal=6)
+    # check values of a single orbital at r=0.0 and n=2
+    orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[2]]), np.array([0.0]))
+    assert_almost_equal(orbital, slater(0.821620, 2, 0.0), decimal=6)
+    # check values of a single orbital at r=0.0 and n=1
+    orbital = be.radial_slater_orbital(np.array([[0.821620]]), np.array([[1]]), np.array([0.0]))
+    assert_almost_equal(orbital, slater(0.821620, 1, 0.0), decimal=6)
     # check value of tow orbitals at r=1.0 & r=2.0
     exps, nums = np.array([[12.683501], [0.821620]]), np.array([[1], [2]])
     orbitals = be.radial_slater_orbital(exps, nums, np.array([1., 2.]))
@@ -72,7 +80,8 @@ def test_derivative_slater_type_orbital_be():
     # check value of single orbital at r = 0.0
     orbital = be.derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]),
                                                        np.array([0.]))
-    assert_almost_equal(orbital, 0.0, decimal=6)
+    actual = np.power(2. * 0.821620, 2) * np.sqrt((2. * 0.821620) / (4 * 3 * 2))
+    assert_almost_equal(orbital, actual, decimal=6)
     # check value of tow orbitals at r=1.0 & r=2.0
     exps, nums = np.array([[12.683501], [0.821620]]), np.array([[1], [2]])
     orbitals = be.derivative_radial_slater_type_orbital(exps, nums, np.array([1., 2.]))

From fe81c2f92ab4d4de8d68e98315c6e67c93d3f6b3 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:28:37 -0400
Subject: [PATCH 06/11] Fix taking second derivative of Slater-orbital

Certain components of the derivatives aren't needed
when n=1,2
---
 bfit/density.py           | 39 +++++++++++++++++++++++++--------------
 bfit/test/test_density.py |  2 +-
 2 files changed, 26 insertions(+), 15 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index 1c2761f..dab12d0 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -437,7 +437,7 @@ def second_derivative_radial_slater_type_orbital(exponent, number, points):
 
         .. math::
             \frac{d^2 R(r)}{dr^2} = \bigg[
-                \bigg(\frac{n-1}{r} - \alpha \bigg)^2 - \frac{n-1}{r^2} \bigg]
+                \frac{(n-1)(n-2)}{r^2} - \frac{2\alpha (n-1)}{r} + \alpha^2 \bigg]
                  N r^{n-1} e^{- \alpha r},
 
         where
@@ -463,26 +463,37 @@ def second_derivative_radial_slater_type_orbital(exponent, number, points):
 
         Notes
         -----
-        - At r = 0, the derivative is undefined and this function returns zero instead.
+        - When :math:`n=1,2` and :math:`r=0`, then the derivative is un-defined and this
+          function returns 0 instead.
 
         References
         ----------
         See wikipedia page on "Slater-Type orbitals".
 
         """
-        slater = SlaterAtoms.radial_slater_orbital(exponent, number, points)
-        # Consider the case when dividing by zero.
+        norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number))
+        slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True)
+        # Calculate the un-normalized Slater with number less than one
         with np.errstate(divide='ignore'):
-            # derivative
-            # [n - 1] / r - alpha
-            deriv_pref = (number.T - 1.) / np.reshape(points, (points.shape[0], 1)) - exponent.T
-            deriv_pref[np.abs(points) < 1e-10, :] = 0.0
-
-            # [n - 1] / r^2
-            deriv_pref2 = (number.T - 1.) / np.reshape(points**2.0, (points.shape[0], 1))
-            deriv_pref2[np.abs(points) < 1e-10, :] = 0.0
-
-        deriv = (deriv_pref**2.0 - deriv_pref2) * slater
+            slater_minus = SlaterAtoms.radial_slater_orbital(
+                exponent, number - 1, points, normalized=False
+            )
+            slater_minus_two = SlaterAtoms.radial_slater_orbital(
+                exponent, number - 2, points, normalized=False
+            )
+        # Compute \alpha^2 * slater
+        deriv = exponent.T**2.0 * slater
+        # Compute 2 \alpha (n-1) * slater / r, note doesn't occur when n=1
+        deriv_pref = 2.0 * exponent.T * norm.T * slater_minus * (number.T - 1.0)
+        i_numb_one = np.where(number == 1)[0]
+        deriv_pref[:, i_numb_one] = 0.0
+        deriv -= deriv_pref
+        # Compute (n-1)(n-2) * slater / r^2, note doesn't occur when n=1,2
+        deriv_pref = norm.T * slater_minus_two * (number.T - 1.0) * (number.T - 2.0)
+        i_numb_one = np.where(number == 1)[0]
+        deriv_pref[:, i_numb_one] = 0.0
+        deriv_pref[:, np.where(number == 2)[0]] = 0.0
+        deriv += deriv_pref
         return deriv
 
     def positive_definite_kinetic_energy(self, points):
diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index 938a1eb..f89b803 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -98,7 +98,7 @@ def test_second_derivative_of_slater_type_orbital_against_finite_difference():
     be = SlaterAtoms("Be")
 
     # check values of a single orbital at various points
-    for r in np.arange(0.1, 10., 0.1):
+    for r in np.arange(0.0, 10., 0.1):
         exponent, number, pt = np.array([[12.683501]]), np.array([[1]]), np.array([r])
         actual = be.second_derivative_radial_slater_type_orbital(exponent, number, pt)
 

From 1799705a6bc3d88b84c3987863cb9f7a97eb27f4 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:29:57 -0400
Subject: [PATCH 07/11] Fix taking kinetic energy

When r=0, n=1, then the kinetic energy should return infinity
---
 bfit/density.py           | 38 ++++++++++++++++++++++++++++++--------
 bfit/test/test_density.py | 28 +++++++++++++++++-----------
 2 files changed, 47 insertions(+), 19 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index dab12d0..ed19387 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -500,6 +500,15 @@ def positive_definite_kinetic_energy(self, points):
         r"""
         Positive definite or Lagrangian kinetic energy density.
 
+        .. math::
+            K(r) = \sum_i n_i \bigg[ \bigg(\frac{d \sum_j N_j c^i_j R_{n_j}}{dr} \bigg)^2
+             + \frac{l_i (l_i + 1}{r^2} \bigg(\sum_j N_j c^i_j R_{n_j} \bigg)^2 \bigg],
+
+        where
+        :math:`n` is the principal quantum number of that orbital,
+        :math:`N` is the normalizing constant, andx
+        :math:`r` is the radial point, distance to the origin
+
         Parameters
         ----------
         points : ndarray,(N,)
@@ -510,6 +519,10 @@ def positive_definite_kinetic_energy(self, points):
         energy : ndarray, (N,)
             The kinetic energy on the grid points.
 
+        Notes
+        -----
+        - When :math:`n=1` and :math:`r=0`, then kinetic energy is undefined and nan is returned.
+
         """
         phi_matrix = np.zeros((len(points), len(self.orbitals)))
         angular = {
@@ -523,14 +536,23 @@ def positive_definite_kinetic_energy(self, points):
             phi_matrix[:, index] = np.ravel(np.dot(deriv_radial, self.orbitals_coeff[orbital])**2.0)
 
             # Calculate del^2 of spherical component
-            slater = SlaterAtoms.radial_slater_orbital(exps, numbers, points)
-            with np.errstate(divide='ignore'):
-                gradient_sph = (
-                    angular[orbital[1]] *
-                    np.ravel(np.dot(slater, self.orbitals_coeff[orbital]))**2.0 /
-                    points**2.0
-                )
-                gradient_sph[np.abs(points) < 1e-10] = 0.0
+            norm = np.power(2. * exps, numbers) * np.sqrt((2. * exps) / factorial(2. * numbers))
+            # Take unnormalized slater with number n-1, this is needed to remove divide by r^2
+            slater_minus_one = SlaterAtoms.radial_slater_orbital(
+                exps, numbers - 1, points, normalized=False
+            )
+            deriv_pref = norm.T * slater_minus_one
+            # When r=0 and n = 1, then the derivative is undefined and this returns infinity.
+            i_r_zero = np.where(np.abs(points) == 0.0)[0]
+            i_numb_one = np.where(numbers[0] == 1)[0]
+            indices = np.array([[x, y] for x in i_r_zero for y in i_numb_one])
+            if len(indices) != 0:  # if-statement needed to remove numpy warning using list
+                deriv_pref[indices] = np.inf
+            # Sum the slater orbital multipled with their coefficients
+            gradient_sph = (
+                angular[orbital[1]] *
+                np.ravel(np.dot(deriv_pref, self.orbitals_coeff[orbital]))**2.0
+            )
             phi_matrix[:, index] += gradient_sph
 
         orb_occs = self.orbitals_occupation
diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index f89b803..acc5818 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -112,21 +112,27 @@ def test_positive_definite_kinetic_energy_he():
     # load he atomic wave function
     he = SlaterAtoms("he")
     # compute density on an equally distant grid
-    grid = ClenshawRadialGrid(4, 30000, 35000)
+    grid = ClenshawRadialGrid(4, 30000, 35000, include_origin=False)
     energ = he.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert np.all(np.abs(integral - he.kinetic_energy) < 1e-5)
+    # Test at r=0, that it returns nan
+    energ = he.positive_definite_kinetic_energy(np.array([0.]))
+    assert np.all(np.isnan(energ[0]))
 
 
 def test_positive_definite_kinetic_energy_li():
     r"""Test integral of kinetic energy density of lithium against actual value."""
     # load be atomic wave function
-    be = SlaterAtoms("li")
+    li = SlaterAtoms("li")
     # compute density on an equally distant grid
-    grid = ClenshawRadialGrid(3, 20000, 35000)
-    energ = be.positive_definite_kinetic_energy(grid.points)
+    grid = ClenshawRadialGrid(3, 20000, 35000, include_origin=False)
+    energ = li.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
-    assert np.all(np.abs(integral - be.kinetic_energy) < 1e-5)
+    assert np.all(np.abs(integral - li.kinetic_energy) < 1e-5)
+    # Test at r=0, that it returns nan
+    energ = li.positive_definite_kinetic_energy(np.array([0.]))
+    assert np.all(np.isnan(energ[0]))
 
 
 def test_positive_definite_kinetic_energy_c():
@@ -134,7 +140,7 @@ def test_positive_definite_kinetic_energy_c():
     # load be atomic wave function
     be = SlaterAtoms("c")
     # compute density on an equally distant grid
-    grid = ClenshawRadialGrid(6, 20000, 35000)
+    grid = ClenshawRadialGrid(6, 20000, 35000, include_origin=False, dtype=np.float64)
     energ = be.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert np.all(np.abs(integral - be.kinetic_energy) < 1e-5)
@@ -145,7 +151,7 @@ def test_positive_definite_kinetic_energy_p():
     # load be atomic wave function
     be = SlaterAtoms("p")
     # compute density on ClenshawCurtis Grid
-    grid = ClenshawRadialGrid(15, 20000, 35000)
+    grid = ClenshawRadialGrid(15, 20000, 35000, include_origin=False, dtype=np.longlong)
     energ = be.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert np.all(np.abs(integral - be.kinetic_energy) < 1e-3)
@@ -156,7 +162,7 @@ def test_positive_definite_kinetic_energy_ag():
     # load c atomic wave function
     adens = SlaterAtoms("ag")
     # compute density on ClenshawCurtis Grid
-    grid = ClenshawRadialGrid(47, 20000, 35000)
+    grid = ClenshawRadialGrid(47, 20000, 35000, include_origin=False, dtype=np.longlong)
     energ = adens.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert np.all(np.abs(integral - adens.kinetic_energy) < 1e-3)
@@ -406,13 +412,13 @@ def test_atomic_density_heavy_rn():
 def test_kinetic_energy_cation_anion_c():
     r"""Test integral of kinetic energy density of cation/anion of carbon."""
     c = SlaterAtoms("c", cation=True)
-    grid = ClenshawRadialGrid(6, 20000, 35000)
+    grid = ClenshawRadialGrid(6, 20000, 35000, include_origin=False)
     energ = c.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert_almost_equal(integral, c.kinetic_energy, decimal=6)
 
     c = SlaterAtoms("c", anion=True)
-    grid = ClenshawRadialGrid(7, 20000, 35000)
+    grid = ClenshawRadialGrid(7, 20000, 35000, include_origin=False)
     energ = c.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert_almost_equal(integral, c.kinetic_energy, decimal=5)
@@ -445,7 +451,7 @@ def test_derivative_electron_density_cr():
 def test_kinetic_energy_heavy_element_ce():
     r"""Test integral of kinetic energy of cesium."""
     c = SlaterAtoms("ce")
-    grid = ClenshawRadialGrid(55, 10000, 30000)
+    grid = ClenshawRadialGrid(55, 10000, 30000, include_origin=False)
     energ = c.positive_definite_kinetic_energy(grid.points)
     integral = 4.0 * np.pi * grid.integrate(energ * grid.points**2.0)
     assert_almost_equal(integral, c.kinetic_energy, decimal=2)

From be835d27c15bc1d75a538f428c0ce0dc9e6ddbb0 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:31:15 -0400
Subject: [PATCH 08/11] Fix taking Laplacian of atomic density

When r=0, n=1, the Laplacian should be infinity
---
 bfit/density.py           | 43 ++++++++++++++++++++++++++++++++++-----
 bfit/test/test_density.py | 35 ++++++++++++++++---------------
 2 files changed, 57 insertions(+), 21 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index ed19387..2e2e519 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -582,6 +582,22 @@ def laplacian_of_atomic_density(self, points):
         r"""
         Return the Laplacian of the atomic density on a set of points.
 
+        The Laplacian in radial coordinates only is:
+
+        .. math::
+            \Delta f = \frac{1}{r^2}\frac{\partial }{\partial r}\bigg(r^2 \frac{df}{dr} \bigg).
+
+        Letting f be the atomic density we have the following equation:
+
+        .. math::
+            \Delta f = 2 \bigg[\sum n_i \big( \frac{d \phi_i}{dr}^2 +
+                    \frac{2 \phi \frac{d\phi_i}{dr}}{r} +
+                    \phi_i \frac{d^2 \phi_i}{dr^2} \big)
+            \bigg],
+
+        where
+        :math:`\phi_i` is the ith molecular orbital with occupation number :math:`n_i`.
+
         Parameters
         ----------
         points : ndarray(N,)
@@ -599,12 +615,29 @@ def laplacian_of_atomic_density(self, points):
 
         # phi_i * phi_i^\prime / r
         with np.errstate(divide='ignore'):
-            factor = 2.0 * molecular_orb * molecular_orb_deriv / \
-                     np.reshape(points, (points.shape[0], 1))
-            factor[np.abs(points) < 1e-10] = 0.0
-
+            # Absorb phi_i / r together
+            phi_i_r = np.zeros((len(points), len(self.orbitals)))
+            for index, orbital in enumerate(self.orbitals):
+                exps, numbers = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]]
+                # Calculate slater divided by r
+                #    Unnormalized slater with number n-1, this is needed to remove divide by r
+                norm = np.power(2. * exps, numbers) * np.sqrt((2. * exps) / factorial(2. * numbers))
+                slater_minus_one = SlaterAtoms.radial_slater_orbital(
+                    exps, numbers - 1, points, normalized=False
+                )
+                slater_r = norm.T * slater_minus_one
+                # When r=0 and n = 1, then slater/r is infinity.
+                i_r_zero = np.where(np.abs(points) == 0.0)[0]
+                i_numb_one = np.where(numbers[0] == 1)[0]
+                indices = np.array([[x, y] for x in i_r_zero for y in i_numb_one])
+                if len(indices) != 0:  # if-statement needed to remove numpy warning using list
+                    slater_r[indices] = np.inf
+                phi_i_r[:, index] += np.ravel(np.dot(slater_r, self.orbitals_coeff[orbital]))
+
+            factor = 2.0 * phi_i_r * molecular_orb_deriv
         factor = (
             molecular_orb_deriv**2.0 + factor + molecular_orb * molecular_orb_sec_deriv
         )
-        laplacian = np.dot(2. * factor, self.orbitals_occupation).ravel() / (4 * np.pi)
+        # Multiply by occupation number to get molecular orbitals.
+        laplacian = np.dot(2.0 * factor, self.orbitals_occupation).ravel() / (4 * np.pi)
         return laplacian
diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index acc5818..f7d4521 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -424,28 +424,31 @@ def test_kinetic_energy_cation_anion_c():
     assert_almost_equal(integral, c.kinetic_energy, decimal=5)
 
 
+def test_derivative_electron_density_h():
+    r"""Test derivative of atomic density of hydrogen."""
+    h = SlaterAtoms("h")
+    pts = np.random.uniform(0, 1, size=(100,))
+    actual = h.derivative_density(pts)
+    finite_diff = (h.atomic_density(pts + 1e-8) - h.atomic_density(pts)) / 1e-8
+    assert_allclose(actual, finite_diff)
+
+
 def test_derivative_electron_density_c():
     r"""Test derivative of atomic density of cation of carbon."""
     c = SlaterAtoms("c", cation=True)
-    eps = 1e-10
-    grid = np.array([0.1, 0.1 + eps, 0.5, 0.5 + eps])
-    dens = c.atomic_density(grid)
-    actual = c.derivative_density(np.array([0.1, 0.5]))
-    desired_0 = (dens[1] - dens[0]) / eps
-    desired_1 = (dens[3] - dens[2]) / eps
-    assert_almost_equal(actual, np.array([desired_0, desired_1]), decimal=4)
+    pts = np.random.uniform(0, 2, size=(100,))
+    actual = c.derivative_density(pts)
+    finite_diff = (c.atomic_density(pts + 1e-8) - c.atomic_density(pts)) / 1e-8
+    assert_allclose(actual, finite_diff, atol=1e-4)
 
 
 def test_derivative_electron_density_cr():
     r"""Test derivative of atomic density of chromium."""
     cr = SlaterAtoms("cr")
-    eps = 2.0e-8
-    grid = np.array([0.1 - eps, 0.1, 0.1 + eps, 1. - eps, 1., 1. + eps])
-    dens = cr.atomic_density(grid)
-    actual = cr.derivative_density(np.array([0.1, 1.]))
-    desired_0 = (dens[2] - dens[0]) / (2. * eps)
-    desired_1 = (dens[5] - dens[3]) / (2. * eps)
-    assert_almost_equal(actual, np.array([desired_0, desired_1]), decimal=4)
+    pts = np.random.uniform(0, 5, size=(100,))
+    actual = cr.derivative_density(pts)
+    finite_diff = (cr.atomic_density(pts + 1e-8) - cr.atomic_density(pts)) / 1e-8
+    assert_allclose(actual, finite_diff, rtol=1e-6, atol=1e-4)
 
 
 def test_kinetic_energy_heavy_element_ce():
@@ -462,13 +465,13 @@ def test_laplacian_of_hydrogen_slater():
     # Slater is N e^(-r)
     h = SlaterAtoms("h")
 
-    pts = np.arange(0.01, 10., 0.1)
+    pts = np.arange(0.0, 10., 0.1)
     lap = h.laplacian_of_atomic_density(pts)
 
     # Analytic derive it
     solution = np.exp(-2.0 * pts) * (pts - 1) * 16.0 / pts
     solution /= (4.0 * np.pi)
-    assert_almost_equal(lap, solution, decimal=4)
+    assert_almost_equal(lap, solution, decimal=7)
 
 
 def test_raises():

From 353def2c47df36e0229d3bfc3d2e21ff453f7550 Mon Sep 17 00:00:00 2001
From: Ali-Tehrani <alirezatehrani24@gmail.com>
Date: Mon, 6 Jun 2022 15:54:51 -0400
Subject: [PATCH 09/11] Improve documentation

---
 bfit/density.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index 2e2e519..d2ba445 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -228,7 +228,7 @@ def radial_slater_orbital(exponent, number, points, normalized=True):
         points : ndarray, (N,)
             The radial :math:`r` grid points.
         normalized : bool
-            If true, returns the normalization constant :math:`N`.
+            If true, adds the normalization constant :math:`N`.
 
         Returns
         -------
@@ -463,8 +463,7 @@ def second_derivative_radial_slater_type_orbital(exponent, number, points):
 
         Notes
         -----
-        - When :math:`n=1,2` and :math:`r=0`, then the derivative is un-defined and this
-          function returns 0 instead.
+        - When :math:`n=1,2` and :math:`r=0`, then the derivative is infinity.
 
         References
         ----------
@@ -521,7 +520,7 @@ def positive_definite_kinetic_energy(self, points):
 
         Notes
         -----
-        - When :math:`n=1` and :math:`r=0`, then kinetic energy is undefined and nan is returned.
+        - When :math:`n=1` and :math:`r=0`, then kinetic energy is either nan or infinity.
 
         """
         phi_matrix = np.zeros((len(points), len(self.orbitals)))
@@ -542,7 +541,7 @@ def positive_definite_kinetic_energy(self, points):
                 exps, numbers - 1, points, normalized=False
             )
             deriv_pref = norm.T * slater_minus_one
-            # When r=0 and n = 1, then the derivative is undefined and this returns infinity.
+            # When r=0 and n = 1, then the derivative is undefined and this returns infinity or nan
             i_r_zero = np.where(np.abs(points) == 0.0)[0]
             i_numb_one = np.where(numbers[0] == 1)[0]
             indices = np.array([[x, y] for x in i_r_zero for y in i_numb_one])

From 0929f4a11ea0ad344ae57083330b2002ee14e5da Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Tue, 7 Jun 2022 13:45:46 -0400
Subject: [PATCH 10/11] Add comments for computing 1st & 2nd deriv for r

---
 bfit/density.py | 64 +++++++++++++++++++++++++++----------------------
 1 file changed, 36 insertions(+), 28 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index d2ba445..e2f64ed 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -378,7 +378,7 @@ def atomic_density(self, points, mode="total"):
     @staticmethod
     def derivative_radial_slater_type_orbital(exponent, number, points):
         r"""
-        Compute the derivative of the radial component of Slater-type orbital on the given points.
+        Compute the first derivative of the radial component of Slater-type orbital.
 
         The derivative of the Slater-type orbital is defined as:
 
@@ -406,23 +406,26 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
         slater : ndarray, (N, M)
             First derivative of Slater-type orbitals evaluated on the :math:`N` grid points.
 
-        References
-        ----------
-        See wikipedia page on "Slater-Type orbitals".
-
         """
         norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number))
+        # compute R(r) = N r^{n - 1} e^{-\alpha r}
         slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True)
-        # Calculate the un-normalized Slater with number less than one
-        slater_minus = SlaterAtoms.radial_slater_orbital(
-            exponent, number - 1, points, normalized=False
-        )
-        # Consider the case when dividing by zero.
+
+        # compute N (-\alpha) e^{-\alpha r} part of derivative which exists for for all n
+        # -------------------------------------------------------------------------------
+        # n=1: R(r) = N e^{-\alpha r}, so dR(r)/dr = N (-\alpha) e^{-\alpha r}
         deriv = np.zeros((len(points), number.shape[0]))
-        # Compute -\alpha * slater
         deriv -= exponent.T * slater
-        # Compute (n-1)*slater/r, note that this only occurs when n!=1.
-        deriv_pref = norm.T * slater_minus * (number.T - 1.0)
+
+        # compute part of the derivative which only exists for n != 1
+        # -----------------------------------------------------------
+        # n!=1: dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r}
+        # calculate the un-normalized Slater with n - 1; i.e., r^{n - 2} e^{-\alpha r}
+        slater_minus_one = SlaterAtoms.radial_slater_orbital(
+            exponent, number - 1, points, normalized=False
+        )
+        # compute N (n - 1) r^{n - 2} e^{-\alpha r} for n != 1
+        deriv_pref = norm.T * slater_minus_one * (number.T - 1.0)
         i_numb_one = np.where(number == 1)[0]
         deriv_pref[:, i_numb_one] = 0.0
         deriv += deriv_pref
@@ -431,7 +434,7 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
     @staticmethod
     def second_derivative_radial_slater_type_orbital(exponent, number, points):
         r"""
-        Compute the derivative of the radial component of Slater-type orbital on the given points.
+        Compute the second derivative of the radial component of Slater-type orbital.
 
         The derivative of the Slater-type orbital is defined as:
 
@@ -461,33 +464,38 @@ def second_derivative_radial_slater_type_orbital(exponent, number, points):
         slater : ndarray, (N, M)
             Second derivative of Slater-type orbitals evaluated on the :math:`N` grid points.
 
-        Notes
-        -----
-        - When :math:`n=1,2` and :math:`r=0`, then the derivative is infinity.
-
-        References
-        ----------
-        See wikipedia page on "Slater-Type orbitals".
-
         """
         norm = np.power(2. * exponent, number) * np.sqrt((2. * exponent) / factorial(2. * number))
+        # compute R(r) = N r^{n - 1} e^{-\alpha r}
         slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True)
-        # Calculate the un-normalized Slater with number less than one
+        # calculate the un-normalized Slater with n - 1 and n - 2
         with np.errstate(divide='ignore'):
-            slater_minus = SlaterAtoms.radial_slater_orbital(
+            slater_minus_one = SlaterAtoms.radial_slater_orbital(
                 exponent, number - 1, points, normalized=False
             )
             slater_minus_two = SlaterAtoms.radial_slater_orbital(
                 exponent, number - 2, points, normalized=False
             )
-        # Compute \alpha^2 * slater
+
+        # General formula for the first derivative:
+        # dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r}
+
+        # compute N \alpha^2 e^{-\alpha r} part of derivative which exists for for all n
+        # ------------------------------------------------------------------------------
+        # when n=1, R(r) = N e^{-\alpha r}, so d^2R(r)/dr^2 = N \alpha^2 e^{-\alpha r}
         deriv = exponent.T**2.0 * slater
-        # Compute 2 \alpha (n-1) * slater / r, note doesn't occur when n=1
-        deriv_pref = 2.0 * exponent.T * norm.T * slater_minus * (number.T - 1.0)
+
+        # compute part of the derivative which only exists for n != 1
+        # -----------------------------------------------------------
+        # calculate 2 * N (n - 1) (-\alpha) r^{n - 2} e^{-\alpha r}
+        deriv_pref = 2.0 * exponent.T * norm.T * slater_minus_one * (number.T - 1.0)
         i_numb_one = np.where(number == 1)[0]
         deriv_pref[:, i_numb_one] = 0.0
         deriv -= deriv_pref
-        # Compute (n-1)(n-2) * slater / r^2, note doesn't occur when n=1,2
+
+        # compute part of the derivative which only exists for n != 1 & n != 2
+        # --------------------------------------------------------------------
+        # calculate N (n - 1) (n - 2) r^{n - 3} e^{-\alpha r}
         deriv_pref = norm.T * slater_minus_two * (number.T - 1.0) * (number.T - 2.0)
         i_numb_one = np.where(number == 1)[0]
         deriv_pref[:, i_numb_one] = 0.0

From 09c4f9b2d4fa0008421173d40ec4b7be295203ec Mon Sep 17 00:00:00 2001
From: Farnaz Heidar-Zadeh <farnaz_chem@yahoo.com>
Date: Tue, 7 Jun 2022 14:21:04 -0400
Subject: [PATCH 11/11] Rename method for computing 1st derivative

---
 bfit/density.py           | 12 +++++++-----
 bfit/test/test_density.py | 14 +++++++-------
 2 files changed, 14 insertions(+), 12 deletions(-)

diff --git a/bfit/density.py b/bfit/density.py
index e2f64ed..04b0c05 100644
--- a/bfit/density.py
+++ b/bfit/density.py
@@ -307,7 +307,7 @@ def phi_matrix(self, points, deriv=None):
         for index, orbital in enumerate(self.orbitals):
             exps, number = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]]
             if deriv == 1:
-                slater = self.derivative_radial_slater_type_orbital(exps, number, points)
+                slater = self.first_derivative_radial_slater_type_orbital(exps, number, points)
             elif deriv == 2:
                 slater = self.second_derivative_radial_slater_type_orbital(exps, number, points)
             else:
@@ -376,7 +376,7 @@ def atomic_density(self, points, mode="total"):
         return dens
 
     @staticmethod
-    def derivative_radial_slater_type_orbital(exponent, number, points):
+    def first_derivative_radial_slater_type_orbital(exponent, number, points):
         r"""
         Compute the first derivative of the radial component of Slater-type orbital.
 
@@ -411,15 +411,17 @@ def derivative_radial_slater_type_orbital(exponent, number, points):
         # compute R(r) = N r^{n - 1} e^{-\alpha r}
         slater = SlaterAtoms.radial_slater_orbital(exponent, number, points, normalized=True)
 
+        # Derivative of R(r) w.r.t. r:
+        # n=1: R(r) = N e^{-\alpha r}, so dR(r)/dr = N (-\alpha) e^{-\alpha r} = -\alpha R(r)
+        # n!=1: dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r}
+
         # compute N (-\alpha) e^{-\alpha r} part of derivative which exists for for all n
         # -------------------------------------------------------------------------------
-        # n=1: R(r) = N e^{-\alpha r}, so dR(r)/dr = N (-\alpha) e^{-\alpha r}
         deriv = np.zeros((len(points), number.shape[0]))
         deriv -= exponent.T * slater
 
         # compute part of the derivative which only exists for n != 1
         # -----------------------------------------------------------
-        # n!=1: dR(r)/dr = N (-\alpha) r^{n - 1} e^{-\alpha r} + N (n - 1) r^{n - 2} e^{-\alpha r}
         # calculate the un-normalized Slater with n - 1; i.e., r^{n - 2} e^{-\alpha r}
         slater_minus_one = SlaterAtoms.radial_slater_orbital(
             exponent, number - 1, points, normalized=False
@@ -539,7 +541,7 @@ def positive_definite_kinetic_energy(self, points):
             exps, numbers = self.orbitals_exp[orbital[1]], self.basis_numbers[orbital[1]]
             # Calculate del^2 of radial component
             # derivative of the radial component
-            deriv_radial = self.derivative_radial_slater_type_orbital(exps, numbers, points)
+            deriv_radial = self.first_derivative_radial_slater_type_orbital(exps, numbers, points)
             phi_matrix[:, index] = np.ravel(np.dot(deriv_radial, self.orbitals_coeff[orbital])**2.0)
 
             # Calculate del^2 of spherical component
diff --git a/bfit/test/test_density.py b/bfit/test/test_density.py
index f7d4521..6790216 100644
--- a/bfit/test/test_density.py
+++ b/bfit/test/test_density.py
@@ -70,21 +70,21 @@ def test_derivative_slater_type_orbital_be():
     # load Be atomic wave function
     be = SlaterAtoms("Be")
     # check values of a single orbital at r=1.0
-    orbital = be.derivative_radial_slater_type_orbital(np.array([[12.683501]]),
-                                                       np.array([[1]]), np.array([1]))
+    orbital = be.first_derivative_radial_slater_type_orbital(np.array([[12.683501]]),
+                                                             np.array([[1]]), np.array([1]))
     assert_almost_equal(orbital, slater(12.683501, 1, 1.0, derivative=True), decimal=6)
     # check values of a single orbital at r=2.0
-    orbital = be.derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]),
-                                                       np.array([2]))
+    orbital = be.first_derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]),
+                                                             np.array([2]))
     assert_almost_equal(orbital, slater(0.821620, 2, 2.0, derivative=True), decimal=6)
     # check value of single orbital at r = 0.0
-    orbital = be.derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]),
-                                                       np.array([0.]))
+    orbital = be.first_derivative_radial_slater_type_orbital(np.array([[0.821620]]), np.array([[2]]),
+                                                             np.array([0.]))
     actual = np.power(2. * 0.821620, 2) * np.sqrt((2. * 0.821620) / (4 * 3 * 2))
     assert_almost_equal(orbital, actual, decimal=6)
     # check value of tow orbitals at r=1.0 & r=2.0
     exps, nums = np.array([[12.683501], [0.821620]]), np.array([[1], [2]])
-    orbitals = be.derivative_radial_slater_type_orbital(exps, nums, np.array([1., 2.]))
+    orbitals = be.first_derivative_radial_slater_type_orbital(exps, nums, np.array([1., 2.]))
     expected = np.array([[slater(exps[0, 0], nums[0, 0], 1., True),
                           slater(exps[1, 0], nums[1, 0], 1., True)],
                          [slater(exps[0, 0], nums[0, 0], 2., True),