From c9d10db4aa2038377b766c31ebe709871f7a629e Mon Sep 17 00:00:00 2001
From: Philip Loche <ploche@physik.fu-berlin.de>
Date: Fri, 12 Apr 2024 14:04:58 +0200
Subject: [PATCH] Add radial integral for custom densities
---
.../rascaline/utils/splines/atomic_density.py | 116 ++++++++++++++----
.../rascaline/utils/splines/splines.py | 84 ++++++++++---
.../rascaline/tests/utils/atomic_density.py | 18 +++
python/rascaline/tests/utils/splines.py | 45 +++++++
4 files changed, 222 insertions(+), 41 deletions(-)
create mode 100644 python/rascaline/tests/utils/atomic_density.py
diff --git a/python/rascaline/rascaline/utils/splines/atomic_density.py b/python/rascaline/rascaline/utils/splines/atomic_density.py
index 610309041..2fa446798 100644
--- a/python/rascaline/rascaline/utils/splines/atomic_density.py
+++ b/python/rascaline/rascaline/utils/splines/atomic_density.py
@@ -52,12 +52,22 @@ class AtomicDensityBase(ABC):
@abstractmethod
def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
- """Compute the atomic density arising from atoms at ``positions``
+ """Compute the atomic density arising from atoms at ``positions``.
:param positions: positions to evaluate the atomic densities
:returns: evaluated atomic density
"""
+ @abstractmethod
+ def compute_derivative(
+ self, positions: Union[float, np.ndarray]
+ ) -> Union[float, np.ndarray]:
+ """Derivative of the atomic density arising from atoms at ``positions``.
+
+ :param positions: positions to evaluate the derivatives atomic densities
+ :returns: evaluated derivative of the atomic density with respect to positions
+ """
+
class DeltaDensity(AtomicDensityBase):
r"""Delta atomic densities of the form :math:`g(r)=\delta(r)`."""
@@ -67,6 +77,14 @@ def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarra
"Compute function of the delta density should never called directly."
)
+ def compute_derivative(
+ self, positions: Union[float, np.ndarray]
+ ) -> Union[float, np.ndarray]:
+ raise ValueError(
+ "Compute derivative function of the delta density should never called "
+ "directly."
+ )
+
class GaussianDensity(AtomicDensityBase):
r"""Gaussian atomic density function.
@@ -83,6 +101,14 @@ class GaussianDensity(AtomicDensityBase):
\|g\|^2 = \int \mathrm{d}^3\boldsymbol{r} |g(r)|^2 = 1\,,
+ The derivatives of the Gaussian atomic density with respect to the position is
+
+ .. math::
+
+ g^\prime(r) =
+ \frac{\partial g(r)}{\partial r} = \frac{-r}{\sigma^2(\pi
+ \sigma^2)^{3/4}}e^{-\frac{r^2}{2\sigma^2}} \,.
+
:param atomic_gaussian_width: Width of the atom-centered gaussian used to create the
atomic density
"""
@@ -90,11 +116,26 @@ class GaussianDensity(AtomicDensityBase):
def __init__(self, atomic_gaussian_width: float):
self.atomic_gaussian_width = atomic_gaussian_width
- def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
+ def _compute(
+ self, positions: Union[float, np.ndarray], derivative: bool = False
+ ) -> Union[float, np.ndarray]:
atomic_gaussian_width_sq = self.atomic_gaussian_width**2
- return np.exp(-0.5 * positions**2 / atomic_gaussian_width_sq) / (
- np.pi * atomic_gaussian_width_sq
- ) ** (3 / 4)
+ x = positions**2 / (2 * atomic_gaussian_width_sq)
+
+ density = np.exp(-x) / (np.pi * atomic_gaussian_width_sq) ** (3 / 4)
+
+ if derivative:
+ density *= -positions / atomic_gaussian_width_sq
+
+ return density
+
+ def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
+ return self._compute(positions=positions, derivative=False)
+
+ def compute_derivative(
+ self, positions: Union[float, np.ndarray]
+ ) -> Union[float, np.ndarray]:
+ return self._compute(positions=positions, derivative=True)
class LodeDensity(AtomicDensityBase):
@@ -108,10 +149,11 @@ class LodeDensity(AtomicDensityBase):
\frac{\gamma\left( \frac{p}{2}, \frac{r^2}{2\sigma^2} \right)}
{r^p},
- where :math:`\Gamma(z)` is the Gamma function and :math:`\gamma(a, x)` is the
- incomplete lower Gamma function. However its evaluation at :math:`r=0` is
- problematic because :math:`g(r)` is of the form :math:`0/0`. For practical
- implementations, it is thus more convenient to rewrite the density as
+ where :math:`p` is the potential exponent, :math:`\Gamma(z)` is the Gamma function
+ and :math:`\gamma(a, x)` is the incomplete lower Gamma function. However its
+ evaluation at :math:`r=0` is problematic because :math:`g(r)` is of the form
+ :math:`0/0`. For practical implementations, it is thus more convenient to rewrite
+ the density as
.. math::
@@ -123,10 +165,10 @@ class LodeDensity(AtomicDensityBase):
& x \geq 10^{-5}
\end{cases}
- It is convenient to use the expression for sufficiently small :math:`x` since the
- relative weight of the first neglected term is on the order of :math:`1/6x^3`.
- Therefore, the threshold :math:`x = 10^{-5}` leads to relative errors on the order
- of the machine epsilon.
+ where :math:`a=p/2`. It is convenient to use the expression for sufficiently small
+ :math:`x` since the relative weight of the first neglected term is on the order of
+ :math:`1/6x^3`. Therefore, the threshold :math:`x = 10^{-5}` leads to relative
+ errors on the order of the machine epsilon.
:param atomic_gaussian_width: Width of the atom-centered gaussian used to create the
atomic density
@@ -141,18 +183,32 @@ def __init__(self, atomic_gaussian_width: float, potential_exponent: int):
if not HAS_SCIPY:
raise ValueError("LodeDensity requires scipy to be installed")
- self.potential_exponent = potential_exponent
self.atomic_gaussian_width = atomic_gaussian_width
+ self.potential_exponent = potential_exponent
- def _short_range(self, a, x):
- return 1 / a - x / (a + 1) + x**2 / (2 * (a + 2))
+ def _short_range(
+ self, a: float, x: Union[float, np.ndarray], derivative: bool = False
+ ):
+ if derivative:
+ return -1 / (a + 1) + x / (a + 2)
+ else:
+ return 1 / a - x / (a + 1) + x**2 / (2 * (a + 2))
- def _long_range(self, a, x):
- return gamma(a) * gammainc(a, x) / x**a
+ def _long_range(
+ self, a: float, x: Union[float, np.ndarray], derivative: bool = False
+ ):
+ if derivative:
+ return (np.exp(-x) - a * gamma(a) * gammainc(a, x) / x**a) / x
+ else:
+ return gamma(a) * gammainc(a, x) / x**a
- def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
+ def _compute(
+ self, positions: Union[float, np.ndarray], derivative: bool = False
+ ) -> Union[float, np.ndarray]:
if self.potential_exponent == 0:
- return GaussianDensity.compute(self, positions=positions)
+ return GaussianDensity._compute(
+ self, positions=positions, derivative=derivative
+ )
else:
atomic_gaussian_width_sq = self.atomic_gaussian_width**2
a = self.potential_exponent / 2
@@ -160,8 +216,22 @@ def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarra
prefac = 1 / gamma(a) / (2 * atomic_gaussian_width_sq) ** a
- return prefac * np.where(
+ density = prefac * np.where(
x < 1e-5,
- self._short_range(a, x),
- self._long_range(a, x),
+ self._short_range(a, x, derivative=derivative),
+ self._long_range(a, x, derivative=derivative),
)
+
+ # add inner derivative: ∂x/∂r
+ if derivative:
+ density *= positions / atomic_gaussian_width_sq
+
+ return density
+
+ def compute(self, positions: Union[float, np.ndarray]) -> Union[float, np.ndarray]:
+ return self._compute(positions=positions, derivative=False)
+
+ def compute_derivative(
+ self, positions: Union[float, np.ndarray]
+ ) -> Union[float, np.ndarray]:
+ return self._compute(positions=positions, derivative=True)
diff --git a/python/rascaline/rascaline/utils/splines/splines.py b/python/rascaline/rascaline/utils/splines/splines.py
index 7aaf9bddf..8eb3f211b 100644
--- a/python/rascaline/rascaline/utils/splines/splines.py
+++ b/python/rascaline/rascaline/utils/splines/splines.py
@@ -43,8 +43,8 @@
try:
- from scipy.integrate import quad, quad_vec
- from scipy.special import spherical_in, spherical_jn
+ from scipy.integrate import dblquad, quad, quad_vec
+ from scipy.special import legendre, spherical_in, spherical_jn
HAS_SCIPY = True
except ImportError:
@@ -593,7 +593,7 @@ def radial_integral_derivative(
elif type(self.density) is GaussianDensity:
return self._radial_integral_gaussian_derivative(n, ell, positions)
else:
- return self._radial_integral_custom_derivative(n, ell, positions)
+ return self._radial_integral_custom(n, ell, positions, derivative=True)
def _radial_integral_delta(
self, n: int, ell: int, positions: np.ndarray
@@ -632,9 +632,10 @@ def integrand(
return (
prefactor
* quad_vec(
- f=lambda x: integrand(x, n, ell, positions),
+ f=integrand,
a=0,
b=self.basis.integration_radius,
+ args=(n, ell, positions),
)[0]
)
@@ -666,26 +667,71 @@ def integrand(
return atomic_gaussian_width_sq**-1 * (
prefactor
* quad_vec(
- f=lambda x: integrand(x, n, ell, positions),
+ f=integrand,
a=0,
b=self.basis.integration_radius,
+ args=(n, ell, positions),
)[0]
- positions * self._radial_integral_gaussian(n, ell, positions)
)
def _radial_integral_custom(
- self, n: int, ell: int, positions: np.ndarray, derivative: bool
+ self, n: int, ell: int, positions: np.ndarray, derivative: bool = False
) -> np.ndarray:
- raise NotImplementedError(
- "Radial integral with custom atomic densities is not implemented yet!"
- )
- def _radial_integral_custom_derivative(
- self, n: int, ell: int, positions: np.ndarray, derivative: bool
- ) -> np.ndarray:
- raise NotImplementedError(
- "Radial integral with custom atomic densities is not implemented yet!"
- )
+ P_ell = legendre(ell)
+
+ if derivative:
+
+ def integrand(
+ u: float, integrand_position: float, n: int, ell: int, position: float
+ ) -> float:
+ arg = np.sqrt(
+ integrand_position**2
+ + position**2
+ - 2 * integrand_position * position * u
+ )
+
+ return (
+ integrand_position**2
+ * self.basis.compute(n, ell, integrand_position)
+ * P_ell(u)
+ * (position - u * integrand_position)
+ * self.density.compute_derivative(arg)
+ / arg
+ )
+
+ else:
+
+ def integrand(
+ u: float, integrand_position: float, n: int, ell: int, position: float
+ ) -> float:
+ arg = np.sqrt(
+ integrand_position**2
+ + position**2
+ - 2 * integrand_position * position * u
+ )
+
+ return (
+ integrand_position**2
+ * self.basis.compute(n, ell, integrand_position)
+ * P_ell(u)
+ * self.density.compute(arg)
+ )
+
+ radial_integral = np.zeros(len(positions))
+
+ for i, position in enumerate(positions):
+ radial_integral[i], _ = dblquad(
+ func=integrand,
+ a=0,
+ b=self.basis.integration_radius,
+ gfun=-1,
+ hfun=1,
+ args=(n, ell, position),
+ )
+
+ return 2 * np.pi * radial_integral
class LodeSpliner(RadialIntegralSplinerBase):
@@ -799,9 +845,10 @@ def integrand(
)
return quad_vec(
- f=lambda x: integrand(x, n, ell, positions),
+ f=integrand,
a=0,
b=self.basis.integration_radius,
+ args=(n, ell, positions),
)[0]
def radial_integral_derivative(
@@ -817,9 +864,10 @@ def integrand(
)
return quad_vec(
- f=lambda x: integrand(x, n, ell, positions),
+ f=integrand,
a=0,
b=self.basis.integration_radius,
+ args=(n, ell, positions),
)[0]
@property
@@ -848,5 +896,5 @@ def integrand(integrand_position: float, n: int) -> np.ndarray:
func=integrand,
a=0,
b=self.basis.integration_radius,
- args=(n),
+ args=(n,),
)[0]
diff --git a/python/rascaline/tests/utils/atomic_density.py b/python/rascaline/tests/utils/atomic_density.py
new file mode 100644
index 000000000..ae2c9f95c
--- /dev/null
+++ b/python/rascaline/tests/utils/atomic_density.py
@@ -0,0 +1,18 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from rascaline.utils import GaussianDensity
+
+
+pytest.importorskip("scipy")
+
+
+@pytest.mark.parametrize("density", [GaussianDensity(atomic_gaussian_width=1.2)])
+def test_derivative(density):
+ ppositions = np.linspace(0, 5, num=1000)
+ dens = density.compute(ppositions)
+ grad_ana = density.compute_derivative(ppositions)
+ grad_numerical = np.gradient(dens, ppositions)
+
+ assert_allclose(grad_numerical, grad_ana, atol=1e-3)
diff --git a/python/rascaline/tests/utils/splines.py b/python/rascaline/tests/utils/splines.py
index 7cb8c5de4..4915e3cfd 100644
--- a/python/rascaline/tests/utils/splines.py
+++ b/python/rascaline/tests/utils/splines.py
@@ -326,3 +326,48 @@ def test_center_contribution_gto_gaussian():
center_contr_analytical *= normalization * 2 * np.pi / np.sqrt(4 * np.pi)
assert_allclose(spliner.center_contribution, center_contr_analytical, rtol=1e-14)
+
+
+def test_custom_radial_integral():
+ cutoff = 2.0
+ max_radial = 3
+ max_angular = 2
+ atomic_gaussian_width = 1.2
+ n_spline_points = 20
+
+ spliner_args = {
+ "max_radial": max_radial,
+ "max_angular": max_angular,
+ "cutoff": cutoff,
+ "basis": GtoBasis(cutoff=cutoff, max_radial=max_radial),
+ }
+
+ spliner_analytical = SoapSpliner(
+ density=GaussianDensity(atomic_gaussian_width),
+ **spliner_args,
+ )
+
+ # Create a custom density that has not type "GaussianDensity" to trigger full
+ # numerical evaluation of the radial integral.
+ class mydensity(GaussianDensity):
+ pass
+
+ spliner_numerical = SoapSpliner(
+ density=mydensity(atomic_gaussian_width),
+ **spliner_args,
+ )
+
+ splines_analytic = spliner_analytical.compute(n_spline_points)
+ splines_numerical = spliner_numerical.compute(n_spline_points)
+
+ n_points = len(splines_analytic["TabulatedRadialIntegral"]["points"])
+
+ # TODO fix 0 component!!
+ for point in range(1, n_points):
+ point_analytical = splines_analytic["TabulatedRadialIntegral"]["points"][point]
+ point_numerical = splines_numerical["TabulatedRadialIntegral"]["points"][point]
+
+ assert point_numerical["position"] == point_analytical["position"]
+
+ for d in ["values", "derivatives"]:
+ assert_allclose(point_analytical[d]["data"], point_numerical[d]["data"])