Mitigating catastrophic cancellation error in density-based gradient calculation

[TheophaneBernhard]

Decreasing the threshold value for sin_theta for the gradient calculation of spherical harmonics seem to increase the accuracy of the computed atom-density spherical expansion coefficient for atoms close in XY (small theta angle). This suggests that the original formula is numerically more accurate than the one used when sin_theta<threshold. We think that the latter suffers from "catastrophic cancellation" of digits when subtracting pairs of associated Legendre polynomials, hitting numerical precision.

I report here the error between the analytical gradient provided by old and new featomic implementations and the gradient calculated using a central finite difference scheme performed with decreasing atomic displacements. The plots show the absolute error between representative spherical harmonics components, taking the mean over radial channels, for l=3,6 angular orders, at different m values, and for the specific case where the displaced atom is close in XY with respect to the central atom.
One can observe that the convergence is vastly improved when decreasing the threshold value.

The python script to perform these tests is attached.

analysis.py

[Luthaf]

Thanks for the PR @TheophaneBernhard! I'll run CI to make sure all tests are passing and then I can merge this.

[Luthaf]

Ok, I know what's causing the CI failure, and it is unrelated to your code. I'll fix it ASAP

[TheophaneBernhard]

Ah ok, thanks!

[Luthaf]

OK, so the CI error is relevant here, the failing test (you can reproduce locally with cargo test --test soap-power-spectrum) is a regression test for the gradients of the SOAP power spectrum.

I'll need to find some time to dig deeper in the actual error to see if it is safe to ignore or not