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
Labels(
o3_lambda o3_sigma first_atom_type second_atom_type
0 1 1 1
1 1 1 1
2 1 1 1
0 1 1 6
1 1 1 6
2 1 1 6
0 1 1 7
1 1 1 7
2 1 1 7
0 1 1 8
1 1 1 8
2 1 1 8
1 -1 1 1
2 -1 1 1
3 1 1 1
1 -1 1 6
2 -1 1 6
3 1 1 6
1 -1 1 7
2 -1 1 7
3 1 1 7
1 -1 1 8
2 -1 1 8
3 1 1 8
3 -1 1 1
4 1 1 1
3 -1 1 6
4 1 1 6
3 -1 1 7
4 1 1 7
3 -1 1 8
4 1 1 8
0 1 6 1
1 1 6 1
2 1 6 1
0 1 6 6
1 1 6 6
2 1 6 6
0 1 6 7
1 1 6 7
2 1 6 7
0 1 6 8
1 1 6 8
2 1 6 8
1 -1 6 1
2 -1 6 1
3 1 6 1
1 -1 6 6
2 -1 6 6
3 1 6 6
1 -1 6 7
2 -1 6 7
3 1 6 7
1 -1 6 8
2 -1 6 8
3 1 6 8
3 -1 6 1
4 1 6 1
3 -1 6 6
4 1 6 6
3 -1 6 7
4 1 6 7
3 -1 6 8
4 1 6 8
0 1 7 1
1 1 7 1
2 1 7 1
0 1 7 6
1 1 7 6
2 1 7 6
0 1 7 7
1 1 7 7
2 1 7 7
0 1 7 8
1 1 7 8
2 1 7 8
1 -1 7 1
2 -1 7 1
3 1 7 1
1 -1 7 6
2 -1 7 6
3 1 7 6
1 -1 7 7
2 -1 7 7
3 1 7 7
1 -1 7 8
2 -1 7 8
3 1 7 8
3 -1 7 1
4 1 7 1
3 -1 7 6
4 1 7 6
3 -1 7 7
4 1 7 7
3 -1 7 8
4 1 7 8
0 1 8 1
1 1 8 1
2 1 8 1
0 1 8 6
1 1 8 6
2 1 8 6
0 1 8 7
1 1 8 7
2 1 8 7
0 1 8 8
1 1 8 8
2 1 8 8
1 -1 8 1
2 -1 8 1
3 1 8 1
1 -1 8 6
2 -1 8 6
3 1 8 6
1 -1 8 7
2 -1 8 7
3 1 8 7
1 -1 8 8
2 -1 8 8
3 1 8 8
3 -1 8 1
4 1 8 1
3 -1 8 6
4 1 8 6
3 -1 8 7
4 1 8 7
3 -1 8 8
4 1 8 8
)
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
Labels(
o3_lambda o3_sigma first_atom_type second_atom_type
0 1 1 1
1 1 1 1
2 1 1 1
0 1 1 6
1 1 1 6
2 1 1 6
0 1 1 7
1 1 1 7
2 1 1 7
0 1 1 8
1 1 1 8
2 1 1 8
1 -1 1 1
2 -1 1 1
1 -1 1 6
2 -1 1 6
1 -1 1 7
2 -1 1 7
1 -1 1 8
2 -1 1 8
0 1 6 1
1 1 6 1
2 1 6 1
0 1 6 6
1 1 6 6
2 1 6 6
0 1 6 7
1 1 6 7
2 1 6 7
0 1 6 8
1 1 6 8
2 1 6 8
1 -1 6 1
2 -1 6 1
1 -1 6 6
2 -1 6 6
1 -1 6 7
2 -1 6 7
1 -1 6 8
2 -1 6 8
0 1 7 1
1 1 7 1
2 1 7 1
0 1 7 6
1 1 7 6
2 1 7 6
0 1 7 7
1 1 7 7
2 1 7 7
0 1 7 8
1 1 7 8
2 1 7 8
1 -1 7 1
2 -1 7 1
1 -1 7 6
2 -1 7 6
1 -1 7 7
2 -1 7 7
1 -1 7 8
2 -1 7 8
0 1 8 1
1 1 8 1
2 1 8 1
0 1 8 6
1 1 8 6
2 1 8 6
0 1 8 7
1 1 8 7
2 1 8 7
0 1 8 8
1 1 8 8
2 1 8 8
1 -1 8 1
2 -1 8 1
1 -1 8 6
2 -1 8 6
1 -1 8 7
2 -1 8 7
1 -1 8 8
2 -1 8 8
)