diff --git a/featomic/CHANGELOG.md b/featomic/CHANGELOG.md index 769b6ffa6..1d5e968ab 100644 --- a/featomic/CHANGELOG.md +++ b/featomic/CHANGELOG.md @@ -17,6 +17,10 @@ a changelog](https://keepachangelog.com/en/1.1.0/) format. This project follows ### Removed --> +### Fixed +- Fixed `featomic.clebsch_gordan.cartesian_to_spherical` transformation for tensors of rank greater than 2 (#371) +- Fixed the calculation of Clebsch-Gordan coefficients, used all across the `featomic.clebsch_gordan` module (#371) + ## [Version 0.6.0](https://github.com/metatensor/featomic/releases/tag/featomic-v0.6.0) - 2024-12-20 ### Added diff --git a/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py b/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py index 32fc7eba2..2c9b678d0 100644 --- a/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py +++ b/python/featomic/featomic/clebsch_gordan/_cartesian_spherical.py @@ -245,6 +245,10 @@ def cartesian_to_spherical( # [..., 3, 5, ...] => [..., 3, ...] (o3_lambda=1, o3_sigma=+1) # => [..., 5, ...] (o3_lambda=2, o3_sigma=-1) # => [..., 7, ...] (o3_lambda=3, o3_sigma=+1) + + # Use the rolled block values as the starting point for higher-rank tensors + tensor = TensorMap(tensor.keys, new_blocks) + if cg_coefficients is None: if torch_jit_is_scripting(): raise ValueError( diff --git a/python/featomic/featomic/clebsch_gordan/_coefficients.py b/python/featomic/featomic/clebsch_gordan/_coefficients.py index e7d94f336..bc134d63f 100644 --- a/python/featomic/featomic/clebsch_gordan/_coefficients.py +++ b/python/featomic/featomic/clebsch_gordan/_coefficients.py @@ -197,23 +197,53 @@ def _build_dense_cg_coeff_dict( device=device, dtype=complex_like.dtype, ) - - real_cg = (r2c[l1] @ complex_cg.reshape(2 * l1 + 1, -1)).reshape( - complex_cg.shape - ) - - real_cg = real_cg.swapaxes(0, 1) - real_cg = (r2c[l2] @ real_cg.reshape(2 * l2 + 1, -1)).reshape( - real_cg.shape - ) - real_cg = real_cg.swapaxes(0, 1) - - real_cg = real_cg @ c2r[o3_lambda] + # We want to create a CG coefficient for a real representation, using + # the expression: + # + # C_{l1m1l2m2}^{lm} + # = \sum_{m'm1'm2'} ... + # ... U_{mm'}^{l} C_{l1m1'm2m2'}^{lm'} ... + # ... U^{\dagger}_{m1'm1}^{l1} U^{\dagger}_{m2'm2}^{l2} + # + # where: + # U is the "c2r" transformation matrix below + # U^{\dagger} is the "r2c" transformation matrix below + # + # 0) We have a CG coefficient of shape: + # complex_cg.shape = (2*l1+1, 2*l2+1, 2*L+1) + # and transormation matrices of shape: + # c2r[L].shape = (2*L+1, 2*L+1) + # r2c[L].shape = (2*L+1, 2*L+1) + # + # 1) take the matrix product: + # \sum_{m1'} C_{l1m1'm2m2'}^{lm'} U^{\dagger}_{m1'm1}^{l1} + # this requires some permuting of axes of the objects involved to + # ensure proper matching for matmul. i.e. permute axes of the complex + # CG coefficient from: (m1, m2, m) -> (m, m2, m1) + first_step = _dispatch.permute(complex_cg, [2, 1, 0]) @ r2c[l1] + # The result `first_step` has shape (2*L+1, 2*l2+1, 2*l1+1) + # + # 2) take the matrix product: + # \sum_{m2'} ... + # ... (first_step)_{m',m1,m2'}^{l,l1} U^{\dagger}_{m2'm2}^{l2} + first_step_swap = first_step.swapaxes(1, 2) + # The result `first_step_swap` has shape (2*L+1, 2*l1+1, 2*l2+1) + second_step = first_step_swap @ r2c[l2] + # The result `second_step` has shape (2*L+1, 2*l1+1, 2*l2+1) + # + # 3) take the matrix product: + # \sum_{m'} U_{mm'}^{l} (second_step)_{m',m1,m2}^{l,l1,l2} + second_step_swap = _dispatch.permute(second_step, [1, 0, 2]) + # The result `second_step_swap` has shape (2*l1+1, 2*L+1, 2*l2+1) + third_step = c2r[o3_lambda] @ second_step_swap + # The result `third_step` has shape (2*l1+1, 2*L+1, 2*l2+1) + third_step_swap = _dispatch.permute(third_step, [0, 2, 1]) + # The result `third_step_swap` has shape (2*l1+1, 2*l2+1, 2*L+1) if (l1 + l2 + o3_lambda) % 2 == 0: - cg_l1l2lam_dense = _dispatch.real(real_cg) + cg_l1l2lam_dense = _dispatch.real(third_step_swap) else: - cg_l1l2lam_dense = _dispatch.imag(real_cg) + cg_l1l2lam_dense = _dispatch.imag(third_step_swap) coeff_dict[(l1, l2, o3_lambda)] = _dispatch.to( cg_l1l2lam_dense, @@ -366,7 +396,7 @@ def _real2complex(o3_lambda: int, like: Array) -> Array: # Positive part result[o3_lambda + m, o3_lambda + m] = inv_sqrt_2 * ((-1) ** m) - return result + return result.T def _complex2real(o3_lambda: int, like) -> Array: @@ -394,10 +424,10 @@ def cg_couple( Go from an uncoupled product basis that behave like a product of spherical harmonics to a coupled basis that behaves like a single spherical harmonic. - The ``array`` shape should be ``(n_samples, 2 * l1 + 1, 2 * l2 + 1, n_q)``. - ``n_samples`` is the number of samples, and ``n_q`` is the number of properties. + The ``array`` shape should be ``(n_s, 2 * l1 + 1, 2 * l2 + 1, n_q)``. + ``n_s`` is the number of samples, and ``n_q`` is the number of properties. - This function will output a list of arrays, whose shape will be ``[n_samples, (2 * + This function will output a list of arrays, whose shape will be ``[n_s, (2 * o3_lambda+1), n_q]``, with the requested ``o3_lambdas``. These arrays will contain the result of projecting from a product of spherical @@ -411,7 +441,7 @@ def cg_couple( The operation is dispatched such that numpy arrays or torch tensors are automatically handled. - :param array: input array with shape ``[n_samples, 2 * l1 + 1, 2 * l2 + 1, n_q]`` + :param array: input array with shape ``[n_s, 2 * l1 + 1, 2 * l2 + 1, n_q]`` :param o3_lambdas: list of degrees of spherical harmonics to compute :param cg_coefficients: CG coefficients as returned by :py:func:`calculate_cg_coefficients` with the same ``cg_backed`` given to this @@ -438,18 +468,24 @@ def cg_couple( for o3_lambda in o3_lambdas ] elif cg_backend == "python-dense": - results = [] - n_samples = array.shape[0] n_properties = array.shape[3] - array = array.swapaxes(1, 3) - array = array.reshape(n_samples * n_properties, 2 * l2 + 1, 2 * l1 + 1) + # Move properties next to samples + array = _dispatch.permute(array, [0, 3, 1, 2]) + array = array.reshape(n_samples * n_properties, 2 * l1 + 1, 2 * l2 + 1) + + results = [] for o3_lambda in o3_lambdas: result = _cg_couple_dense(array, o3_lambda, cg_coefficients) + + # Separate samples and properties result = result.reshape(n_samples, n_properties, -1) + + # Back to TensorBlock-like axes result = result.swapaxes(1, 2) + results.append(result) return results @@ -499,7 +535,11 @@ def _cg_couple_sparse( m2 = int(m1m2mu[1]) mu = int(m1m2mu[2]) # Broadcast arrays, multiply together and with CG coeff - output[mu, :, :] += arrays[str((m1, m2))] * cg_l1l2lam.values[i, 0] + output[mu, :, :] += ( + arrays[str((m1, m2))] + * (-1) ** (l1 + l2 + o3_lambda) + * cg_l1l2lam.values[i, 0] + ) return output.swapaxes(0, 1) @@ -514,26 +554,23 @@ def _cg_couple_dense( degree ``o3_lambda``) using CG coefficients. This is a "dense" implementation, using all CG coefficients at the same time. - :param array: input array, we expect a shape of ``[samples, 2*l1 + 1, 2*l2 + 1]`` + :param array: input array, we expect a shape of ``[N, 2*l1 + 1, 2*l2 + 1]`` :param o3_lambda: value of lambda for the output spherical harmonic :param cg_coefficients: CG coefficients as returned by :py:func:`calculate_cg_coefficients` with ``cg_backed="python-dense"`` + + :return: array of shape ``[N, 2 * o3_lambda + 1]`` """ assert len(array.shape) == 3 l1 = (array.shape[1] - 1) // 2 l2 = (array.shape[2] - 1) // 2 - cg_l1l2lam = cg_coefficients.block({"l1": l1, "l2": l2, "lambda": o3_lambda}).values - - # [samples, l1, l2] => [samples, (l1 l2)] - array = array.reshape(-1, (2 * l1 + 1) * (2 * l2 + 1)) - - # [l1, l2, lambda] -> [(l1 l2), lambda] - cg_l1l2lam = cg_l1l2lam.reshape(-1, 2 * o3_lambda + 1) + cg_l1l2lam = (-1) ** (l1 + l2 + o3_lambda) * cg_coefficients.block( + {"l1": l1, "l2": l2, "lambda": o3_lambda} + ).values - # [samples, (l1 l2)] @ [(l1 l2), lambda] => [samples, lambda] - return array @ cg_l1l2lam + return _dispatch.tensordot(array, cg_l1l2lam[0, ..., 0], axes=([2, 1], [1, 0])) # ======================================================================= # diff --git a/python/featomic/featomic/clebsch_gordan/_dispatch.py b/python/featomic/featomic/clebsch_gordan/_dispatch.py index dc097df56..43d9a8d01 100644 --- a/python/featomic/featomic/clebsch_gordan/_dispatch.py +++ b/python/featomic/featomic/clebsch_gordan/_dispatch.py @@ -63,6 +63,34 @@ def concatenate(arrays: List[TorchTensor], axis: int): raise TypeError(UNKNOWN_ARRAY_TYPE) +def permute(array, axes: List[int]): + """ + Permute the axes of the given ``array``. For numpy and torch, ``transpose`` and + ``permute`` are equivalent operations. + """ + if isinstance(array, TorchTensor): + return torch.permute(array, axes) + elif isinstance(array, np.ndarray): + return np.transpose(array, axes) + else: + raise TypeError(UNKNOWN_ARRAY_TYPE) + + +def tensordot(array_1, array_2, axes: List[List[int]]): + """ + Compute tensor dot product along specified axes for arrays. + + This function has the same behavior as ``np.tensordot(array_1, array_2, axes)`` and + ``torch.tensordot(array_1, array_2, axes)``. + """ + if isinstance(array_1, TorchTensor): + return torch.tensordot(array_1, array_2, axes) + elif isinstance(array_1, np.ndarray): + return np.tensordot(array_1, array_2, axes) + else: + raise TypeError(UNKNOWN_ARRAY_TYPE) + + def empty_like(array, shape: Optional[List[int]] = None, requires_grad: bool = False): """ Create an uninitialized array, with the given ``shape``, and similar dtype, device diff --git a/python/featomic/tests/clebsch_gordan/cartesian_spherical.py b/python/featomic/tests/clebsch_gordan/cartesian_spherical.py index cffdfc685..16d805d74 100644 --- a/python/featomic/tests/clebsch_gordan/cartesian_spherical.py +++ b/python/featomic/tests/clebsch_gordan/cartesian_spherical.py @@ -1,9 +1,12 @@ +import metatensor as mts import numpy as np import pytest from metatensor import Labels, TensorBlock, TensorMap from featomic.clebsch_gordan import cartesian_to_spherical +from .rotations import WignerDReal, cartesian_rotation + @pytest.fixture def cartesian(): @@ -198,3 +201,278 @@ def test_cartesian_to_spherical_errors(cartesian): components=["xyz_1", "xyz_2", "xyz_3"], keep_l_in_keys=False, ) + + +def _l1_components_from_matrix(A): + """ + *The* reference equations for projecting a cartesian rank 2 tensor to the + irreducible spherical component with lambda = 1 and sigma = -1. + """ + A = A.reshape(3, 3) + + l1_A = np.empty((3,)) + l1_A[0] = A[2, 0] - A[0, 2] # y + l1_A[1] = A[0, 1] - A[1, 0] # z + l1_A[2] = A[1, 2] - A[2, 1] # x + + return -l1_A / np.sqrt(2) + + +def _l2_components_from_matrix(A): + """ + *The* reference equations for projecting a cartesian rank 2 tensor to the + irreducible spherical component with lambda = 2. + """ + A = A.reshape(3, 3) + + l2_A = np.empty((5,)) + l2_A[0] = (A[0, 1] + A[1, 0]) / 2.0 + l2_A[1] = (A[1, 2] + A[2, 1]) / 2.0 + l2_A[2] = (2.0 * A[2, 2] - A[0, 0] - A[1, 1]) / ((2.0) * np.sqrt(3.0)) + l2_A[3] = (A[0, 2] + A[2, 0]) / 2.0 + l2_A[4] = (A[0, 0] - A[1, 1]) / 2.0 + + return l2_A * np.sqrt(2) + + +def _l3_components_from_matrix(A): + """ + *The* reference equations for projecting a cartesian rank 3 tensor to the + irreducible spherical component with lambda = 3. + """ + A = A.reshape(3, 3, 3) + + l3_A = np.empty((7,)) + + # o3_mu = -3 + l3_A[0] = (A[0, 1, 0] + A[1, 0, 0] + A[0, 0, 1] - A[1, 1, 1]) / 2.0 + + # o3_mu = -2 + l3_A[1] = ( + A[0, 1, 2] + A[1, 0, 2] + A[1, 2, 0] + A[2, 1, 0] + A[0, 2, 1] + A[2, 0, 1] + ) / np.sqrt(6) + + # o3_mu = -1 + l3_A[2] = ( + 4.0 * A[1, 2, 2] + + 4.0 * A[2, 1, 2] + + 4.0 * A[2, 2, 1] + - 3.0 * A[1, 1, 1] + - A[0, 0, 1] + - A[0, 1, 0] + - A[1, 0, 0] + ) / np.sqrt(60) + + # o3_mu = 0 + l3_A[3] = ( + 2.0 * A[2, 2, 2] + - A[0, 2, 0] + - A[2, 0, 0] + - A[0, 0, 2] + - A[1, 2, 1] + - A[2, 1, 1] + - A[1, 1, 2] + ) / np.sqrt(10) + + # o3_mu = 1 + l3_A[4] = ( + 4.0 * A[0, 2, 2] + + 4.0 * A[2, 0, 2] + + 4.0 * A[2, 2, 0] + - 3.0 * A[0, 0, 0] + - A[1, 1, 0] + - A[0, 1, 1] + - A[1, 0, 1] + ) / np.sqrt(60) + + # o3_mu = 2 + l3_A[5] = ( + A[0, 0, 2] + A[0, 2, 0] + A[2, 0, 0] - A[1, 1, 2] - A[1, 2, 1] - A[2, 1, 1] + ) / np.sqrt(6) + + # o3_mu = 3 + l3_A[6] = (A[0, 0, 0] - A[1, 1, 0] - A[0, 1, 1] - A[1, 0, 1]) / 2.0 + + return l3_A + + +@pytest.mark.parametrize("cg_backend", ["python-dense", "python-sparse"]) +def test_cartesian_to_spherical_rank_2_by_equation(cg_backend): + """ + Tests cartesian_to_spherical for a random rank-2 tensor by comparing the result to + the result from *the* reference equation. + """ + # Build the reference (lambda, sigma) = (1, -1) and (lambda, sigma) = (2, 1) + # components + random_rank_2_arr = np.random.rand(100, 3, 3, 1) + l1_and_l2_reference = TensorMap( + keys=Labels(["o3_lambda", "o3_sigma"], np.array([[1, -1], [2, 1]])), + blocks=[ + TensorBlock( + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[Labels(["o3_mu"], np.arange(-1, 2).reshape(-1, 1))], + properties=Labels(["_"], np.array([[0]])), + values=np.stack( + [_l1_components_from_matrix(A[..., 0]) for A in random_rank_2_arr] + ).reshape(random_rank_2_arr.shape[0], 3, 1), + ), + TensorBlock( + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[Labels(["o3_mu"], np.arange(-2, 3).reshape(-1, 1))], + properties=Labels(["_"], np.array([[0]])), + values=np.stack( + [_l2_components_from_matrix(A[..., 0]) for A in random_rank_2_arr] + ).reshape(random_rank_2_arr.shape[0], 5, 1), + ), + ], + ) + + # Build the cartesian tensor and do cartesian to spherical + rank_2_input_cart = TensorMap( + keys=Labels.single(), + blocks=[ + TensorBlock( + values=random_rank_2_arr, + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[ + Labels([a], np.arange(3).reshape(-1, 1)) for a in ["xyz1", "xyz2"] + ], + properties=Labels(["_"], np.array([[0]])), + ) + ], + ) + rank_2_input_sph = cartesian_to_spherical( + rank_2_input_cart, ["xyz1", "xyz2"], cg_backend=cg_backend + ) + + # Extract the lambda = 1 and lambda = 2 components + l1_and_l2_input = mts.drop_blocks( + mts.remove_dimension(rank_2_input_sph, "keys", "_"), + keys=Labels(["o3_lambda"], np.array([[0]])), + ) + + assert mts.equal_metadata(l1_and_l2_input, l1_and_l2_reference) + assert mts.allclose(l1_and_l2_input, l1_and_l2_reference) + + +@pytest.mark.parametrize("cg_backend", ["python-dense", "python-sparse"]) +def test_cartesian_to_spherical_rank_3_by_equation(cg_backend): + """ + Tests cartesian_to_spherical for a random rank-2 tensor by comparing the result to + the result from *the* reference equation. + """ + # Build the reference lambda = 2 component + random_rank_3_arr = np.random.rand(100, 3, 3, 3, 1) + l3_reference = TensorMap( + keys=Labels(["o3_lambda", "o3_sigma"], np.array([[3, 1]])), + blocks=[ + TensorBlock( + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[Labels(["o3_mu"], np.arange(-3, 4).reshape(-1, 1))], + properties=Labels(["_"], np.array([[0]])), + values=np.stack( + [_l3_components_from_matrix(A[..., 0]) for A in random_rank_3_arr] + ).reshape(random_rank_3_arr.shape[0], 7, 1), + ) + ], + ) + + # Build the cartesian tensor and do cartesian to spherical + rank_3_input_cart = TensorMap( + keys=Labels.single(), + blocks=[ + TensorBlock( + values=random_rank_3_arr, + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[ + Labels([a], np.arange(3).reshape(-1, 1)) + for a in ["xyz1", "xyz2", "xyz3"] + ], + properties=Labels(["_"], np.array([[0]])), + ) + ], + ) + rank_3_input_sph = cartesian_to_spherical( + rank_3_input_cart, ["xyz1", "xyz2", "xyz3"], cg_backend=cg_backend + ) + + # Extract the lambda = 3 component + l3_input = mts.drop_blocks( + mts.remove_dimension(rank_3_input_sph, "keys", "_"), + keys=Labels(["o3_lambda"], np.array([[0], [1], [2]])), + ) + for dim in ["l_3", "k_1", "l_2", "l_1"]: + l3_input = mts.remove_dimension(l3_input, "keys", dim) + + assert mts.equal_metadata(l3_input, l3_reference) + assert mts.allclose(l3_input, l3_reference) + + +@pytest.mark.parametrize("cg_backend", ["python-dense", "python-sparse"]) +def test_cartesian_to_spherical_equivariance(cg_backend): + # Define a random cartesian rank 2 tensor + random_rank_2_arr = np.random.rand(100, 3, 3, 1) + + # Define some rotation angles, the cartesian rotation matrix, and the Wigner + # matrices + angles = np.random.randn(3) * np.pi + R = cartesian_rotation(angles) + wigner = WignerDReal(max_angular=2, angles=angles) + + # Rotate the cartesian tensor + random_rank_2_arr_rot = ( + ((random_rank_2_arr.copy().reshape(100, 3, 3) @ R.T).transpose(0, 2, 1) @ R.T) + .transpose(0, 2, 1) + .reshape(100, 3, 3, 1) + ) + + # Build the Cartesian TMs and do cartesian to spherical + # Build the cartesian tensor and do cartesian to spherical + rank_2_input_cart = TensorMap( + keys=Labels.single(), + blocks=[ + TensorBlock( + values=random_rank_2_arr, + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[ + Labels([a], np.arange(3).reshape(-1, 1)) for a in ["xyz1", "xyz2"] + ], + properties=Labels(["_"], np.array([[0]])), + ) + ], + ) + rank_2_input_sph = cartesian_to_spherical( + rank_2_input_cart, ["xyz1", "xyz2"], cg_backend=cg_backend + ) + rank_2_input_cart_rot = TensorMap( + keys=Labels.single(), + blocks=[ + TensorBlock( + values=random_rank_2_arr_rot, + samples=Labels(["system"], np.arange(100).reshape(-1, 1)), + components=[ + Labels([a], np.arange(3).reshape(-1, 1)) for a in ["xyz1", "xyz2"] + ], + properties=Labels(["_"], np.array([[0]])), + ) + ], + ) + rank_2_input_sph_rot = cartesian_to_spherical( + rank_2_input_cart_rot, ["xyz1", "xyz2"], cg_backend=cg_backend + ) + + # Extract the lambda = 2 components + l2_input = mts.drop_blocks( + mts.remove_dimension(rank_2_input_sph, "keys", "_"), + keys=Labels(["o3_lambda"], np.array([[0], [1]])), + ) + l2_input_rot = mts.drop_blocks( + mts.remove_dimension(rank_2_input_sph_rot, "keys", "_"), + keys=Labels(["o3_lambda"], np.array([[0], [1]])), + ) + + # Rotate the unrotated L=2 component + l2_input_original_rot = wigner.transform_tensormap_so3(l2_input) + + assert mts.equal_metadata(l2_input_rot, l2_input_original_rot) + assert mts.allclose(l2_input_rot, l2_input_original_rot) diff --git a/python/featomic/tests/clebsch_gordan/cg_product.py b/python/featomic/tests/clebsch_gordan/cg_product.py index 59743c7f9..33071ed29 100644 --- a/python/featomic/tests/clebsch_gordan/cg_product.py +++ b/python/featomic/tests/clebsch_gordan/cg_product.py @@ -86,13 +86,16 @@ def spherical_expansion_by_pair(frames: List[ase.Atoms]): return calculator.compute(frames) -def test_keys_are_matched(): +@pytest.mark.parametrize("cg_backend", ["python-sparse", "python-dense"]) +def test_keys_are_matched(cg_backend): """ Tests that key dimensions named the same in two tensors are matched. """ # Set up frames = h2o_isolated() - calculator = ClebschGordanProduct(max_angular=MAX_ANGULAR * 2) + calculator = ClebschGordanProduct( + max_angular=MAX_ANGULAR * 2, cg_backend=cg_backend + ) # Compute lambda-SOAP density = spherical_expansion(frames) @@ -431,3 +434,35 @@ def test_device_dtype(dtype, device): o3_lambda_1_new_name="l_1", o3_lambda_2_new_name="l_2", ) + + +def test_dense_sparse_agree(): + """ + Tests that the max_angular is sufficient to correlate the two tensors when not large + enough to cover MAX_ANGULAR, but when an angular cutoff is applied. + """ + frames = h2o_isolated() + density = spherical_expansion(frames) + + results = [] + for cg_backend in ["python-sparse", "python-dense"]: + # max_angular to be twice as big here if not using an angular cutoff + calculator = ClebschGordanProduct( + max_angular=MAX_ANGULAR, + cg_backend=cg_backend, + ) + results.append( + calculator.compute( + metatensor.rename_dimension(density, "properties", "n", "n_1"), + metatensor.rename_dimension(density, "properties", "n", "n_2"), + o3_lambda_1_new_name="l_1", + o3_lambda_2_new_name="l_2", + selected_keys=Labels( + names=["o3_lambda"], + values=np.arange(MAX_ANGULAR + 1).reshape(-1, 1), + ), + ) + ) + + assert metatensor.equal_metadata(results[0], results[1]) + assert metatensor.allclose(results[0], results[1]) diff --git a/python/featomic/tests/clebsch_gordan/coefficients.py b/python/featomic/tests/clebsch_gordan/coefficients.py index 3c36a6d19..44fa37249 100644 --- a/python/featomic/tests/clebsch_gordan/coefficients.py +++ b/python/featomic/tests/clebsch_gordan/coefficients.py @@ -29,9 +29,10 @@ def complex_to_real_matrix(sph): matrix = _complex2real(ell, sph) - real = sph @ matrix + real = matrix @ sph.T + real = real.T - assert np.linalg.norm(np.imag(real)) < 1e-15 + assert np.allclose(np.imag(real), 0) return np.real(real) diff --git a/python/featomic/tests/clebsch_gordan/rotations.py b/python/featomic/tests/clebsch_gordan/rotations.py index db67bc79d..ba0cc83f2 100644 --- a/python/featomic/tests/clebsch_gordan/rotations.py +++ b/python/featomic/tests/clebsch_gordan/rotations.py @@ -96,9 +96,12 @@ def __init__(self, max_angular: int, angles: Sequence[float] = None): :param angles: Sequence[float], the alpha, beta, gamma Euler angles, in radians. """ self.max_angular = max_angular - # Randomly generate Euler angles between 0 and 2 pi if none are provided + # Randomly generate Euler angles between if none are provided if angles is None: - angles = np.random.uniform(size=(3)) * 2 * np.pi + psi = 2 * np.pi * np.random.uniform() # psi is between 0 and 2pi + theta = np.pi * np.random.uniform() # theta is between 0 and pi + phi = 2 * np.pi * np.random.uniform() # phi is between 0 and 2pi + angles = np.array([psi, theta, phi]) self.angles = angles self.rotation = cartesian_rotation(angles)