Spherical coordinates and/or complex harmonics

The algorithms implemented in sphericart are designed to work with Cartesian input positions and real spherical (or solid) harmonics. However, depending on the use case, it might be more convenient to compute harmonics from spherical coordinates and/or work with complex harmonics.

Below, we provide a series of Python examples that illustrate how to use the sphericart library in such cases. The examples can be easily adapted to the other languages supported by the library. All the examples are consistent with the definitions of the spherical harmonics from Wikipedia.

The use of these (and similar) adaptors is not recommended for applications where performance is critical, as they introduce some computational overhead.

Computing harmonics from spherical coordinates

This is a simple class that computes the spherical harmonics from spherical coordinates.

import numpy as np
from sphericart import SphericalHarmonics as CartesianSphericalHarmonics


class SphericalHarmonics:
    def __init__(self, l_max):
        self.l_max = l_max
        self.cartesian_spherical_harmonics = CartesianSphericalHarmonics(l_max, normalized=True)

    def compute(self, theta, phi):
        assert theta.shape == phi.shape
        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)
        xyz = np.stack([x, y, z], axis=-1)
        return self.cartesian_spherical_harmonics.compute(xyz)

Computing complex harmonics

This is a simple class that computes complex spherical harmonics.

import numpy as np
from sphericart import SphericalHarmonics as RealSphericalHarmonics


class SphericalHarmonics:
    def __init__(self, l_max):
        self.l_max = l_max
        self.real_spherical_harmonics = RealSphericalHarmonics(l_max, normalized=True)

    def compute(self, xyz):
        if xyz.dtype == np.float32:
            complex_dtype = np.complex64
        elif xyz.dtype == np.float64:
            complex_dtype = np.complex128
        
        real_spherical_harmonics = self.real_spherical_harmonics.compute(xyz)

        one_over_sqrt_2 = 1.0 / np.sqrt(2.0)

        complex_spherical_harmonics = np.zeros((self.l_max + 1) ** 2, dtype=complex_dtype)
        for l in range(self.l_max + 1):
            for m in range(-l, l + 1):
                l_m_index = l ** 2 + l + m
                l_minus_m_index = l ** 2 + l - m
                if m < 0:
                    complex_spherical_harmonics[..., l_m_index] = (
                        one_over_sqrt_2 * (real_spherical_harmonics[..., l_minus_m_index] - 1j * real_spherical_harmonics[..., l_m_index])
                    )
                elif m == 0:
                    complex_spherical_harmonics[..., l_m_index] = real_spherical_harmonics[..., l_m_index]
                else:  # m > 0
                    complex_spherical_harmonics[..., l_m_index] = (
                        (-1)**m * one_over_sqrt_2 * (real_spherical_harmonics[..., l_m_index] + 1j * real_spherical_harmonics[..., l_minus_m_index])
                    )

Computing complex harmonics from spherical coordinates

This is a simple class that computes complex spherical harmonics from spherical coordinates. Its correctness can be verified against the scipy implementation of the spherical harmonics.

import numpy as np
from sphericart import SphericalHarmonics as RealCartesianSphericalHarmonics


class SphericalHarmonics:
    def __init__(self, l_max):
        self.l_max = l_max
        self.cartesian_spherical_harmonics = RealCartesianSphericalHarmonics(l_max, normalized=True)

    def compute(self, theta, phi):
        assert theta.shape == phi.shape

        x = np.sin(theta) * np.cos(phi)
        y = np.sin(theta) * np.sin(phi)
        z = np.cos(theta)
        xyz = np.stack([x, y, z], axis=-1)
        real_spherical_harmonics = self.cartesian_spherical_harmonics.compute(xyz)

        one_over_sqrt_2 = 1.0 / np.sqrt(2.0)

        if xyz.dtype == np.float32:
            complex_dtype = np.complex64
        elif xyz.dtype == np.float64:
            complex_dtype = np.complex128
        else:
            raise ValueError("Unsupported dtype")
        
        complex_spherical_harmonics = np.zeros(theta.shape + ((self.l_max + 1) ** 2,), dtype=complex_dtype)

        for l in range(self.l_max + 1):
            for m in range(-l, l + 1):
                l_m_index = l ** 2 + l + m
                l_minus_m_index = l ** 2 + l - m
                if m < 0:
                    complex_spherical_harmonics[..., l_m_index] = (
                        one_over_sqrt_2 * (real_spherical_harmonics[..., l_minus_m_index] - 1j * real_spherical_harmonics[..., l_m_index])
                    )
                elif m == 0:
                    complex_spherical_harmonics[..., l_m_index] = real_spherical_harmonics[..., l_m_index]
                else:  # m > 0
                    complex_spherical_harmonics[..., l_m_index] = (
                        (-1)**m * one_over_sqrt_2 * (real_spherical_harmonics[..., l_m_index] + 1j * real_spherical_harmonics[..., l_minus_m_index])
                    )

        return complex_spherical_harmonics