Add how-to to construct pure long range descriptor

docs/src/how-to/index.rst

@@ -16,3 +16,4 @@ are a total beginner, you can go to :ref:`userdoc-get-started` section.
     keys-selection
     le-basis
     splined-radial-integral
+    long-range

docs/src/how-to/long-range.rst

@@ -0,0 +1,45 @@
+.. _userdoc-how-to-long-range-descriptor:
+
+Pure LODE long-range descriptor
+===============================
+
+The :py:class:`rascaline.LodeSphericalExpansion` allows the calculation of a descriptor that
+includes all atoms within the system and projects them onto a spherical expansion/
+fingerprint within a given ``cutoff``. This is very useful if long-range interactions
+between atoms are important to describe the physics and chemistry of a collection of
+atoms. However, as stated the descriptor contains **ALL** atoms of the system and
+sometimes it can be desired to only have a pure long-range/exterior descriptor that only
+includes the atoms outside a given cutoff.
+
+.. figure:: ../../static/images/long-range-descriptor.*
+    :align: center
+    :class: only-light
+
+.. figure:: ../../static/images/long-range-descriptor-dark.*
+    :align: center
+    :class: only-dark
+
+This descriptor can be particular useful when already having a good descriptor for the
+short-range density like (SOAP) and the long-range descriptor (far field) should not
+have the information the short-range descriptor already contains. 
+
+In this example will construct such a descriptor using the :ref:`radial integral
+splining <python-utils-splines>` tools of ``rascaline``.
+
+.. tabs::
+
+    .. group-tab:: Python
+
+        .. container:: sphx-glr-footer sphx-glr-footer-example
+
+            .. container:: sphx-glr-download sphx-glr-download-python
+
+                :download:`Download Python source code for this example: long-range-descriptor.py <../examples/long-range-descriptor.py>`
+
+            .. container:: sphx-glr-download sphx-glr-download-jupyter
+
+                :download:`Download Jupyter notebook for this example: long-range-descriptor.ipynb <../examples/long-range-descriptor.ipynb>`
+
+        .. include:: ../examples/long-range-descriptor.rst
+            :start-after: start-body
+            :end-before: end-body

docs/src/references/api/python/utils/splines.rst

@@ -1 +1,3 @@
+.. _python-utils-splines:
+
 .. automodule:: rascaline.utils.splines.splines

python/rascaline/examples/long-range-descriptor.py

