Projected Radial Distribution Function¶

The projected RDF decomposes the standard pair correlation \(g(r)\) into orientationally-resolved components using spherical harmonics. It answers the question: do pairs at distance \(r\) have a preferred direction, and in which direction?

Theory¶

Expansion onto spherical harmonics¶

The full angular pair correlation can be written as:

\[ g(\mathbf{r}) = \sum_{l,m} g_{lm}(r) \, Y_l^m(\hat{\mathbf{r}}) \]

For amorphous materials the \(l=0\) term recovers the standard isotropic \(g(r)\). The \(l=2\) terms carry the lowest-order anisotropy signal and are the ones accessible from glasses under deformation. The spherical harmonic decomposition of the pair correlation has been used to quantify shear- and deformation-induced structural anisotropy in model glass formers [1] and oxide glasses [2,3].

\(l=2\) components¶

Using real spherical harmonics and Cartesian unit-vector components \(\hat{u}_x, \hat{u}_y, \hat{u}_z = \mathbf{r}_{ij}/r_{ij}\):

Component	Cartesian form	Physical coupling
\(Y_2^0\)	\(\sqrt{\tfrac{5}{4\pi}}\!\left(\tfrac{3}{2}\hat{u}_z^2 - \tfrac{1}{2}\right)\)	Uniaxial anisotropy along \(z\)
\(\text{Re}\,Y_2^1\)	\(\sqrt{\tfrac{15}{8\pi}}\,\hat{u}_x\hat{u}_z\)	Shear coupling in \(xz\)
\(\text{Im}\,Y_2^1\)	\(\sqrt{\tfrac{15}{8\pi}}\,\hat{u}_y\hat{u}_z\)	Shear coupling in \(yz\)
\(\text{Re}\,Y_2^2\)	\(\sqrt{\tfrac{15}{32\pi}}\!\left(\hat{u}_x^2 - \hat{u}_y^2\right)\)	Uniaxial anisotropy in \(xy\)
\(\text{Im}\,Y_2^2\)	\(\sqrt{\tfrac{15}{8\pi}}\,\hat{u}_x\hat{u}_y\)	Shear coupling in \(xy\)

Uniaxial signal¶

\(g_{20}(r)\) is accumulated with \(z\) as the internal reference axis. For deformation along a different axis the signal is reconstructed analytically:

\[ g_\text{uniaxial}(r) = \begin{cases} g_{20}(r) & \text{axis} = z \\ g_{20}(r) + \sqrt{3/5}\,\text{Re}\,g_{22}(r) & \text{axis} = x \\ g_{20}(r) - \sqrt{3/5}\,\text{Re}\,g_{22}(r) & \text{axis} = y \end{cases} \]

Interpretation:

\(g_\text{uniaxial}(r) > 0\) at the first peak: pairs at that bond length preferentially align along the deformation axis.
\(g_\text{uniaxial}(r) < 0\): pairs preferentially align perpendicular to the axis (transverse).
\(g_\text{uniaxial}(r) \approx 0\) everywhere: the structure is isotropic at that length scale.

Shear signal¶

Each shear plane maps to the Y_lm component whose Cartesian form contains the off-diagonal product of that plane's two axes:

\[ g_\text{shear}(r) = \begin{cases} \text{Im}\,g_{22}(r) & \text{plane} = xy \\ \text{Re}\,g_{21}(r) & \text{plane} = xz \\ \text{Im}\,g_{21}(r) & \text{plane} = yz \end{cases} \]

A nonzero value at the first-shell distance means pairs at that distance develop a net preference for the shear-plane direction.

Computing the projected RDF¶

from amorphouspy import compute_projected_rdf

# Uniaxial only (z-compression)
r, uniaxial, _ = compute_projected_rdf(
    structure,
    deformation_axis="z",
    r_max=8.0,
    n_bins=500,
)

# Shear only (xy plane)
r, _, shear = compute_projected_rdf(
    structure,
    shear_plane="xy",
    r_max=8.0,
    n_bins=500,
)

# Both at once
r, uniaxial, shear = compute_projected_rdf(
    structure,
    deformation_axis="z",
    shear_plane="xy",
    r_max=8.0,
    n_bins=500,
)

Interactive tutorial¶

See the Projected RDF Tutorial notebook for:

Interactive 3D visualisations of the \(Y_2^0\) and \(Y_2^2\) spherical harmonic lobes with simulation-box overlays
Worked example on a strained Si supercell showing a nonzero uniaxial signal
Side-by-side comparison of isotropic vs deformed \(g_\text{uniaxial}(r)\)

API reference¶

`amorphouspy.properties.structural.projected_rdf.compute_projected_rdf(frames: Atoms | list[Atoms], deformation_axis: Literal['x', 'y', 'z'] | None = None, shear_plane: Literal['xy', 'xz', 'yz'] | None = None, r_max: float = 10.0, n_bins: int = 500, type_pairs: list[tuple[int, int]] | None = None) -> tuple[np.ndarray, dict[tuple[int, int], np.ndarray] | None, dict[tuple[int, int], np.ndarray] | None]` ¶

