pyextal.jaxMultipole.harmonics
==============================

.. py:module:: pyextal.jaxMultipole.harmonics

.. autoapi-nested-parse::

   Density-normalised real spherical harmonics and their site-symmetry constraints.

   Convention
   ----------

   The angular functions are XD's (Hansen & Coppens 1978; Paturle & Coppens 1988):

   .. math::

       c_{lm\pm}(\hat u) = \frac{d^m P_l}{dz^m}(z)\;
           \{\mathrm{Re}, \mathrm{Im}\}\,(x + iy)^m, \qquad
       d_{lm\pm} = c_{lm\pm} / M_{lm},

   i.e. associated Legendre functions *without* the Condon-Shortley phase times
   ``cos m phi`` (``+``) or ``sin m phi`` (``-``), written as polynomials in the
   Cartesian components of a unit vector. ``M_lm`` makes them density-normalised:
   ``int |d_lm| dOmega = 2`` for ``l > 0`` and ``d_00 = 1/(4 pi)``. A population
   ``P_lm`` is then roughly the number of electrons moved from the negative lobes
   to the positive ones, and the cubic combination comes out as the familiar
   ``P44+ = 0.74045 P40``.

   ``M_lm`` is computed once, from 1-D integrals split at the roots of the
   Legendre polynomial, so it carries no grid error.

   Coefficients are ordered ``l = 0..lmax``, within ``l`` as ``m = 0, 1+, 1-, 2+,
   2-, ...``, and labelled ``P00, P10, P11+, P11-, P20, ...``.

   Symmetry
   --------

   A rotation (proper or improper) ``R`` of the density, ``rho'(r) = rho(R^T r)``,
   mixes the coefficients of one ``l`` by a ``(2l+1)`` square matrix ``Q(R)``
   (:func:`real_wigner`). It is obtained by least squares from the functions
   evaluated at random directions rather than from closed-form Wigner matrices, so
   improper operations are handled without special cases (inversion gives
   ``(-1)^l``). A site symmetry ``S`` allows exactly the populations with ``Q(R) P =
   P`` for every ``R`` in ``S``, and :func:`allowed_multipole_basis` returns the
   null space of the stacked ``I - Q(R)``, canonicalised to reduced row-echelon
   form in the order above. In a frame chosen as Kurki-Suonio's index-picking rules
   choose it (z along the main axis, x along a 2-fold axis or mirror normal) that
   basis is the usual list of individual ``P_lm``, plus ties like the cubic one.



Functions
---------

.. autoapisummary::

   pyextal.jaxMultipole.harmonics.n_functions
   pyextal.jaxMultipole.harmonics.labels
   pyextal.jaxMultipole.harmonics.index
   pyextal.jaxMultipole.harmonics.normalisation
   pyextal.jaxMultipole.harmonics.dlm
   pyextal.jaxMultipole.harmonics.real_wigner
   pyextal.jaxMultipole.harmonics.allowed_multipole_basis
   pyextal.jaxMultipole.harmonics.rotation_angle
   pyextal.jaxMultipole.harmonics.character_count


Module Contents
---------------

.. py:function:: n_functions(lmax)

.. py:function:: labels(lmax)

   ``('P00', 'P10', 'P11+', 'P11-', ...)`` up to ``lmax``.


.. py:function:: index(l, m=0, sign='+')

   Position of ``P_lm+/-`` in the flat coefficient vector.


.. py:function:: normalisation(l, m)

   ``M_lm``, the divisor that density-normalises ``c_lm``.


.. py:function:: dlm(u, lmax=4, xp=jnp)

   Density-normalised real harmonics at directions ``u``.

   :param u: ``(..., 3)`` vectors in the local Cartesian frame; they
             need not be unit length. The zero vector gives ``d_00`` and zeros,
             which is harmless because every ``l > 0`` radial transform vanishes
             at ``k = 0``.
   :type u: array-like
   :param lmax: Highest order.
   :type lmax: int
   :param xp: ``jax.numpy`` (default, differentiable) or ``numpy`` --
              for constants (symmetry matrices, precomputed geometry), where
              NumPy is exact float64 and avoids dispatching hundreds of tiny
              operations to an accelerator.
   :type xp: module

   :returns: ``(..., (lmax+1)^2)`` in :func:`labels` order.
   :rtype: array


.. py:function:: real_wigner(R, l)

   The matrix ``Q`` with ``P' = Q P`` for the density rotated by ``R``.

   :param R: (3, 3) orthogonal matrix, proper or improper, acting
             on Cartesian vectors in the frame the harmonics are defined in.
   :type R: np.ndarray
   :param l: Order.
   :type l: int

   :returns: ``(2l+1, 2l+1)``, orthogonal.
   :rtype: np.ndarray


.. py:function:: allowed_multipole_basis(rots, l, tol=1e-08)

   Populations of order ``l`` a site symmetry leaves free.

   :param rots: ``(m, 3, 3)`` site-symmetry operations as Cartesian
                matrices in the atom's *local* frame.
   :type rots: np.ndarray
   :param l: Order.
   :type l: int
   :param tol: Singular-value cutoff.
   :type tol: float

   :returns: ``(2l+1, nfree)``. Columns are in reduced row-echelon form:
             each starts with a 1 at its own leading coefficient, and a tie shows up
             as further non-zero entries (``[1, 0, ..., 0.74045, 0, ...]`` for the
             cubic l = 4 function).
   :rtype: np.ndarray


.. py:function:: rotation_angle(R)

   ``(angle, proper)`` of an orthogonal matrix; improper ones as ``-R``.


.. py:function:: character_count(rots, l)

   Number of invariant functions of order ``l``, from group characters.

   ``(1/|G|) sum chi_l(R)`` with ``chi_l = sin((l+1/2) phi) / sin(phi/2)``,
   times ``(-1)^l`` for improper operations. Independent of the projector
   route, so it serves as a check of :func:`allowed_multipole_basis`.