Tutorials for high-level calculators from the clebsch_gordan module

docs/src/how-to/computing-lambda-soap.rst

@@ -0,0 +1,22 @@
+.. _userdoc-how-to-computing-lambda-soap:
+
+Computing :math:`\lambda`-SOAP features
+=======================================
+
+In a `previous example
+<https://metatensor.github.io/featomic/latest/how-to/computing-soap.html>`_, we
+computed the SOAP **scalar** descriptor. This tutorial focuses on the
+equivariant generalization of SOAP to **tensorial** objects, commonly known as
+:math:`\lambda`-SOAP. The key difference is that the SOAP scalar is the inner product
+of two spherical expansions, whereas :math:`\lambda`-SOAP contracts two spherical
+expansions using a Clebsch-Gordan coefficient:
+
+.. math::
+
+  \rho^{\otimes 2}_{z_1 z_2, n_1 n_2, l_1, l_2, \lambda \mu} (\{\mathbf{r}_i\}) =
+  \sum_{m_1, m_2} C_{l_1 l_2, m_1 m_2}^{\lambda\mu}
+  \rho_{z_1,n_1 l_1 m_1}(\{\mathbf{r}_i\}) \rho_{z_2,n_2 l_2 m_2}(\{\mathbf{r}_i\})
+
+.. include:: ../examples/compute-lambda-soap.rst
+    :start-after: start-body
+    :end-before: end-body

docs/src/how-to/density-correlations.rst

@@ -0,0 +1,18 @@
+.. _userdoc-how-to-density-correlations:
+
+Computing density correlations
+==============================
+
+In a `previous example
+<https://metatensor.github.io/featomic/latest/how-to/computing-lambda-soap.html>`_,
+we computed the :math:`\lambda`-SOAP **equivariant** descriptor. This tutorial
+focuses on the computation of arbitrary body-order descriptors by means of an
+auto-correlations of spherical expansion (or 'density').
+
+We will explore both :math:`\lambda`-SOAP --- a version of auto-correlations
+with body-order 3 (taking into account one central atom and two neighbors) ---
+an the bispectrum which is a representation with body-order 4.
+
+.. include:: ../examples/density-correlations.rst
+    :start-after: start-body
+    :end-before: end-body

docs/src/how-to/index.rst

@@ -11,8 +11,11 @@ are a total beginner, you can go to :ref:`userdoc-get-started` section.
     :maxdepth: 1
 
     computing-soap
+    computing-lambda-soap
+    density-correlations
     sample-selection
     property-selection
     keys-selection
     splined-radial-integral
     long-range
+    per-pair-equivariant-features

docs/src/how-to/per-pair-equivariant-features.rst

@@ -0,0 +1,15 @@
+.. _userdoc-how-to-per-pair-equivariant-features:
+
+Computing per-pair equivariant features
+=======================================
+
+In a `previous example
+<https://metatensor.github.io/featomic/latest/how-to/computing-lambda-soap.html>`_, we
+computed the :math:`\lambda`-SOAP equivariant descriptor. This tutorial focuses on the
+generalization of :math:`\lambda`-SOAP to pairs of atoms rather than single centers. 
+Per-pair equivariant features are useful as descriptors for `two-center quantities such 
+as the Hamiltonian matrix <https://doi.org/10.1063/5.0072784>`_.
+
+.. include:: ../examples/per-pair-equivariant-features.rst
+    :start-after: start-body
+    :end-before: end-body

python/featomic/examples/compute-lambda-soap.py

@@ -0,0 +1,111 @@
+"""
+.. _compute-lambda-soap:
+
+Computing  λ-SOAP features
+==========================
+
+.. start-body
+"""
+
+import chemfiles
+import numpy as np
+from metatensor import Labels
+
+from featomic import LodeSphericalExpansion, SphericalExpansion
+from featomic.clebsch_gordan import EquivariantPowerSpectrum
+
+
+# %%
+#
+# Let's see how to compute the :math:`\lambda`-SOAP descriptor using featomic.
+#
+# First we can read the input systems using chemfiles. You can download the dataset for
+# this example from our :download:`website <../../static/dataset.xyz>`.
+
+with chemfiles.Trajectory("dataset.xyz") as trajectory:
+    systems = [s for s in trajectory]
+
+# %%
+#
+# Featomic also handles systems read by `ASE <https://wiki.fysik.dtu.dk/ase/>`_:
+#
+# ``systems = ase.io.read("dataset.xyz", ":")``.
+#
+# Next, define the hyperparameters for the spherical expansion:
+
+HYPERPARAMETERS = {
+    "cutoff": {
+        "radius": 5.0,
+        "smoothing": {"type": "ShiftedCosine", "width": 0.5},
+    },
+    "density": {
+        "type": "Gaussian",
+        "width": 0.3,
+    },
+    "basis": {
+        "type": "TensorProduct",
+        "max_angular": 2,
+        "radial": {"type": "Gto", "max_radial": 2},
+    },
+}
+
+# %%
+# Create the spherical expansion calculator. The :class:`~featomic.SphericalExpansion`
+# class uses the hyperparameters above. Then, wrap it with
+# :class:`~featomic.clebsch_gordan.EquivariantPowerSpectrum` to compute the
+# Clebsch-Gordan contraction for :math:`\lambda`-SOAP.
+
+spex_calculator = SphericalExpansion(**HYPERPARAMETERS)
+calculator = EquivariantPowerSpectrum(spex_calculator)
+
+# %%
+# Run the actual calculation
+
+power_spectrum = calculator.compute(systems, neighbors_to_properties=True)
+
+# %%
+# The result is a :class:`~metatensor.TensorMap` whose keys encode symmetry:
+
+power_spectrum.keys
+
+# %%
+# Often, you only need specific :math:`\lambda` values. For example, if the
+# `target property is the polarizability tensor
+# <https://atomistic-cookbook.org/examples/polarizability/polarizability.html>`_,
+# (a rank-2 symmetric Cartesian tensor), you can restrict the output to
+# :math:`\lambda=0` and :math:`\lambda=2` (with :math:`\sigma=1` to discard
+# `pseudotensors <https://en.wikipedia.org/wiki/Pseudotensor>`_) using the
+# ``selected_keys`` parameter:
+
+power_spectrum_0_2 = calculator.compute(
+    systems,
+    neighbors_to_properties=True,
+    selected_keys=Labels(["o3_lambda", "o3_sigma"], np.array([[0, 1], [2, 1]])),
+)
+power_spectrum_0_2.keys
+
+# %%
+# You can also compute a :math:`\lambda`-SOAP-like descriptor using two different
+# expansions. For instance, combine a standarad spherical expansion with a long-range
+# :class:`~featomic.LodeSphericalExpansion`:
+
+LODE_HYPERPARAMETERS = {
+    "density": {
+        "type": "SmearedPowerLaw",
+        "smearing": 0.3,
+        "exponent": 1,
+    },
+    "basis": {
+        "type": "TensorProduct",
+        "max_angular": 2,
+        "radial": {"type": "Gto", "max_radial": 3, "radius": 1.0},
+    },
+}
+lode_calculator = LodeSphericalExpansion(**LODE_HYPERPARAMETERS)
+calculator = EquivariantPowerSpectrum(spex_calculator, lode_calculator)
+power_spectrum = calculator.compute(systems, neighbors_to_properties=True)
+power_spectrum.keys
+
+# %%
+#
+# .. end-body