Compute the uniaxial and shear projected radial distribution functions.

Expands the pair correlation function onto the l=2 spherical harmonic components that carry the deformation signal:

uniaxial_rdf: the Y20-based component along deformation_axis. Vanishes for isotropic structures; grows with uniaxial strain.
shear_rdf: the Y21 or Y22 component for shear_plane. Vanishes for isotropic structures; grows with shear.

At least one of deformation_axis or shear_plane must be provided. Omit the one you do not need, its slot in the return tuple will be None.

For multiple frames the per-frame normalised values are averaged, so NPT trajectories with varying cell geometries are handled correctly.

Performance: Numba-JIT kernels accumulate all five needed Y*_lm components in a single i<j pass. On first call Numba compiles the kernels (~5 s); subsequent calls use cached bytecode.

Parameters:

Name	Type	Description	Default
`frames`	`Atoms \| list[Atoms]`	Single ASE Atoms object or list of frames for trajectory averaging.	required
`deformation_axis`	`Literal['x', 'y', 'z'] \| None`	Axis along which uniaxial deformation was applied ('x', 'y', or 'z'). Determines which linear combination of Y20 and Re Y22 is returned as the uniaxial signal. `None` skips uniaxial extraction and returns `None` for `uniaxial_rdf`.	`None`
`shear_plane`	`Literal['xy', 'xz', 'yz'] \| None`	Plane in which shear deformation was applied ('xy', 'xz', or 'yz'). Determines which Y21/Y22 component is returned as the shear signal. `None` skips shear extraction and returns `None` for `shear_rdf`.	`None`
`r_max`	`float`	Maximum distance in Angstrom (default 10.0). Automatically clamped to floor(min_perpendicular_height / 2) with a warning if exceeded.	`10.0`
`n_bins`	`int`	Number of radial bins (default 500).	`500`
`type_pairs`	`list[tuple[int, int]] \| None`	List of (atomic_number_a, atomic_number_b) pairs to compute. None computes all unique unordered combinations present in the first frame.	`None`

Returns:

Name	Type	Description
`r`	`ndarray`	Radial bin centres in Angstrom, shape (n_bins,).
`uniaxial_rdf`	`dict[tuple[int, int], ndarray] \| None`	Dict keyed by (type1, type2) -> real array of shape (n_bins,), or `None` if `deformation_axis` was not provided.
`shear_rdf`	`dict[tuple[int, int], ndarray] \| None`	Dict keyed by (type1, type2) -> real array of shape (n_bins,), or `None` if `shear_plane` was not provided.

Raises:

Type	Description
`ValueError`	If both `deformation_axis` and `shear_plane` are `None`.
`ValueError`	If r_max after clamping would be <= 0.

Example

from ase.build import bulk atoms = bulk("Si", "diamond", a=5.43) * (4, 4, 4)

Uniaxial only -- no shear plane needed¶

r, uniaxial, _ = compute_projected_rdf( ... atoms, deformation_axis="z", r_max=6.0 ... ) signal_SiSi = uniaxial[(14, 14)] # near zero for undeformed bulk

Source code in amorphouspy/src/amorphouspy/properties/structural/projected_rdf.py

