## Use sparse force constants and exploit the label symmetries.

Right now, calculating the contribution from a cluster to the total energy (or forces) takes roughly `\mathcal{O}(n3^n)`

flops where `n`

is the order. However, consider the cluster `(iiii)`

. This cluster has only 15 independent components compared to 81 for a four body 4th order cluster `ijkl`

. These internal symmetries are not exploited in the code at the moment (we do however use the fact that e.g. cluster `(ijk)`

is equivalent to `(ikj)`

).

## Example

Since we are often interested in high order - low n-body interactions this would benefit both the MD and fit-matrix calculation as well as save a lot of memory. The difference is of course even larger for sixth order etc.

Consider the cluster `iiiijjjj`

. It is a two body, 8th order interaction. The raw format would take around 50kB while the actual number of independent components is only 2kB *before* any symmetry constraint. And as the computational work is roughly equivalent to the memory footprint I think this could increase the performance with about one to two orders of magnitude.

Note that the above is true for a P1 symmetry system and the benefit could be even greater for e.g. Titanium.

## Implementation

Depending on the level of ambition some tricks can be used to speed up things. To me, they are all not trivial though. With the sparse format in place, changing the actual contraction function should be easy. For the simple `(iii)`

cluster the data would look like:

```
cluster = (i, i, i)
indices = [
[x, x, x],
[x, x, y],
[x, x, z],
[x, y, y],
[x, y, z],
[x, z, z],
[y, y, y],
[y, y, z],
[y, z, z],
[z, z, z],
]
multiplicities = [1, 3, 3, 3, 6, 3, 1, 3, 3, 1]
```

and the corresponding list of fcis

` \text{fcis = [} \Phi^{xxx}, \Phi^{xxy}, ..., \Phi^{yzz}, \Phi^{zzz} \text{]}`

and the contraction would be

```
E = 0
for multi_index, multiplicity, fci in zip(indices, multiplicities, fcis):
u = 1.0
for index, atom in zip(multi_index, cluster):
u *= displacements[atom, index]
E += fci * multiplicity * u
```

The force calculation is a bit more complicated depending on level of ambition but goes along the same idea.