Computing per-pair equivariant features¶
In a previous example, we computed the \(\lambda\)-SOAP equivariant descriptor. This tutorial focuses on the generalization of \(\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.
import chemfiles
from metatensor import Labels
from featomic import SphericalExpansion, SphericalExpansionByPair
from featomic.clebsch_gordan import EquivariantPowerSpectrumByPair
Let’s see how to compute the per-pair \(\lambda\)-SOAP descriptor using Featomic.
Read systems using Chemfiles. You can download the dataset for this example from our
website.
with chemfiles.Trajectory("dataset.xyz") as trajectory:
systems = [s for s in trajectory]
Featomic also handles systems read by 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 a spherical expansion and a spherical expansion by pair calculator.
The SphericalExpansion and
SphericalExpansionByPair classes use the hyperparameters above.
Then, wrap them with EquivariantPowerSpectrumByPair
to compute the Clebsch-Gordan contraction for the per-pair \(\lambda\)-SOAP.
spex_calculator = SphericalExpansion(**HYPERPARAMETERS)
per_pair_spex_calculator = SphericalExpansionByPair(**HYPERPARAMETERS)
calculator = EquivariantPowerSpectrumByPair(spex_calculator, per_pair_spex_calculator)
Run the actual calculation
per_pair_power_spectrum = calculator.compute(systems, neighbors_to_properties=True)
The result is a TensorMap whose keys encode symmetry and the
species of the atoms involved:
per_pair_power_spectrum.keys
Often, you only need specific \(\lambda\) values. For example, if the
target property is the Hamiltonian matrix on a minimal basis, you can restrict the output to \(\lambda\) values
up to \(\lambda=2\) using the selected_keys parameter:
per_pair_power_spectrum_minimal_basis = calculator.compute(
systems,
neighbors_to_properties=True,
selected_keys=Labels.range("o3_lambda", 3),
)
per_pair_power_spectrum_minimal_basis.keys