_radii

Vectorised sphere-cast root finders.

Given a centre origin and a batch of unit directions, find the radial distance r at which a scalar field along each ray equals a target isovalue. The shape-descriptor pipeline samples a Gauss-Legendre / equispaced angular grid of directions and uses the resulting radii as the input to a spherical-harmonic transform.

Both the promolecule density and the Hirshfeld stockholder weight are monotonically decreasing along outward rays in the regime of interest, so the unique root in [lower, upper] is found via:

1. one batched evaluation of the field at K points along each ray to locate a tight bracket where the sign flips (K = 64 by default);
2. refine_iter iterations of bisection inside that tight bracket.

Step 1 amortises the rho-evaluation overhead across all (ray, sample) pairs in a single numpy call; step 2 then bisects within the ~1.3-Å initial bracket. The defaults (n_samples=16, max_iter=10) give ~6e-4 Å precision per ray on a 20-Å bracket, well below typical surface mesh accuracy, and run ~3× faster than plain bisection over the full bracket.

sphere_promolecule_radii(promol, origin, directions, lower, upper, tol=1e-07, max_iter=10, isovalue=0.0002, n_samples=16)

Find the per-ray distance at which the promolecule density equals isovalue.

Parameters:

Name	Type	Description	Default
promol	PromoleculeDensity	A PromoleculeDensity (the python wrapper).	required
origin		(3,) array, the ray start point.	required
directions		(N, 3) array of ray directions (typically unit vectors).	required
lower, upper		bracket the root in radial distance.	required
tol	float	kept for parity with the legacy signature; ignored.	1e-07
max_iter	int	number of bisection iterations applied within the tight bracket found by the coarse search. Default 10 gives an absolute tolerance of (upper - lower) / (n_samples - 1) / 2**max_iter, ~1e-3 Å for the default 20-Å bracket and 16 samples.	10
isovalue	float	target density.	0.0002
n_samples	int	number of coarse samples per ray. Larger values find a tighter initial bracket at the cost of one bigger batched eval.	16

Returns:

Type	Description
	(N,) array of radial distances. Rays whose coarse search finds
	no sign change return -1.0.

Source code in chmpy/interpolate/_radii.py

def sphere_promolecule_radii(
    promol: PromoleculeDensity,
    origin,
    directions,
    lower: float,
    upper: float,
    tol: float = 1e-7,
    max_iter: int = 10,
    isovalue: float = 0.0002,
    n_samples: int = 16,
):
    """Find the per-ray distance at which the promolecule density equals ``isovalue``.

    Args:
        promol: A ``PromoleculeDensity`` (the python wrapper).
        origin: ``(3,)`` array, the ray start point.
        directions: ``(N, 3)`` array of ray directions (typically unit vectors).
        lower, upper: bracket the root in radial distance.
        tol: kept for parity with the legacy signature; ignored.
        max_iter: number of bisection iterations applied within the tight
            bracket found by the coarse search. Default 10 gives an absolute
            tolerance of ``(upper - lower) / (n_samples - 1) / 2**max_iter``,
            ~1e-3 Å for the default 20-Å bracket and 16 samples.
        isovalue: target density.
        n_samples: number of coarse samples per ray. Larger values find a
            tighter initial bracket at the cost of one bigger batched eval.

    Returns:
        ``(N,)`` array of radial distances. Rays whose coarse search finds
        no sign change return ``-1.0``.
    """
    del tol
    origin = np.asarray(origin, dtype=np.float32)
    directions = np.asarray(directions, dtype=np.float32)
    iso = np.float32(isovalue)
    positions = promol.positions
    rho_data = promol.rho_data

    def eval_pts(pts):
        return _backends.rho(positions, rho_data, _DOMAIN, pts) - iso

    return _root_along_rays(
        eval_pts, origin, directions, lower, upper, n_samples, max_iter
    )

sphere_stockholder_radii(stock, origin, directions, lower, upper, tol=1e-07, max_iter=10, isovalue=0.5, n_samples=16)

Find the per-ray distance at which the stockholder weight equals isovalue.

Source code in chmpy/interpolate/_radii.py

def sphere_stockholder_radii(
    stock: StockholderWeight,
    origin,
    directions,
    lower: float,
    upper: float,
    tol: float = 1e-7,
    max_iter: int = 10,
    isovalue: float = 0.5,
    n_samples: int = 16,
):
    """Find the per-ray distance at which the stockholder weight equals ``isovalue``."""
    del tol
    origin = np.asarray(origin, dtype=np.float32)
    directions = np.asarray(directions, dtype=np.float32)
    iso = np.float32(isovalue)
    pos_a = stock.dens_a.positions
    rho_a = stock.dens_a.rho_data
    pos_b = stock.dens_b.positions
    rho_b = stock.dens_b.rho_data
    bg = stock.background

    def eval_pts(pts):
        return _backends.weights(pos_a, rho_a, pos_b, rho_b, _DOMAIN, pts, bg) - iso

    return _root_along_rays(
        eval_pts, origin, directions, lower, upper, n_samples, max_iter
    )