Merge pull request #131 from theochem/add_laplacian

Laplacian method for Species class

atomdb/datasets/numeric/run.py

@@ -31,13 +31,28 @@
 
 import atomdb
 
-from atomdb.utils import MULTIPLICITIES
+from atomdb.utils import MULTIPLICITIES, DEFAULT_DATAPATH
 
 from atomdb.periodic import Element
 
 
-def load_numerical_hf_data():
-    """Load data from desnity.out file into a `SpeciesTable`."""
+def load_numerical_hf_data(data_path):
+    """Load data from desnity.out file into a `SpeciesTable`.
+    
+    Parameters
+    ----------
+    data_path : str
+        Path to the directory containing a folder named `raw` where the desnity.out file is stored.
+    
+    Returns
+    -------
+    species : dict
+        Dictionary of atomic species containing the information from the numeric Hartree-Fock calculation.
+        This is energy components, grid, density, gradient, and laplacian values.
+    
+    """
+    # set the path to the raw data
+    data_path = os.path.join(data_path, "numeric", "raw")
 
     from io import StringIO
 
@@ -56,7 +71,7 @@ def helper_data():
         return data[:, 0], data[:, 1], data[:, 2], data[:, 3]
 
     species = {}
-    with open(os.path.join(os.path.dirname(__file__), "raw/density.out"), "r") as f:
+    with open(os.path.join(data_path, "density.out"), "r") as f:
         line = f.readline()
         while line:
             if line.startswith(" 1st line is atomic no"):
@@ -99,6 +114,48 @@ def helper_data():
     return species
 
 
+def eval_radial_dd_density(gradient, laplacian, points, err='ignore', tol=1e-10):
+    """Helper function to compute the radial second derivative of the density.
+
+    From a set of radial points :math:`r`, the gradient of the density, :math:`df/dr`, and the
+    Laplacian of the density, :math:`\nabla^2 f`, the radial second derivative of the density is
+    computed as:
+    
+    .. math::
+        d/dr (df/dr) = \nabla^2 f - 2/r * df/dr
+    
+    Parameters
+    ----------
+    gradient : np.ndarray
+        Gradient of the density.
+    laplacian : np.ndarray
+        Laplacian of the density.
+    points : np.ndarray
+        Radial points where the density gradient and Laplacian are evaluated.
+    err : str, optional
+        Error handling for division by zero.
+    tol : float, optional
+        Tolerance for the points close to zero.
+    
+    Returns
+    -------
+    d2dens : np.ndarray
+        Radial second derivative of the density.
+    
+    Notes
+    -----
+    When the points are close to zero, the radial second derivative of the density tends to infinity.
+    In this case, this function returns zero.
+
+    """
+    # Handle the case when the points are close to zero
+    with np.errstate(divide=err):
+        # Compute the radial second derivative of the density
+        d2dens =  laplacian - 2 * gradient / points
+        d2dens = np.where(points < tol, 0.0, d2dens)
+    return d2dens
+
+
 DOCSTRING = """Numeric Dataset
 
 Load data from desnity.out file into a `SpeciesTable`.
@@ -134,7 +191,7 @@ def run(elem, charge, mult, nexc, dataset, datapath):
             f"Expected multiplicity is {expected_mult}."
         )
 
-    species_table = load_numerical_hf_data()
+    species_table = load_numerical_hf_data(datapath)
     data = species_table[(atnum, nelec)]
 
     #
@@ -165,6 +222,9 @@ def run(elem, charge, mult, nexc, dataset, datapath):
     lapl_tot = data["laplacian"]
     ked_tot = None
 
+    # Compute the second derivative of the density
+    dd_dens_tot = eval_radial_dd_density(d_dens_tot, lapl_tot, points)
+
     # Return Species instance
     fields = dict(
         elem=elem,
@@ -183,7 +243,7 @@ def run(elem, charge, mult, nexc, dataset, datapath):
         rs=points,
         dens_tot=dens_tot,
         d_dens_tot=d_dens_tot,
-        dd_dens_tot=lapl_tot,
+        dd_dens_tot=dd_dens_tot,
         ked_tot=ked_tot,
     )
     return atomdb.Species(dataset, fields)

atomdb/species.py

@@ -666,7 +666,7 @@ def d_dens_func(self):
     @spline
     def dd_dens_func(self):
         r"""
-        Return a cubic spline of the electronic density Laplacian.
+        Return a cubic spline of the second derivative of the electronic density.
 
         Parameters
         ----------
@@ -683,14 +683,60 @@ def dd_dens_func(self):
 
         Returns
         -------