@@ -0,0 +1,211 @@
+"""
+.. _userdoc-tutorials-long-range-descriptor:
+
+Pure LODE long-range descriptor
+===============================
+
+.. start-body
+
+We start the example by loading the required packages
+"""
+
+import ase
+import ase.visualize.plot
+import matplotlib.pyplot as plt
+import numpy as np
+from ase.build import molecule
+from numpy.testing import assert_allclose
+
+from rascaline import LodeSphericalExpansion, SphericalExpansion
+from rascaline.utils import GtoBasis, LodeDensity, LodeSpliner, SoapSpliner
+
+
+# %%
+# Our test system is a single water molecule with a large vacuum layer around it.
+
+atoms = molecule("H2O", vacuum=10, pbc=True)
+
+ase.visualize.plot.plot_atoms(atoms)
+
+plt.xlabel("Å")
+plt.ylabel("Å")
+
+plt.show()
+
+
+# %%
+# We choose a ``cutoff`` that includes all atoms of the system for the projection of the
+# spherical expansion and the neighbor search of the real space spherical expansion.
+#
+
+cutoff = 5
+
+
+# %%
+# The idea of the combination of the test system and the ``cutoff`` is to in the end
+# be able to show that the full atomic fingerprint is contained within the ``cutoff``.
+# By later subtracting the LODE density and the short-range density we will see that the
+# difference between them is zero - meaning that this is the short-range system.
+#
+# On top of that we choose a large potential exponent to have a quickly decaying LODE
+# density.
+
+potential_exponent = 6
+
+
+# %%
+# We also define some usual hyperparameters to compute the spherical expansions.
+#
+
+max_radial = 5
+max_angular = 1
+atomic_gaussian_width = 1.2
+recalc_splines = True
+center_atom_weight = 1.0
+
+
+# %%
+# We choose a relatively low spline accuracy (default is ``1e-8``) to achieve a quick
+# compuyattion of the spline points. You can increase spline accuracy if required,
+# computation time will increase!
+
+spline_accuracy = 1e-2
+
+
+# %%
+# As a projection basis we choose the usual ``GTO`` Basis
+
+basis = GtoBasis(cutoff=cutoff, max_radial=max_radial)
+
+
+# %%
+# And as the basis we choose a smeared power law, as used in LODE. Which does not decay
+# exponentially as for example a :py:class:`Gaussian density
+# <rascaline.utils.GaussianDensity>` and therefore is suited to describe long range
+# interactions between atoms.
+
+density = LodeDensity(
+    atomic_gaussian_width=atomic_gaussian_width,
+    potential_exponent=potential_exponent,
+)
+
+
+# %%
+# Now we have all the building blocks to construct the spline points for the real and
+# Fourier space Spherical expansions.
+
+rs_splines = SoapSpliner(
+    cutoff=cutoff,
+    max_radial=max_radial,
+    max_angular=max_angular,
+    basis=basis,
+    density=density,
+    accuracy=spline_accuracy,
+).compute()
+
+
+# Usually this value for `k_cutoff` gives good convergences for the k-space version
+k_cutoff = 1.2 * np.pi / atomic_gaussian_width
+
+# Fourier space splines
+fs_splines = LodeSpliner(
+    k_cutoff=k_cutoff,
+    max_radial=max_radial,
+    max_angular=max_angular,
+    basis=basis,
+    density=density,
+    accuracy=spline_accuracy,
+).compute()
+
+
+# %%
+# You might want to save the spline points using :py:func:`json.dump` to a file to reuse
+# them later without recalculating them again. Saving them is especially useful if the
+# spline calculations are expensive i.e. if you increase the ``spline_accuracy``.
+#
+# With the spline points on hand we now compute the real space spherical expansion
+
+rs_se = SphericalExpansion(
+    cutoff=cutoff,
+    max_radial=max_radial,
+    max_angular=max_angular,
+    atomic_gaussian_width=atomic_gaussian_width,
+    radial_basis=rs_splines,
+    center_atom_weight=center_atom_weight,
+    cutoff_function={"Step": {}},
+).compute(atoms)
+
+
+# %%
+# and the Fourier Space / Lode Spherical expansion.
+
+fs_se = LodeSphericalExpansion(
+    cutoff=cutoff,
+    max_radial=max_radial,
+    max_angular=max_angular,
+    atomic_gaussian_width=atomic_gaussian_width,
+    center_atom_weight=center_atom_weight,
+    potential_exponent=potential_exponent,
+    radial_basis=fs_splines,
+    k_cutoff=k_cutoff,
+).compute(atoms)
+
+
+# %%
+# As described in the beginning we now subtract real space LODE contributions from
+# fourier space.
+
+
+subtracted_se = fs_se - rs_se
+
+
+# %%
+# We densify the TensorMap to have only one block per ``"center_type"`` to visualize our
+# result.
+
+fs_se = fs_se.components_to_properties("o3_mu")
+fs_se = fs_se.keys_to_samples("neighbor_type")
+fs_se = fs_se.keys_to_properties(["o3_lambda", "o3_sigma"])
+
+subtracted_se = subtracted_se.components_to_properties("o3_mu")
+subtracted_se = subtracted_se.keys_to_samples("neighbor_type")
+subtracted_se = subtracted_se.keys_to_properties(["o3_lambda", "o3_sigma"])
+
+
+# %%
+# Finally, we plot the values of each block for the Fourier Space spherical expansion
+# the upper panel and the difference between the Fourier Space and the real space in
+# the lower panel.
+
+for key in fs_se.keys:
+    fig, ax = plt.subplots(2, layout="constrained")
+
+    values_subtracted = subtracted_se[key].values
+    values_fs = fs_se[key].values
+
+    # maximal and minimal values for same scaling of axes
+    vmin = min(np.min(values_subtracted), np.min(values_fs))
+    vmax = max(np.max(values_subtracted), np.max(values_fs))
+
+    ax[0].set_title(f"center_type={key.values[0]}\n Fourier space sph. expansion")
+    im = ax[0].matshow(values_fs, vmin=vmin, vmax=vmax)
+    ax[0].set_ylabel("sample index")
+
+    ax[1].set_title("Difference between Fourier and real space sph. expansion")
+    ax[1].matshow(values_subtracted, vmin=vmin, vmax=vmax)
+    ax[1].set_ylabel("sample index")
+    ax[1].set_xlabel("property index")
+
+    fig.colorbar(im, ax=ax[0], orientation="horizontal", fraction=0.1, label="values")
+
+    assert_allclose(values_subtracted, 0, atol=spline_accuracy)
+
+
+# %%
+# The plot shows that the spherical expansion for the Fourier Space is non-zero while
+# the difference between the two spherical expansions is zero within the given
+# tolerance. In addition to the visual inspection, we also added an assertion test using
+# :py:func:`assert_allclose <numpy.testing.assert_allclose>` to verify that the values
+# are close actually to zero.
+#
+# .. end-body