Assertion error in `stevens_basis_for_symmetry_allowed_anisotropies` triggered by certain crystals

[ddahlbom]

Sometimes the analysis of allowed onsite anisotropies fails, for example when calling print_symmetry_table. Specifically, one may get AssertionError: norm(imag(B)) < 1.0e-12, which is thrown here.

A minimal example which produces this error is given by an FCC lattice with rhombohedral distortion:

a = 6.22
distortion = 0.15
latvecs = lattice_vectors(a, a, a, 90+distortion, 90+distortion, 90+distortion)
positions = Sunny.fcc_crystal().positions
crystal = Crystal(latvecs, positions; types = ["A", "B", "B", "B"])
print_symmetry_table(crystal, 7.0)

The values of distortion for which this happens are hard to predict. For example, distortion=5.0 works fine, as does distortion=1. distortion=1.1 fails.

@kbarros checked that the rotation matrix derived from the lattice vectors is consistently orthogonal to high precision, as is the unitary constructed from that rotation matrix.

[kbarros]

Given the crystal above, the error can be directly triggered with this line:

Sunny.stevens_basis_for_symmetry_allowed_anisotropies(crystal, 1; k=2, R=Sunny.Mat3(I))

It looks like all the formal math in Sunny is correct, but that 12 digits of floating point accuracy cannot be maintained in the current implementation. Analysis of floating point error is difficult.

Maybe a quick workaround is to reduce the threshold to something like:

@assert norm(imag(B)) < 1e-8

and to warn the user whenever the imaginary part exceeds 1e-10?

[ddahlbom]

That seems like a reasonable solution for the moment.

I must admit that I don't have a feeling for how "off" things can get and at what point the symbolic result reported in print_symmetry_table would start to appear different due to rounding issues. It would take some experimentation.