-        DensitySpline
-            A DensitySpline instance for the density and its derivatives.
-            Given a set of radial points, it can evaluate densities and
-            derivatives up to order 2.
+        Callable[[np.ndarray(N,), int] -> np.ndarray(N,)]
+            a callable function evaluating the second derivative of the density given a set of radial
+            points (1-D array).
 
         """
         pass
 
+    def dd_dens_lapl_func(self, spin="t", index=None, log=False):
+        r"""
+        Return the function for the electronic density Laplacian.
+
+        .. math::
+            
+            \nabla^2 \rho(\mathbf{r}) = \frac{d^2 \rho(r)}{dr^2} + \frac{2}{r} \frac{d \rho(r)}{dr}
+
+        Parameters
+        ----------
+        spin : str, default="t"
+            Type of occupied spin orbitals.
+            Can be either "t" (for alpha + beta), "a" (for alpha),
+            "b" (for beta), or "m" (for alpha - beta).
+        index : sequence of int, optional
+            Sequence of integers representing the spin orbitals.
+            These are indexed from 0 to the number of basis functions.
+            By default, all orbitals of the given spin(s) are included.
+        log : bool, default=False
+            Whether the logarithm of the density is used for interpolation.
+        
+        Returns
+        -------
+        Callable[np.ndarray(N,) -> np.ndarray(N,)]
+            a callable function evaluating the Laplacian of the density given a set of radial
+            points (1-D array).
+        
+        Notes
+        -----
+        When this function is evaluated at a point close to zero, the Laplacian becomes undefined.
+        In this case, this function returns zero.
+
+        """
+        # Obtain cubic spline functions for the first and second derivatives of the density
+        d_dens_sp_spline = self.d_dens_func(spin=spin, index=index, log=log)
+        dd_dens_spline = self.dd_dens_func(spin=spin, index=index, log=log)
+
+        # Define the Laplacian function
+        def densityspline_like_func(rs):
+            # Avoid division by zero and handle small values of r
+            with np.errstate(divide='ignore'):
+                laplacian = dd_dens_spline(rs) + 2 * d_dens_sp_spline(rs) / rs
+                laplacian = np.where(rs < 1e-10, 0.0, laplacian)
+            return laplacian
+
+        return densityspline_like_func
+
     @spline
     def ked_func(self):
         r"""

atomdb/test/test_numeric.py

@@ -64,18 +64,13 @@ def test_numerical_hf_data_h():
     assert np.allclose(sp._data.dens_tot[-2:], [0.0, 0.0], atol=1e-10)
 
     # evaluate radial density gradient (first derivative of density spline)
-    dens = sp.dens_func(spin="t", log=True)
-    gradient = dens(sp._data.rs, deriv=1)
+    gradient = sp.d_dens_func(log=False)(sp._data.rs)
 
-    # load gradient reference values from numerical HF raw files
-    fname = "001_q000_m02_numeric_gradient.npy"
-    ref_grad = np.load(f"{TEST_DATAPATH}/numeric/db/{fname}")
-
-    # check interpolated gradient values against reference
-    # close to the nuclei the spline derivative does not describe the gradient well
-    assert np.allclose(gradient[:2], ref_grad[:2], atol=1e-3)
-    # away from the nuclei the spline derivative describes the gradient well
-    assert np.allclose(gradient[-2:], ref_grad[-2:], atol=1e-10)
+    # check interpolated gradient values against reference values from numerical HF raw files
+    # close to the nuclei
+    assert np.allclose(gradient[:2], [-0.636619761671399, -0.613581739284137], atol=1e-10)
+    # away from the nuclei
+    assert np.allclose(gradient[-2:], [0.0, 0.0], atol=1e-10)
 
 
 def test_numerical_hf_data_h_anion():
@@ -116,18 +111,12 @@ def test_numerical_hf_data_h_anion():
     assert np.allclose(sp._data.dens_tot[-20:], 0.0, atol=1e-10)
 
     # evaluate radial density gradient (first derivative of density spline)
-    dens = sp.dens_func(spin="t", log=True)
-    gradient = dens(sp._data.rs, deriv=1)
-
-    # load gradient reference values from numerical HF raw files
-    fname = "001_q-01_m01_numeric_gradient.npy"
-    ref_grad = np.load(f"{TEST_DATAPATH}/numeric/db/{fname}")
+    gradient = sp.d_dens_func(log=False)(sp._data.rs)
 
