Fix the definition of real2complex

We where using the complex conjugate of the wikipedia convention

python/featomic/featomic/clebsch_gordan/_coefficients.py

@@ -353,9 +353,9 @@ def _real2complex(o3_lambda: int, like: Array) -> Array:
     for m in range(-o3_lambda, o3_lambda + 1):
         if m < 0:
             # Positive part
-            result[o3_lambda + m, o3_lambda + m] = i_sqrt_2
+            result[o3_lambda + m, o3_lambda + m] = -i_sqrt_2
             # Negative part
-            result[o3_lambda + m, o3_lambda - m] = -i_sqrt_2 * ((-1) ** m)
+            result[o3_lambda + m, o3_lambda - m] = i_sqrt_2 * ((-1) ** m)
 
         if m == 0:
             result[o3_lambda, o3_lambda] = 1.0

python/featomic/tests/clebsch_gordan/coefficients.py

@@ -0,0 +1,50 @@
+import numpy as np
+import pytest
+
+from featomic.clebsch_gordan._coefficients import _complex2real
+
+
+scipy = pytest.importorskip("scipy")
+
+
+def complex_to_real_manual(sph):
+    # following https://en.wikipedia.org/wiki/Spherical_harmonics#Real_form
+    ell = (sph.shape[1] - 1) // 2
+
+    real = np.zeros(sph.shape)
+    for m in range(-ell, ell + 1):
+        if m < 0:
+            real[:, ell + m] = np.sqrt(2) * (-1) ** m * np.imag(sph[:, ell + abs(m)])
+        elif m == 0:
+            assert np.all(np.imag(sph[:, ell + m]) == 0)
+            real[:, ell + m] = np.real(sph[:, ell + m])
+        else:
+            real[:, ell + m] = np.sqrt(2) * (-1) ** m * np.real(sph[:, ell + m])
+
+    return real
+
+
+def complex_to_real_matrix(sph):
+    ell = (sph.shape[1] - 1) // 2
+
+    matrix = _complex2real(ell, sph)
+
+    real = sph @ matrix
+
+    assert np.linalg.norm(np.imag(real)) < 1e-15
+    return np.real(real)
+
+
+def test_complex_to_real():
+    theta = 2 * np.pi * np.random.rand(10)
+    phi = np.pi * np.random.rand(10)
+
+    for ell in range(4):
+        values = np.zeros((10, 2 * ell + 1), dtype=np.complex128)
+        for m in range(-ell, ell + 1):
+            values[:, ell + m] = scipy.special.sph_harm(m, ell, theta, phi)
+
+        real_manual = complex_to_real_manual(values)
+        real_matrix = complex_to_real_matrix(values)
+
+        assert np.allclose(real_manual, real_matrix)