diff --git a/atomdb/datasets/numeric/run.py b/atomdb/datasets/numeric/run.py index f1ac7354..c24ef5f7 100644 --- a/atomdb/datasets/numeric/run.py +++ b/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) diff --git a/atomdb/species.py b/atomdb/species.py index c74f3692..3ace0942 100644 --- a/atomdb/species.py +++ b/atomdb/species.py @@ -658,7 +658,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 ---------- @@ -675,14 +675,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""" diff --git a/atomdb/test/data/numeric/db/Be_000_001_000.msg b/atomdb/test/data/numeric/db/Be_000_001_000.msg index 51b553be..63172c5b 100644 Binary files a/atomdb/test/data/numeric/db/Be_000_001_000.msg and b/atomdb/test/data/numeric/db/Be_000_001_000.msg differ diff --git a/atomdb/test/data/numeric/db/Cl_000_002_000.msg b/atomdb/test/data/numeric/db/Cl_000_002_000.msg index 086baf7b..ff6fea77 100644 Binary files a/atomdb/test/data/numeric/db/Cl_000_002_000.msg and b/atomdb/test/data/numeric/db/Cl_000_002_000.msg differ diff --git a/atomdb/test/data/numeric/db/H_-01_001_000.msg b/atomdb/test/data/numeric/db/H_-01_001_000.msg index b61afa8b..b65a27da 100644 Binary files a/atomdb/test/data/numeric/db/H_-01_001_000.msg and b/atomdb/test/data/numeric/db/H_-01_001_000.msg differ diff --git a/atomdb/test/data/numeric/db/H_000_002_000.msg b/atomdb/test/data/numeric/db/H_000_002_000.msg index 308f63ab..63e18740 100644 Binary files a/atomdb/test/data/numeric/db/H_000_002_000.msg and b/atomdb/test/data/numeric/db/H_000_002_000.msg differ diff --git a/atomdb/test/data/numeric/db/Ne_000_001_000.msg b/atomdb/test/data/numeric/db/Ne_000_001_000.msg index 4973b63a..e0acef9d 100644 Binary files a/atomdb/test/data/numeric/db/Ne_000_001_000.msg and b/atomdb/test/data/numeric/db/Ne_000_001_000.msg differ diff --git a/atomdb/test/test_numeric.py b/atomdb/test/test_numeric.py index e25debda..ee3d73ab 100644 --- a/atomdb/test/test_numeric.py +++ b/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)