-    # check interpolated gradient values against reference
-    # close to the nuclei the spline derivative does not describe the gradient well
-    assert np.allclose(gradient[:2], ref_grad[:2], atol=1e-3)
-    # away from the nuclei the spline derivative describes the gradient well
-    assert np.allclose(gradient[-2:], ref_grad[-2:], atol=1e-10)
+    # check interpolated gradient values against reference values from numerical HF raw files
+    assert np.allclose(gradient[:2], [-0.618386750431843, -0.594311093621533], atol=1e-10)
+    assert np.allclose(gradient[20:22], [-0.543476018733641, -0.538979599233911], atol=1e-10)
+    assert np.allclose(gradient[-20:], [0.0] * 20, atol=1e-10)
 
 
 @pytest.mark.parametrize(
@@ -163,7 +152,11 @@ def test_numerical_hf_atomic_density(atom, mult, npoints, nelec):
     assert all(sp._data.dens_tot >= 0.0)
 
     # check the density integrates to the correct number of electrons
-    assert_almost_equal(4 * np.pi * np.trapz(grid**2 * dens, grid), nelec, decimal=2)
+    if hasattr(np, "trapezoid"):
+        int_density = 4.0 * np.pi * np.trapezoid(grid**2 * dens, grid)
+    else:
+        int_density = 4.0 * np.pi * np.trapz(grid**2 * dens, grid)
+    assert_almost_equal(int_density, nelec, decimal=2)
 
     # get density spline and check its values
     spline = sp.dens_func(spin="t", log=True)
@@ -195,26 +188,48 @@ def test_numerical_hf_density_gradient(atom, charge, mult):
     assert np.allclose(gradient[indx_radii], ref_grad[indx_radii], atol=1e-3)
 
 
-@pytest.mark.xfail(reason="High errors in spline derivative of order 2")
+@pytest.mark.xfail(reason="High errors in spline derivative of order 2 at intermediate distances")
 @pytest.mark.parametrize(
-    "atom, charge, mult", [("H", 0, 2), ("H", -1, 1), ("Be", 0, 1), ("Cl", 0, 3), ("Ne", 0, 1)]
+    "atom, charge, mult", [("H", 0, 2), ("H", -1, 1), ("Be", 0, 1), ("Cl", 0, 2), ("Ne", 0, 1)]
 )
-def test_numerical_hf_density_laplacian(atom, charge, mult):
+def test_numerical_hf_dd_density(atom, charge, mult):
     # load atomic and density data
     sp = load(atom, charge, mult, dataset="numeric", datapath=TEST_DATAPATH)
 
-    # evaluate density and laplacian (second derivative of density spline) on the radial grid
-    grid = sp._data.rs
-    spline = sp.dens_func(spin="t", log=False)
-    lapl = spline(grid, deriv=2)
+    # evaluate the second derivative of the density on the radial grid
+    dd_dens = sp.dd_dens_func(log=False)(sp._data.rs)
 
     # check shape of arrays
-    assert lapl.shape == grid.shape
+    assert dd_dens.shape == sp._data.rs.shape
+
+    # check interpolated density derivative values against reference values
+    # far away from the nuclei, the second derivative of the density is close to zero
+    assert np.allclose(dd_dens[-10:], [0.0] * 10, atol=1e-10)
+    # for r=0, the second derivative of the density is set to zero
+    assert np.allclose(dd_dens[0], [0.0], atol=1e-10)
+
+    # WARNING: The values of the second order derivative of the density at intermediate r distances
+    # are not tested. Comparisong agains deriv=2 of the density spline:
+    # ref_dd_dens = sp.dens_func(log=True)(sp._data.rs, deriv=2)
+    # results in high errors, rendering this test case as unreliable.
+
+
+@pytest.mark.parametrize(
+    "atom, charge, mult", [("H", 0, 2), ("H", -1, 1), ("Be", 0, 1), ("Cl", 0, 2), ("Ne", 0, 1)]
+)
+def test_numerical_hf_density_laplacian(atom, charge, mult):
+    # load atomic and density data
+    sp = load(atom, charge, mult, dataset="numeric", datapath=TEST_DATAPATH)
+
+    # evaluate the Laplacian of the density on the radial grid
+    laplacian_dens = sp.dd_dens_lapl_func(log=False)(sp._data.rs)
 
     # load reference values from numerical HF raw files
     id = f"{str(sp.atnum).zfill(3)}_q{str(charge).zfill(3)}_m{mult:02d}"
     fname = f"{id}_numeric_laplacian.npy"
     ref_lapl = np.load(f"{TEST_DATAPATH}/numeric/db/{fname}")
 
-    # interpolated laplacian values against reference
-    assert np.allclose(lapl, ref_lapl, atol=1e-6)
+    # check interpolated Laplacian of density values against reference values
+    assert np.allclose(laplacian_dens, ref_lapl, atol=1e-10)
+    # for r=0, the Laplacian function in not well defined and is set to zero
+    assert np.allclose(laplacian_dens[0], [0.0], atol=1e-10)