def compute_projected_rdf(
    frames: Atoms | list[Atoms],
    deformation_axis: Literal["x", "y", "z"] | None = None,
    shear_plane: Literal["xy", "xz", "yz"] | None = None,
    r_max: float = 10.0,
    n_bins: int = 500,
    type_pairs: list[tuple[int, int]] | None = None,
) -> tuple[
    np.ndarray,
    dict[tuple[int, int], np.ndarray] | None,
    dict[tuple[int, int], np.ndarray] | None,
]:
    r"""Compute the uniaxial and shear projected radial distribution functions.

    Expands the pair correlation function onto the l=2 spherical harmonic
    components that carry the deformation signal:

    - **uniaxial_rdf**: the Y20-based component along ``deformation_axis``.
      Vanishes for isotropic structures; grows with uniaxial strain.
    - **shear_rdf**: the Y21 or Y22 component for ``shear_plane``.
      Vanishes for isotropic structures; grows with shear.

    At least one of ``deformation_axis`` or ``shear_plane`` must be provided.
    Omit the one you do not need, its slot in the return tuple will be ``None``.

    For multiple frames the per-frame normalised values are averaged, so NPT
    trajectories with varying cell geometries are handled correctly.

    Performance: Numba-JIT kernels accumulate all five needed Y*_lm components
    in a single i<j pass. On first call Numba compiles the kernels (~5 s);
    subsequent calls use cached bytecode.

    Args:
        frames: Single ASE Atoms object or list of frames for trajectory averaging.
        deformation_axis: Axis along which uniaxial deformation was applied
            ('x', 'y', or 'z'). Determines which linear combination of Y20
            and Re Y22 is returned as the uniaxial signal. ``None`` skips
            uniaxial extraction and returns ``None`` for ``uniaxial_rdf``.
        shear_plane: Plane in which shear deformation was applied
            ('xy', 'xz', or 'yz'). Determines which Y21/Y22 component is
            returned as the shear signal. ``None`` skips shear extraction and
            returns ``None`` for ``shear_rdf``.
        r_max: Maximum distance in Angstrom (default 10.0). Automatically clamped to
            floor(min_perpendicular_height / 2) with a warning if exceeded.
        n_bins: Number of radial bins (default 500).
        type_pairs: List of (atomic_number_a, atomic_number_b) pairs to compute.
            None computes all unique unordered combinations present in the first frame.

    Returns:
        r: Radial bin centres in Angstrom, shape (n_bins,).
        uniaxial_rdf: Dict keyed by (type1, type2) -> real array of shape (n_bins,),
            or ``None`` if ``deformation_axis`` was not provided.
        shear_rdf: Dict keyed by (type1, type2) -> real array of shape (n_bins,),
            or ``None`` if ``shear_plane`` was not provided.

    Raises:
        ValueError: If both ``deformation_axis`` and ``shear_plane`` are ``None``.
        ValueError: If r_max after clamping would be <= 0.

    Example:
        >>> from ase.build import bulk
        >>> atoms = bulk("Si", "diamond", a=5.43) * (4, 4, 4)
        >>> # Uniaxial only -- no shear plane needed
        >>> r, uniaxial, _ = compute_projected_rdf(
        ...     atoms, deformation_axis="z", r_max=6.0
        ... )
        >>> signal_SiSi = uniaxial[(14, 14)]   # near zero for undeformed bulk

    """
    if deformation_axis is None and shear_plane is None:
        msg = "At least one of deformation_axis or shear_plane must be provided."
        raise ValueError(msg)

    if isinstance(frames, Atoms):
        frames = [frames]

    r_max = _resolve_r_max(frames[0], r_max)
    canonical_pairs, pair_lut = _build_pair_lut(frames[0], type_pairs)

    n_pairs = len(canonical_pairs)
    bin_edges = np.linspace(0.0, r_max, n_bins + 1)
    bin_centers = 0.5 * (bin_edges[1:] + bin_edges[:-1])
    bin_width = float(bin_edges[1] - bin_edges[0])
    shell_volumes = 4.0 * np.pi * bin_centers**2 * bin_width

    accumulated_Y20 = np.zeros((n_pairs, n_bins), dtype=np.float64)
    accumulated_Y21_real = np.zeros((n_pairs, n_bins), dtype=np.float64)
    accumulated_Y21_imag = np.zeros((n_pairs, n_bins), dtype=np.float64)
    accumulated_Y22_real = np.zeros((n_pairs, n_bins), dtype=np.float64)
    accumulated_Y22_imag = np.zeros((n_pairs, n_bins), dtype=np.float64)

    for atoms in frames:
        _accumulate_frame(
            atoms,
            canonical_pairs,
            pair_lut,
            r_max,
            bin_width,
            n_bins,
            shell_volumes,
            accumulated_Y20,
            accumulated_Y21_real,
            accumulated_Y21_imag,
            accumulated_Y22_real,
            accumulated_Y22_imag,
        )

    n_frames = len(frames)
    uniaxial_rdf: dict[tuple[int, int], np.ndarray] | None = {} if deformation_axis is not None else None
    shear_rdf: dict[tuple[int, int], np.ndarray] | None = {} if shear_plane is not None else None

    for pair_idx, (type1, type2) in enumerate(canonical_pairs):
        g20 = accumulated_Y20[pair_idx] / n_frames
        g21_real = accumulated_Y21_real[pair_idx] / n_frames
        g21_imag = accumulated_Y21_imag[pair_idx] / n_frames
        g22_real = accumulated_Y22_real[pair_idx] / n_frames
        g22_imag = accumulated_Y22_imag[pair_idx] / n_frames

        if deformation_axis is not None and uniaxial_rdf is not None:
            uniaxial_rdf[(type1, type2)] = _extract_uniaxial_signal(g20, g22_real, deformation_axis)
        if shear_plane is not None and shear_rdf is not None:
            shear_rdf[(type1, type2)] = _extract_shear_signal(g21_real, g21_imag, g22_imag, shear_plane)

    return bin_centers, uniaxial_rdf, shear_rdf

References¶

J. Zausch and J. Horbach, "The build-up and relaxation of stresses in a glass-forming soft-sphere mixture under shear: A computer simulation study," EPL 88, 60001 (2009). https://doi.org/10.1209/0295-5075/88/60001
S. Ganisetti, A. Atila, J. Guénolé, A. Prakash, J. Horbach, L. Wondraczek, and E. Bitzek, "The origin of deformation induced topological anisotropy in silica glass," Acta Mater. 257, 119108 (2023). https://doi.org/10.1016/j.actamat.2023.119108
A. Atila and E. Bitzek, "Atomistic origins of deformation-induced structural anisotropy in metaphosphate glasses and its influence on mechanical properties," J. Non-Cryst. Solids 627, 122822 (2024). https://doi.org/10.1016/j.jnoncrysol.2024.122822