Choosing a practical neutral-proatom interatomic-surface (IAS) proxy¶

This notebook documents the numerical study used to select practical parameters for the IAS estimators in atomref. It compares two pairwise quantities derived from the packaged neutral spherical proatomic densities:

1. a stable proatomic boundary, based on equal proatomic contributions or a midpoint between fixed low-density contours;
2. an optional promolecular line-density minimum, searched only where both proatomic contributions remain above the cutoff and only to a declared practical spatial resolution.

The calculations below form a standalone numerical sensitivity study rather than the package API and do not claim to locate an exact QTAIM zero-flux surface. Local helper functions isolate comparisons behind the documented cutoff and resolution policy.

Prerequisites. This supporting study assumes familiarity with radial electron densities, one-dimensional minimization, and the distinction between neutral-proatom models and molecular-density QTAIM analysis.

Numerical policy examined here¶

The per-atom neutral-tail cutoff is

$$ \rho_c = 10^{-4}\ \mathrm{electron/bohr^3}. $$

The default boundary mode uses homonuclear symmetry, equal proatomic contributions in the overlap region, and the midpoint between cutoff contours when a low-density gap opens. The optional minimum mode searches only the interval where both contributions are at least $\rho_c$, and coalesces features below 0.01 bohr (about 0.0053 Å).

Scientific context and references¶

The two modes answer different questions.

• QTAIM defines atomic basins through zero-flux surfaces of the molecular electron density; a critical point on an internuclear path is therefore a molecular-density object, not something that neutral spherical proatoms can reproduce exactly. See R. F. W. Bader, Chemical Reviews 91 (1991), 893–928, doi:10.1021/cr00005a013.
• Hirshfeld's stockholder construction assigns molecular density in proportion to free-atom reference densities. In the present two-proatom line model, equal neutral-proatom densities therefore give equal pairwise stockholder contributions. See F. L. Hirshfeld, Theoretica Chimica Acta 44 (1977), 129–138, doi:10.1007/BF00549096.
• Radial-density constructions can be useful atom/bond models, but they remain model definitions that need comparison with molecular-density results. See P. L. Warburton, R. A. Poirier, and D. Nippard, J. Phys. Chem. A 115 (2011), 852–867, doi:10.1021/jp1093417.

These references motivate the interpretation. The numerical experiments below support only the documented numerical policy for the supplied neutral-proatom data; they do not validate either mode as an exact QTAIM surface.

In [1]:

from __future__ import annotations

from bisect import bisect_right
from dataclasses import dataclass
from functools import lru_cache
import math
import random
import time

from pathlib import Path
import sys

ROOT = Path.cwd().resolve()
if not (ROOT / "src" / "atomref").is_dir():
    for candidate in (ROOT.parent, ROOT.parent.parent):
        if (candidate / "src" / "atomref").is_dir():
            ROOT = candidate
            break
if not (ROOT / "src" / "atomref").is_dir():
    raise RuntimeError("Run this notebook from the atomref repository checkout")
sys.path.insert(0, str(ROOT / "src"))

import atomref as ar

try:
    import matplotlib.pyplot as plt
except ImportError:  # The package itself has no plotting dependency.
    plt = None

ASSET_DIR = ROOT / "docs" / "assets" / "ias-method-study"
ASSET_DIR.mkdir(parents=True, exist_ok=True)

TAIL_CUTOFF = 1.0e-4
PRACTICAL_RESOLUTION = 0.01
INITIAL_STEP = 0.02
CONFIRM_STEP = 0.01
FALLBACK_STEP = 0.005
COMPETITIVE_RELATIVE_GAP = 1.0e-4


@dataclass(frozen=True)
class StudyEstimate:
    method: str
    status: str
    position: float | None
    rho_a: float | None
    rho_b: float | None
    rho_sum: float | None
    contour_separation: float
    search_resolution: float | None = None
    search_converged: bool | None = None
    alternative_position: float | None = None
    relative_depth_gap: float | None = None
    probes: int = 0


@lru_cache(maxsize=None)
def profile(symbol: str):
    result = ar.get_proatomic_density_profile(symbol)
    if result is None:
        raise ValueError(f"missing profile: {symbol}")
    return result


@lru_cache(maxsize=None)
def cutoff_radius(symbol: str) -> float:
    p = profile(symbol)
    right = next(i for i, value in enumerate(p.densities) if value <= TAIL_CUTOFF)
    if p.densities[right] == TAIL_CUTOFF:
        return p.radii[right]
    left = right - 1
    lr0 = math.log(p.radii[left])
    lr1 = math.log(p.radii[right])
    ld0 = math.log(p.densities[left])
    ld1 = math.log(p.densities[right])
    fraction = (math.log(TAIL_CUTOFF) - ld0) / (ld1 - ld0)
    return math.exp(lr0 + fraction * (lr1 - lr0))


def rho(symbol: str, radius: float) -> float:
    return profile(symbol)._evaluate_bohr(radius)


def objective(atom_a: str, atom_b: str, x: float, distance: float) -> float:
    return rho(atom_a, x) + rho(atom_b, distance - x)


def continuous_log_density(symbol: str, radius: float) -> float:
    p = profile(symbol)
    if radius <= p.radii[0]:
        return math.log(p.densities[0])
    index = min(bisect_right(p.radii, radius) - 1, len(p.radii) - 2)
    lr0 = math.log(p.radii[index])
    lr1 = math.log(p.radii[index + 1])
    ld0 = math.log(p.densities[index])
    ld1 = math.log(p.densities[index + 1])
    fraction = (math.log(radius) - lr0) / (lr1 - lr0)
    return ld0 + fraction * (ld1 - ld0)


def equal_contribution_position(atom_a: str, atom_b: str, distance: float) -> float | None:
    if atom_a == atom_b:
        return distance / 2.0

    def difference(x: float) -> float:
        return continuous_log_density(atom_a, x) - continuous_log_density(
            atom_b, distance - x
        )

    left_value = difference(0.0)
    right_value = difference(distance)
    if left_value == 0.0:
        return 0.0
    if right_value == 0.0:
        return distance
    if left_value * right_value > 0.0:
        return None

    left = 0.0
    right = distance
    for _ in range(100):
        midpoint = (left + right) / 2.0
        value = difference(midpoint)
        if right - left <= 1.0e-10:
            break
        if value > 0.0:
            left = midpoint
        else:
            right = midpoint
    return (left + right) / 2.0


def component_values(atom_a: str, atom_b: str, distance: float, x: float):
    value_a = rho(atom_a, x)
    value_b = rho(atom_b, distance - x)
    return value_a, value_b, value_a + value_b


def boundary_estimate(atom_a: str, atom_b: str, distance: float) -> StudyEstimate:
    radius_a = cutoff_radius(atom_a)
    radius_b = cutoff_radius(atom_b)
    separation = distance - radius_a - radius_b

    if atom_a == atom_b:
        x = distance / 2.0
        value_a, value_b, total = component_values(atom_a, atom_b, distance, x)
        status = "low_density_gap" if separation > 0.0 else "ok"
        return StudyEstimate(
            "homonuclear_midpoint", status, x, value_a, value_b, total, separation
        )

    if separation > 0.0:
        x = (radius_a + distance - radius_b) / 2.0
        value_a, value_b, total = component_values(atom_a, atom_b, distance, x)
        return StudyEstimate(
            "cutoff_gap_midpoint",
            "low_density_gap",
            x,
            value_a,
            value_b,
            total,
            separation,
        )

    x = equal_contribution_position(atom_a, atom_b, distance)
    if x is None or not (0.0 < x < distance):
        return StudyEstimate(
            "none", "one_atom_dominates", None, None, None, None, separation
        )
    value_a, value_b, total = component_values(atom_a, atom_b, distance, x)
    return StudyEstimate(
        "equal_proatom_density", "ok", x, value_a, value_b, total, separation
    )


def overlap_interval(atom_a: str, atom_b: str, distance: float) -> tuple[float, float]:
    return (
        max(0.0, distance - cutoff_radius(atom_b)),
        min(distance, cutoff_radius(atom_a)),
    )


def golden_minimize(function, left: float, right: float):
    ratio = (math.sqrt(5.0) - 1.0) / 2.0
    x1 = right - ratio * (right - left)
    x2 = left + ratio * (right - left)
    f1 = function(x1)
    f2 = function(x2)
    for _ in range(70):
        if f1 <= f2:
            right, x2, f2 = x2, x1, f1
            x1 = right - ratio * (right - left)
            f1 = function(x1)
        else:
            left, x1, f1 = x1, x2, f2
            x2 = left + ratio * (right - left)
            f2 = function(x2)
    return (x1, f1) if f1 <= f2 else (x2, f2)


def grid_valleys(atom_a: str, atom_b: str, distance: float, step: float):
    left, right = overlap_interval(atom_a, atom_b, distance)
    if right <= left:
        return [], 0

    count = max(3, math.ceil((right - left) / step) + 1)
    coordinates = [left + (right - left) * i / (count - 1) for i in range(count)]
    equal = equal_contribution_position(atom_a, atom_b, distance)
    if equal is not None and left < equal < right:
        coordinates.append(equal)
    coordinates = sorted(set(coordinates))

    values = [objective(atom_a, atom_b, x, distance) for x in coordinates]
    function = lambda x: objective(atom_a, atom_b, x, distance)
    candidates = []
    for index in range(1, len(coordinates) - 1):
        if values[index] <= values[index - 1] and values[index] <= values[index + 1]:
            candidates.append(
                golden_minimize(
                    function, coordinates[index - 1], coordinates[index + 1]
                )
            )
    return candidates, len(coordinates)


def coalesce(candidates, resolution: float = PRACTICAL_RESOLUTION):
    if not candidates:
        return []
    ordered = sorted(candidates, key=lambda item: item[0])
    groups: list[list[tuple[float, float]]] = []
    for candidate in ordered:
        if not groups or candidate[0] - groups[-1][-1][0] > resolution:
            groups.append([candidate])
        else:
            groups[-1].append(candidate)
    return [min(group, key=lambda item: item[1]) for group in groups]


def selections_compatible(first, second) -> bool:
    if first is None or second is None:
        return first is second
    position_ok = abs(first[0] - second[0]) <= PRACTICAL_RESOLUTION
    density_ok = abs(first[1] - second[1]) <= (
        COMPETITIVE_RELATIVE_GAP * min(first[1], second[1])
    )
    return position_ok and density_ok


def practical_minimum_estimate(atom_a: str, atom_b: str, distance: float) -> StudyEstimate:
    separation = distance - cutoff_radius(atom_a) - cutoff_radius(atom_b)
    if atom_a == atom_b:
        x = distance / 2.0
        value_a, value_b, total = component_values(atom_a, atom_b, distance, x)
        status = "low_density_gap" if separation > 0.0 else "ok"
        return StudyEstimate(
            "homonuclear_midpoint",
            status,
            x,
            value_a,
            value_b,
            total,
            separation,
            search_resolution=PRACTICAL_RESOLUTION,
            search_converged=True,
        )

    left, right = overlap_interval(atom_a, atom_b, distance)
    if right <= left:
        return StudyEstimate(
            "none", "low_density_gap", None, None, None, None, separation
        )

    all_candidates = []
    pass_best = []
    probes = 0
    for step in (INITIAL_STEP, CONFIRM_STEP):
        candidates, count = grid_valleys(atom_a, atom_b, distance, step)
        probes += count
        all_candidates.extend(candidates)
        pass_best.append(min(candidates, key=lambda item: item[1]) if candidates else None)

    need_fallback = not selections_compatible(pass_best[0], pass_best[1])
    if need_fallback:
        candidates, count = grid_valleys(atom_a, atom_b, distance, FALLBACK_STEP)
        probes += count
        all_candidates.extend(candidates)
        pass_best.append(min(candidates, key=lambda item: item[1]) if candidates else None)

    resolved = coalesce(all_candidates)
    if not resolved:
        return StudyEstimate(
            "none",
            "no_resolved_interior_minimum",
            None,
            None,
            None,
            None,
            separation,
            search_resolution=FALLBACK_STEP if need_fallback else CONFIRM_STEP,
            search_converged=not need_fallback,
            probes=probes,
        )

    ordered = sorted(resolved, key=lambda item: item[1])
    selected = ordered[0]
    alternative = ordered[1] if len(ordered) > 1 else None
    converged = (
        selections_compatible(pass_best[-2], pass_best[-1])
        if need_fallback
        else True
    )

    status = "ok"
    relative_gap = None
    if alternative is not None:
        relative_gap = (alternative[1] - selected[1]) / selected[1]
        if relative_gap <= COMPETITIVE_RELATIVE_GAP:
            status = "ambiguous_competing_minima"
    if not converged:
        status = "search_unstable"

    boundary_value = min(
        objective(atom_a, atom_b, 0.0, distance),
        objective(atom_a, atom_b, distance, distance),
    )
    if boundary_value < selected[1] * (1.0 - 2.0e-14):
        status = "boundary_dominated"

    value_a, value_b, total = component_values(
        atom_a, atom_b, distance, selected[0]
    )
    return StudyEstimate(
        "promolecular_density_minimum",
        status,
        selected[0],
        value_a,
        value_b,
        total,
        separation,
        search_resolution=FALLBACK_STEP if need_fallback else CONFIRM_STEP,
        search_converged=converged,
        alternative_position=None if alternative is None else alternative[0],
        relative_depth_gap=relative_gap,
        probes=probes,
    )


def dense_reference_minimum(atom_a: str, atom_b: str, distance: float, step=0.001):
    if atom_a == atom_b:
        x = distance / 2.0
        return x, objective(atom_a, atom_b, x, distance)
    candidates, _ = grid_valleys(atom_a, atom_b, distance, step)
    if not candidates:
        return None
    return min(candidates, key=lambda item: item[1])


def deterministic_cases(count: int, seed: int = 20260713):
    symbols = [element.symbol for element in ar.iter_elements() if element.z <= 103]
    generator = random.Random(seed)
    cases = []
    for index in range(count):
        atom_a = generator.choice(symbols)
        atom_b = generator.choice(symbols)
        if index % 2:
            distance = generator.uniform(0.2, 12.0)
        else:
            distance = 10.0 ** generator.uniform(math.log10(0.2), math.log10(12.0))
        cases.append((atom_a, atom_b, distance))
    return cases


def run_comparison(count=300):
    cases = deterministic_cases(count)
    practical = []
    reference = []
    start = time.perf_counter()
    for case in cases:
        practical.append(practical_minimum_estimate(*case))
    practical_seconds = time.perf_counter() - start

    start = time.perf_counter()
    for case in cases:
        reference.append(dense_reference_minimum(*case))
    reference_seconds = time.perf_counter() - start

    status_mismatch = 0
    errors = []
    for case, estimate, expected in zip(cases, practical, reference):
        found = estimate.position is not None
        expected_found = expected is not None
        if found != expected_found:
            # A sub-resolution dense-reference minimum can intentionally be absent.
            status_mismatch += 1
            continue
        if found:
            dx = abs(estimate.position - expected[0])
            relative_density = abs(estimate.rho_sum - expected[1]) / expected[1]
            errors.append((dx, relative_density, case, estimate, expected))

    return {
        "cases": cases,
        "practical": practical,
        "reference": reference,
        "practical_seconds": practical_seconds,
        "reference_seconds": reference_seconds,
        "status_mismatch": status_mismatch,
        "errors": errors,
    }

from __future__ import annotations from bisect import bisect_right from dataclasses import dataclass from functools import lru_cache import math import random import time from pathlib import Path import sys ROOT = Path.cwd().resolve() if not (ROOT / "src" / "atomref").is_dir(): for candidate in (ROOT.parent, ROOT.parent.parent): if (candidate / "src" / "atomref").is_dir(): ROOT = candidate break if not (ROOT / "src" / "atomref").is_dir(): raise RuntimeError("Run this notebook from the atomref repository checkout") sys.path.insert(0, str(ROOT / "src")) import atomref as ar try: import matplotlib.pyplot as plt except ImportError: # The package itself has no plotting dependency. plt = None ASSET_DIR = ROOT / "docs" / "assets" / "ias-method-study" ASSET_DIR.mkdir(parents=True, exist_ok=True) TAIL_CUTOFF = 1.0e-4 PRACTICAL_RESOLUTION = 0.01 INITIAL_STEP = 0.02 CONFIRM_STEP = 0.01 FALLBACK_STEP = 0.005 COMPETITIVE_RELATIVE_GAP = 1.0e-4 @dataclass(frozen=True) class StudyEstimate: method: str status: str position: float | None rho_a: float | None rho_b: float | None rho_sum: float | None contour_separation: float search_resolution: float | None = None search_converged: bool | None = None alternative_position: float | None = None relative_depth_gap: float | None = None probes: int = 0 @lru_cache(maxsize=None) def profile(symbol: str): result = ar.get_proatomic_density_profile(symbol) if result is None: raise ValueError(f"missing profile: {symbol}") return result @lru_cache(maxsize=None) def cutoff_radius(symbol: str) -> float: p = profile(symbol) right = next(i for i, value in enumerate(p.densities) if value <= TAIL_CUTOFF) if p.densities[right] == TAIL_CUTOFF: return p.radii[right] left = right - 1 lr0 = math.log(p.radii[left]) lr1 = math.log(p.radii[right]) ld0 = math.log(p.densities[left]) ld1 = math.log(p.densities[right]) fraction = (math.log(TAIL_CUTOFF) - ld0) / (ld1 - ld0) return math.exp(lr0 + fraction * (lr1 - lr0)) def rho(symbol: str, radius: float) -> float: return profile(symbol)._evaluate_bohr(radius) def objective(atom_a: str, atom_b: str, x: float, distance: float) -> float: return rho(atom_a, x) + rho(atom_b, distance - x) def continuous_log_density(symbol: str, radius: float) -> float: p = profile(symbol) if radius <= p.radii[0]: return math.log(p.densities[0]) index = min(bisect_right(p.radii, radius) - 1, len(p.radii) - 2) lr0 = math.log(p.radii[index]) lr1 = math.log(p.radii[index + 1]) ld0 = math.log(p.densities[index]) ld1 = math.log(p.densities[index + 1]) fraction = (math.log(radius) - lr0) / (lr1 - lr0) return ld0 + fraction * (ld1 - ld0) def equal_contribution_position(atom_a: str, atom_b: str, distance: float) -> float | None: if atom_a == atom_b: return distance / 2.0 def difference(x: float) -> float: return continuous_log_density(atom_a, x) - continuous_log_density( atom_b, distance - x ) left_value = difference(0.0) right_value = difference(distance) if left_value == 0.0: return 0.0 if right_value == 0.0: return distance if left_value * right_value > 0.0: return None left = 0.0 right = distance for _ in range(100): midpoint = (left + right) / 2.0 value = difference(midpoint) if right - left <= 1.0e-10: break if value > 0.0: left = midpoint else: right = midpoint return (left + right) / 2.0 def component_values(atom_a: str, atom_b: str, distance: float, x: float): value_a = rho(atom_a, x) value_b = rho(atom_b, distance - x) return value_a, value_b, value_a + value_b def boundary_estimate(atom_a: str, atom_b: str, distance: float) -> StudyEstimate: radius_a = cutoff_radius(atom_a) radius_b = cutoff_radius(atom_b) separation = distance - radius_a - radius_b if atom_a == atom_b: x = distance / 2.0 value_a, value_b, total = component_values(atom_a, atom_b, distance, x) status = "low_density_gap" if separation > 0.0 else "ok" return StudyEstimate( "homonuclear_midpoint", status, x, value_a, value_b, total, separation ) if separation > 0.0: x = (radius_a + distance - radius_b) / 2.0 value_a, value_b, total = component_values(atom_a, atom_b, distance, x) return StudyEstimate( "cutoff_gap_midpoint", "low_density_gap", x, value_a, value_b, total, separation, ) x = equal_contribution_position(atom_a, atom_b, distance) if x is None or not (0.0 < x < distance): return StudyEstimate( "none", "one_atom_dominates", None, None, None, None, separation ) value_a, value_b, total = component_values(atom_a, atom_b, distance, x) return StudyEstimate( "equal_proatom_density", "ok", x, value_a, value_b, total, separation ) def overlap_interval(atom_a: str, atom_b: str, distance: float) -> tuple[float, float]: return ( max(0.0, distance - cutoff_radius(atom_b)), min(distance, cutoff_radius(atom_a)), ) def golden_minimize(function, left: float, right: float): ratio = (math.sqrt(5.0) - 1.0) / 2.0 x1 = right - ratio * (right - left) x2 = left + ratio * (right - left) f1 = function(x1) f2 = function(x2) for _ in range(70): if f1 <= f2: right, x2, f2 = x2, x1, f1 x1 = right - ratio * (right - left) f1 = function(x1) else: left, x1, f1 = x1, x2, f2 x2 = left + ratio * (right - left) f2 = function(x2) return (x1, f1) if f1 <= f2 else (x2, f2) def grid_valleys(atom_a: str, atom_b: str, distance: float, step: float): left, right = overlap_interval(atom_a, atom_b, distance) if right <= left: return [], 0 count = max(3, math.ceil((right - left) / step) + 1) coordinates = [left + (right - left) * i / (count - 1) for i in range(count)] equal = equal_contribution_position(atom_a, atom_b, distance) if equal is not None and left < equal < right: coordinates.append(equal) coordinates = sorted(set(coordinates)) values = [objective(atom_a, atom_b, x, distance) for x in coordinates] function = lambda x: objective(atom_a, atom_b, x, distance) candidates = [] for index in range(1, len(coordinates) - 1): if values[index] <= values[index - 1] and values[index] <= values[index + 1]: candidates.append( golden_minimize( function, coordinates[index - 1], coordinates[index + 1] ) ) return candidates, len(coordinates) def coalesce(candidates, resolution: float = PRACTICAL_RESOLUTION): if not candidates: return [] ordered = sorted(candidates, key=lambda item: item[0]) groups: list[list[tuple[float, float]]] = [] for candidate in ordered: if not groups or candidate[0] - groups[-1][-1][0] > resolution: groups.append([candidate]) else: groups[-1].append(candidate) return [min(group, key=lambda item: item[1]) for group in groups] def selections_compatible(first, second) -> bool: if first is None or second is None: return first is second position_ok = abs(first[0] - second[0]) <= PRACTICAL_RESOLUTION density_ok = abs(first[1] - second[1]) <= ( COMPETITIVE_RELATIVE_GAP * min(first[1], second[1]) ) return position_ok and density_ok def practical_minimum_estimate(atom_a: str, atom_b: str, distance: float) -> StudyEstimate: separation = distance - cutoff_radius(atom_a) - cutoff_radius(atom_b) if atom_a == atom_b: x = distance / 2.0 value_a, value_b, total = component_values(atom_a, atom_b, distance, x) status = "low_density_gap" if separation > 0.0 else "ok" return StudyEstimate( "homonuclear_midpoint", status, x, value_a, value_b, total, separation, search_resolution=PRACTICAL_RESOLUTION, search_converged=True, ) left, right = overlap_interval(atom_a, atom_b, distance) if right <= left: return StudyEstimate( "none", "low_density_gap", None, None, None, None, separation ) all_candidates = [] pass_best = [] probes = 0 for step in (INITIAL_STEP, CONFIRM_STEP): candidates, count = grid_valleys(atom_a, atom_b, distance, step) probes += count all_candidates.extend(candidates) pass_best.append(min(candidates, key=lambda item: item[1]) if candidates else None) need_fallback = not selections_compatible(pass_best[0], pass_best[1]) if need_fallback: candidates, count = grid_valleys(atom_a, atom_b, distance, FALLBACK_STEP) probes += count all_candidates.extend(candidates) pass_best.append(min(candidates, key=lambda item: item[1]) if candidates else None) resolved = coalesce(all_candidates) if not resolved: return StudyEstimate( "none", "no_resolved_interior_minimum", None, None, None, None, separation, search_resolution=FALLBACK_STEP if need_fallback else CONFIRM_STEP, search_converged=not need_fallback, probes=probes, ) ordered = sorted(resolved, key=lambda item: item[1]) selected = ordered[0] alternative = ordered[1] if len(ordered) > 1 else None converged = ( selections_compatible(pass_best[-2], pass_best[-1]) if need_fallback else True ) status = "ok" relative_gap = None if alternative is not None: relative_gap = (alternative[1] - selected[1]) / selected[1] if relative_gap <= COMPETITIVE_RELATIVE_GAP: status = "ambiguous_competing_minima" if not converged: status = "search_unstable" boundary_value = min( objective(atom_a, atom_b, 0.0, distance), objective(atom_a, atom_b, distance, distance), ) if boundary_value < selected[1] * (1.0 - 2.0e-14): status = "boundary_dominated" value_a, value_b, total = component_values( atom_a, atom_b, distance, selected[0] ) return StudyEstimate( "promolecular_density_minimum", status, selected[0], value_a, value_b, total, separation, search_resolution=FALLBACK_STEP if need_fallback else CONFIRM_STEP, search_converged=converged, alternative_position=None if alternative is None else alternative[0], relative_depth_gap=relative_gap, probes=probes, ) def dense_reference_minimum(atom_a: str, atom_b: str, distance: float, step=0.001): if atom_a == atom_b: x = distance / 2.0 return x, objective(atom_a, atom_b, x, distance) candidates, _ = grid_valleys(atom_a, atom_b, distance, step) if not candidates: return None return min(candidates, key=lambda item: item[1]) def deterministic_cases(count: int, seed: int = 20260713): symbols = [element.symbol for element in ar.iter_elements() if element.z <= 103] generator = random.Random(seed) cases = [] for index in range(count): atom_a = generator.choice(symbols) atom_b = generator.choice(symbols) if index % 2: distance = generator.uniform(0.2, 12.0) else: distance = 10.0 ** generator.uniform(math.log10(0.2), math.log10(12.0)) cases.append((atom_a, atom_b, distance)) return cases def run_comparison(count=300): cases = deterministic_cases(count) practical = [] reference = [] start = time.perf_counter() for case in cases: practical.append(practical_minimum_estimate(*case)) practical_seconds = time.perf_counter() - start start = time.perf_counter() for case in cases: reference.append(dense_reference_minimum(*case)) reference_seconds = time.perf_counter() - start status_mismatch = 0 errors = [] for case, estimate, expected in zip(cases, practical, reference): found = estimate.position is not None expected_found = expected is not None if found != expected_found: # A sub-resolution dense-reference minimum can intentionally be absent. status_mismatch += 1 continue if found: dx = abs(estimate.position - expected[0]) relative_density = abs(estimate.rho_sum - expected[1]) / expected[1] errors.append((dx, relative_density, case, estimate, expected)) return { "cases": cases, "practical": practical, "reference": reference, "practical_seconds": practical_seconds, "reference_seconds": reference_seconds, "status_mismatch": status_mismatch, "errors": errors, }

Cutoff radii across H–Lr¶

Every packaged profile is monotone over its stored positive segments, so the cutoff radius is unique and can be obtained by directly inverting one log–log segment. The plot shows that the chosen contour is deliberately diffuse: the radii span roughly 1.8–4.2 Å.

In [2]:

cutoff_rows = []
for element in ar.iter_elements():
    if element.z <= 103:
        radius_bohr = cutoff_radius(element.symbol)
        cutoff_rows.append(
            (element.z, element.symbol, radius_bohr, radius_bohr * ar.proatoms.BOHR_TO_ANGSTROM)
        )

smallest = min(cutoff_rows, key=lambda row: row[2])
largest = max(cutoff_rows, key=lambda row: row[2])
print(
    f"Smallest cutoff radius: {smallest[1]} = {smallest[2]:.6f} bohr "
    f"({smallest[3]:.6f} Å)"
)
print(
    f"Largest cutoff radius:  {largest[1]} = {largest[2]:.6f} bohr "
    f"({largest[3]:.6f} Å)"
)
print("Representative radii:")
for symbol in ("H", "He", "Li", "C", "O", "Fe", "U", "Lr"):
    row = next(item for item in cutoff_rows if item[1] == symbol)
    print(f"  {symbol:>2}: {row[2]:.6f} bohr = {row[3]:.6f} Å")

figure = None
if plt is not None:
    figure, axis = plt.subplots()
    axis.plot([row[0] for row in cutoff_rows], [row[3] for row in cutoff_rows])
    axis.set_xlabel("Atomic number")
    axis.set_ylabel("Cutoff radius (Å)")
    axis.set_title(r"Neutral-proatom radius at $\rho=10^{-4}$ electron/bohr$^3$")
    figure.tight_layout()
    figure.savefig(ASSET_DIR / "cutoff-radii.png", dpi=160)
    plt.close(figure)
else:
    print("Matplotlib is unavailable; the cutoff-radius plot was skipped.")

figure

cutoff_rows = [] for element in ar.iter_elements(): if element.z <= 103: radius_bohr = cutoff_radius(element.symbol) cutoff_rows.append( (element.z, element.symbol, radius_bohr, radius_bohr * ar.proatoms.BOHR_TO_ANGSTROM) ) smallest = min(cutoff_rows, key=lambda row: row[2]) largest = max(cutoff_rows, key=lambda row: row[2]) print( f"Smallest cutoff radius: {smallest[1]} = {smallest[2]:.6f} bohr " f"({smallest[3]:.6f} Å)" ) print( f"Largest cutoff radius: {largest[1]} = {largest[2]:.6f} bohr " f"({largest[3]:.6f} Å)" ) print("Representative radii:") for symbol in ("H", "He", "Li", "C", "O", "Fe", "U", "Lr"): row = next(item for item in cutoff_rows if item[1] == symbol) print(f" {symbol:>2}: {row[2]:.6f} bohr = {row[3]:.6f} Å") figure = None if plt is not None: figure, axis = plt.subplots() axis.plot([row[0] for row in cutoff_rows], [row[3] for row in cutoff_rows]) axis.set_xlabel("Atomic number") axis.set_ylabel("Cutoff radius (Å)") axis.set_title(r"Neutral-proatom radius at $\rho=10^{-4}$ electron/bohr$^3$") figure.tight_layout() figure.savefig(ASSET_DIR / "cutoff-radii.png", dpi=160) plt.close(figure) else: print("Matplotlib is unavailable; the cutoff-radius plot was skipped.") figure

Smallest cutoff radius: He = 3.373985 bohr (1.785436 Å)
Largest cutoff radius:  Ba = 7.966363 bohr (4.215618 Å)
Representative radii:
   H: 4.128012 bohr = 2.184450 Å
  He: 3.373985 bohr = 1.785436 Å
  Li: 6.469571 bohr = 3.423549 Å
   C: 4.895923 bohr = 2.590811 Å
   O: 4.333035 bohr = 2.292943 Å
  Fe: 6.132435 bohr = 3.245145 Å
   U: 7.398370 bohr = 3.915049 Å
  Lr: 7.736203 bohr = 4.093822 Å

Out[2]:

Representative pairwise results¶

The same profiles can support two useful but different answers. The boundary mode is the recommended geometry-facing default. The minimum mode is retained as a Bader-oriented neutral-promolecular proxy and may explicitly decline to return a resolved minimum.

In [3]:

examples = [
    ("Li", "Li", 5.0),
    ("H", "U", 3.8),
    ("C", "O", 1.5),
    ("H", "O", 12.0),
    ("H", "U", 1.58),
]

header = (
    "pair", "R/bohr", "boundary method", "boundary x", "minimum status", "minimum x"
)
print(f"{header[0]:<8} {header[1]:>9}  {header[2]:<25} {header[3]:>12}  {header[4]:<31} {header[5]:>12}")
for atom_a, atom_b, distance in examples:
    boundary = boundary_estimate(atom_a, atom_b, distance)
    minimum = practical_minimum_estimate(atom_a, atom_b, distance)
    bx = "None" if boundary.position is None else f"{boundary.position:.9f}"
    mx = "None" if minimum.position is None else f"{minimum.position:.9f}"
    print(
        f"{atom_a+'-'+atom_b:<8} {distance:9.4f}  {boundary.method:<25} {bx:>12}  "
        f"{minimum.status:<31} {mx:>12}"
    )

examples = [ ("Li", "Li", 5.0), ("H", "U", 3.8), ("C", "O", 1.5), ("H", "O", 12.0), ("H", "U", 1.58), ] header = ( "pair", "R/bohr", "boundary method", "boundary x", "minimum status", "minimum x" ) print(f"{header[0]:<8} {header[1]:>9} {header[2]:<25} {header[3]:>12} {header[4]:<31} {header[5]:>12}") for atom_a, atom_b, distance in examples: boundary = boundary_estimate(atom_a, atom_b, distance) minimum = practical_minimum_estimate(atom_a, atom_b, distance) bx = "None" if boundary.position is None else f"{boundary.position:.9f}" mx = "None" if minimum.position is None else f"{minimum.position:.9f}" print( f"{atom_a+'-'+atom_b:<8} {distance:9.4f} {boundary.method:<25} {bx:>12} " f"{minimum.status:<31} {mx:>12}" )

pair        R/bohr  boundary method             boundary x  minimum status                     minimum x
Li-Li       5.0000  homonuclear_midpoint       2.500000000  ok                               2.500000000
H-U         3.8000  equal_proatom_density      1.184253035  ok                               1.191871522
C-O         1.5000  equal_proatom_density      0.556215252  ok                               0.611760239
H-O        12.0000  cutoff_gap_midpoint        5.897488625  low_density_gap                         None
H-U         1.5800  equal_proatom_density      0.017084468  no_resolved_interior_minimum            None

Why homonuclear symmetry overrides raw shell minima¶

For Li–Li at 5 bohr, a dense local-minimum scan finds symmetric off-centre minima slightly below the midpoint. Returning either off-centre point would be a poor pairwise separator for identical atoms. Both estimator modes therefore return $R/2$ and preserve the raw profile behaviour only as a diagnostic.

In [4]:

atom_a, atom_b, distance = "Li", "Li", 5.0
xs = [1.75 + (3.25 - 1.75) * i / 1600 for i in range(1601)]
ys = [objective(atom_a, atom_b, x, distance) for x in xs]
raw_valleys, _ = grid_valleys(atom_a, atom_b, distance, 0.001)

print("Raw fine-grid valley coordinates near the internuclear region:")
for x, value in raw_valleys:
    if 1.75 <= x <= 3.25:
        print(f"  x={x:.12f} bohr, rho_sum={value:.15g}")
print(f"Symmetry-selected coordinate: {distance/2:.12f} bohr")

figure = None
if plt is not None:
    figure, axis = plt.subplots()
    axis.plot(xs, ys, label=r"$\rho_{Li}(x)+\rho_{Li}(R-x)$")
    shown = [(x, value) for x, value in raw_valleys if 1.75 <= x <= 3.25]
    axis.scatter([item[0] for item in shown], [item[1] for item in shown], label="raw local valleys")
    axis.axvline(distance / 2.0, linestyle="--", label="symmetry midpoint")
    axis.set_xlabel("Position from atom A (bohr)")
    axis.set_ylabel(r"Promolecular line density (electron/bohr$^3$)")
    axis.set_title("Li–Li at R = 5 bohr")
    axis.legend()
    figure.tight_layout()
    figure.savefig(ASSET_DIR / "li-li-symmetry.png", dpi=160)
    plt.close(figure)
else:
    print("Matplotlib is unavailable; the Li–Li plot was skipped.")

figure

atom_a, atom_b, distance = "Li", "Li", 5.0 xs = [1.75 + (3.25 - 1.75) * i / 1600 for i in range(1601)] ys = [objective(atom_a, atom_b, x, distance) for x in xs] raw_valleys, _ = grid_valleys(atom_a, atom_b, distance, 0.001) print("Raw fine-grid valley coordinates near the internuclear region:") for x, value in raw_valleys: if 1.75 <= x <= 3.25: print(f" x={x:.12f} bohr, rho_sum={value:.15g}") print(f"Symmetry-selected coordinate: {distance/2:.12f} bohr") figure = None if plt is not None: figure, axis = plt.subplots() axis.plot(xs, ys, label=r"$\rho_{Li}(x)+\rho_{Li}(R-x)$") shown = [(x, value) for x, value in raw_valleys if 1.75 <= x <= 3.25] axis.scatter([item[0] for item in shown], [item[1] for item in shown], label="raw local valleys") axis.axvline(distance / 2.0, linestyle="--", label="symmetry midpoint") axis.set_xlabel("Position from atom A (bohr)") axis.set_ylabel(r"Promolecular line density (electron/bohr$^3$)") axis.set_title("Li–Li at R = 5 bohr") axis.legend() figure.tight_layout() figure.savefig(ASSET_DIR / "li-li-symmetry.png", dpi=160) plt.close(figure) else: print("Matplotlib is unavailable; the Li–Li plot was skipped.") figure

Raw fine-grid valley coordinates near the internuclear region:
  x=2.118430478685 bohr, rho_sum=0.00654500903481724
  x=2.472517881659 bohr, rho_sum=0.00658890668816022
  x=2.500000037765 bohr, rho_sum=0.00658920833374835
  x=2.527482056015 bohr, rho_sum=0.00658890668816022
  x=2.881569505107 bohr, rho_sum=0.00654500903481724
Symmetry-selected coordinate: 2.500000000000 bohr

Out[4]:

Boundary and minimum are different scientific definitions¶

For C–O at 1.5 bohr, the equal-contribution point and the selected promolecular minimum differ visibly. This is not a solver failure: the first balances the two reference-atom contributions, while the second minimizes their sum.

In [5]:

atom_a, atom_b, distance = "C", "O", 1.5
boundary = boundary_estimate(atom_a, atom_b, distance)
minimum = practical_minimum_estimate(atom_a, atom_b, distance)
xs = [distance * i / 1000 for i in range(1001)]
rho_a_values = [rho(atom_a, x) for x in xs]
rho_b_values = [rho(atom_b, distance - x) for x in xs]
total_values = [a + b for a, b in zip(rho_a_values, rho_b_values)]

print(f"Boundary coordinate: {boundary.position:.12f} bohr")
print(f"Minimum coordinate:  {minimum.position:.12f} bohr")
print(f"Coordinate difference: {abs(boundary.position-minimum.position):.12f} bohr")
print(f"Boundary component densities: {boundary.rho_a:.12g}, {boundary.rho_b:.12g}")
print(f"Minimum component densities:  {minimum.rho_a:.12g}, {minimum.rho_b:.12g}")

figure = None
if plt is not None:
    figure, axis = plt.subplots()
    axis.plot(xs, rho_a_values, label=r"$\rho_C(x)$")
    axis.plot(xs, rho_b_values, label=r"$\rho_O(R-x)$")
    axis.plot(xs, total_values, label="sum")
    axis.axvline(boundary.position, linestyle="--", label="proatomic boundary")
    axis.axvline(minimum.position, linestyle=":", label="promolecular minimum")
    axis.set_xlim(0.35, 0.90)
    axis.set_ylim(0.15, 2.8)
    axis.set_xlabel("Position from C (bohr)")
    axis.set_ylabel(r"Density (electron/bohr$^3$)")
    axis.set_title("C–O interatomic region at R = 1.5 bohr")
    axis.legend(loc="upper right")
    figure.tight_layout()
    figure.savefig(ASSET_DIR / "c-o-method-comparison.png", dpi=160)
    plt.close(figure)
else:
    print("Matplotlib is unavailable; the C–O plot was skipped.")

figure

atom_a, atom_b, distance = "C", "O", 1.5 boundary = boundary_estimate(atom_a, atom_b, distance) minimum = practical_minimum_estimate(atom_a, atom_b, distance) xs = [distance * i / 1000 for i in range(1001)] rho_a_values = [rho(atom_a, x) for x in xs] rho_b_values = [rho(atom_b, distance - x) for x in xs] total_values = [a + b for a, b in zip(rho_a_values, rho_b_values)] print(f"Boundary coordinate: {boundary.position:.12f} bohr") print(f"Minimum coordinate: {minimum.position:.12f} bohr") print(f"Coordinate difference: {abs(boundary.position-minimum.position):.12f} bohr") print(f"Boundary component densities: {boundary.rho_a:.12g}, {boundary.rho_b:.12g}") print(f"Minimum component densities: {minimum.rho_a:.12g}, {minimum.rho_b:.12g}") figure = None if plt is not None: figure, axis = plt.subplots() axis.plot(xs, rho_a_values, label=r"$\rho_C(x)$") axis.plot(xs, rho_b_values, label=r"$\rho_O(R-x)$") axis.plot(xs, total_values, label="sum") axis.axvline(boundary.position, linestyle="--", label="proatomic boundary") axis.axvline(minimum.position, linestyle=":", label="promolecular minimum") axis.set_xlim(0.35, 0.90) axis.set_ylim(0.15, 2.8) axis.set_xlabel("Position from C (bohr)") axis.set_ylabel(r"Density (electron/bohr$^3$)") axis.set_title("C–O interatomic region at R = 1.5 bohr") axis.legend(loc="upper right") figure.tight_layout() figure.savefig(ASSET_DIR / "c-o-method-comparison.png", dpi=160) plt.close(figure) else: print("Matplotlib is unavailable; the C–O plot was skipped.") figure

Boundary coordinate: 0.556215252014 bohr
Minimum coordinate:  0.611760239435 bohr
Coordinate difference: 0.055544987420 bohr
Boundary component densities: 0.414069174671, 0.414069174631
Minimum component densities:  0.325849689041, 0.475793723795

Out[5]:

Practical minimum search against a slower reference¶

The public minimum mode is not intended to enumerate every mathematical local minimum. It asks whether a valley is resolved at 0.01 bohr and checks it at two grid scales, with one finer fallback. The comparison below uses 300 deterministic H–Lr cases, plus seven targeted numerical edge cases, against a slower 0.001-bohr fine-grid search over the same cutoff-bounded interval.

A dense-reference minimum may be intentionally rejected when it is narrower than the documented resolution or lies immediately beside a nucleus at an extreme short distance.

In [6]:

adversarial_cases = [
    ("Ni", "Te", 1.5897045597171517),
    ("Al", "Pb", 1.3645306063699802),
    ("Pr", "Re", 0.7592738817632778),
    ("C", "Lu", 0.2217301367970158),
    ("He", "Ru", 0.6552909074349161),
    ("H", "U", 1.58),
    ("C", "O", 1.5),
]
cases = deterministic_cases(300) + adversarial_cases

start = time.perf_counter()
practical_results = [practical_minimum_estimate(*case) for case in cases]
practical_seconds = time.perf_counter() - start

start = time.perf_counter()
reference_results = [dense_reference_minimum(*case) for case in cases]
reference_seconds = time.perf_counter() - start

joint_errors = []
intentional_nonmatches = []
both_absent = 0
for case, practical, reference in zip(cases, practical_results, reference_results):
    if practical.position is None and reference is None:
        both_absent += 1
    elif (practical.position is None) != (reference is None):
        intentional_nonmatches.append((case, practical, reference))
    else:
        joint_errors.append(
            (
                abs(practical.position - reference[0]),
                abs(practical.rho_sum - reference[1]) / reference[1],
                case,
            )
        )

max_position = max(joint_errors, key=lambda item: item[0])
max_density = max(joint_errors, key=lambda item: item[1])
within_resolution = sum(error[0] <= PRACTICAL_RESOLUTION for error in joint_errors)

print(f"Cases: {len(cases)}")
print(f"Both methods returned a minimum: {len(joint_errors)}")
print(f"Both returned no resolved minimum: {both_absent}")
print(f"Practical-resolution nonmatches: {len(intentional_nonmatches)}")
print(f"Joint results within 0.01 bohr: {within_resolution}/{len(joint_errors)}")
print(
    f"Largest position difference: {max_position[0]:.6g} bohr "
    f"for {max_position[2][0]}-{max_position[2][1]} at R={max_position[2][2]:.6g} bohr"
)
print(f"Largest relative density difference: {max_density[1]:.6g}")
print(f"Practical search time: {practical_seconds:.3f} s")
print(f"Slower reference time: {reference_seconds:.3f} s")

if intentional_nonmatches:
    print("\nIntentionally unresolved dense-reference minima:")
    for case, practical, reference in intentional_nonmatches:
        print(
            f"  {case[0]}-{case[1]} R={case[2]:.12g}: "
            f"status={practical.status}, reference_x={reference[0]:.12g} bohr"
        )

adversarial_cases = [ ("Ni", "Te", 1.5897045597171517), ("Al", "Pb", 1.3645306063699802), ("Pr", "Re", 0.7592738817632778), ("C", "Lu", 0.2217301367970158), ("He", "Ru", 0.6552909074349161), ("H", "U", 1.58), ("C", "O", 1.5), ] cases = deterministic_cases(300) + adversarial_cases start = time.perf_counter() practical_results = [practical_minimum_estimate(*case) for case in cases] practical_seconds = time.perf_counter() - start start = time.perf_counter() reference_results = [dense_reference_minimum(*case) for case in cases] reference_seconds = time.perf_counter() - start joint_errors = [] intentional_nonmatches = [] both_absent = 0 for case, practical, reference in zip(cases, practical_results, reference_results): if practical.position is None and reference is None: both_absent += 1 elif (practical.position is None) != (reference is None): intentional_nonmatches.append((case, practical, reference)) else: joint_errors.append( ( abs(practical.position - reference[0]), abs(practical.rho_sum - reference[1]) / reference[1], case, ) ) max_position = max(joint_errors, key=lambda item: item[0]) max_density = max(joint_errors, key=lambda item: item[1]) within_resolution = sum(error[0] <= PRACTICAL_RESOLUTION for error in joint_errors) print(f"Cases: {len(cases)}") print(f"Both methods returned a minimum: {len(joint_errors)}") print(f"Both returned no resolved minimum: {both_absent}") print(f"Practical-resolution nonmatches: {len(intentional_nonmatches)}") print(f"Joint results within 0.01 bohr: {within_resolution}/{len(joint_errors)}") print( f"Largest position difference: {max_position[0]:.6g} bohr " f"for {max_position[2][0]}-{max_position[2][1]} at R={max_position[2][2]:.6g} bohr" ) print(f"Largest relative density difference: {max_density[1]:.6g}") print(f"Practical search time: {practical_seconds:.3f} s") print(f"Slower reference time: {reference_seconds:.3f} s") if intentional_nonmatches: print("\nIntentionally unresolved dense-reference minima:") for case, practical, reference in intentional_nonmatches: print( f" {case[0]}-{case[1]} R={case[2]:.12g}: " f"status={practical.status}, reference_x={reference[0]:.12g} bohr" )

Cases: 307
Both methods returned a minimum: 298
Both returned no resolved minimum: 7
Practical-resolution nonmatches: 2
Joint results within 0.01 bohr: 298/298
Largest position difference: 0.0012907 bohr for C-O at R=1.5 bohr
Largest relative density difference: 5.96583e-06
Practical search time: 0.244 s
Slower reference time: 1.076 s

Intentionally unresolved dense-reference minima:
  C-Lu R=0.221730136797: status=no_resolved_interior_minimum, reference_x=0.00284120268288 bohr
  He-Ru R=0.655290907435: status=no_resolved_interior_minimum, reference_x=0.0112558333887 bohr

Cached timing comparison¶

Timing is machine-dependent, so these figures are evidence rather than portable thresholds. The important point is the scale difference: the boundary mode needs one cached contour calculation plus a monotone bisection, while the optional minimum mode remains below about a millisecond per ordinary pair in this standalone pure-Python study.

In [7]:

import statistics

timing_cases = deterministic_cases(1000)
for case in timing_cases[:30]:
    boundary_estimate(*case)
    practical_minimum_estimate(*case)

for label, function in (
    ("boundary", boundary_estimate),
    ("minimum", practical_minimum_estimate),
):
    samples = []
    for _ in range(5):
        start = time.perf_counter()
        [function(*case) for case in timing_cases]
        samples.append((time.perf_counter() - start) * 1000.0 / len(timing_cases))
    print(f"{label:>8}: median {statistics.median(samples):.4f} ms/pair; runs={samples}")

import statistics timing_cases = deterministic_cases(1000) for case in timing_cases[:30]: boundary_estimate(*case) practical_minimum_estimate(*case) for label, function in ( ("boundary", boundary_estimate), ("minimum", practical_minimum_estimate), ): samples = [] for _ in range(5): start = time.perf_counter() [function(*case) for case in timing_cases] samples.append((time.perf_counter() - start) * 1000.0 / len(timing_cases)) print(f"{label:>8}: median {statistics.median(samples):.4f} ms/pair; runs={samples}")

boundary: median 0.0541 ms/pair; runs=[0.05746927700010929, 0.055642846000182544, 0.05410611999991488, 0.05177689599986479, 0.050906101000009585]

 minimum: median 0.7887 ms/pair; runs=[0.8100914009999087, 0.7820286119999764, 0.7892370699998992, 0.7887017269999888, 0.7813305360000413]

Why exhaustive local-minimum enumeration is out of scope¶

The exhaustive analysis showed that an exact “return every mathematical local minimum” contract really does require disproportionate numerical machinery:

Case	Exact-topology observation	Consequence for the numerical policy
H–F near a knot	A minimum occurred 27 binary64 coordinates below a knot and needed high-precision derivative signs.	Important only when every microminimum is binding.
C–O exact knot	The stored knot and its immediate neighbor could differ by one density ULP.	Exact-knot topology is below the declared 0.01-bohr resolution.
Li–Li, 5 bohr	Five local minima were found; the two off-centre minima were slightly deepest.	A raw minimum can violate the desired homonuclear separator symmetry.
H–U, 1.58 bohr	Fourteen isolated minima were found while a boundary value was lower.	Counting minima does not produce a useful IAS estimate.
At–Bk event collision	Two exact-distinct profile events rounded to one float and caused orientation-dependent candidate counts in an event-based enumerator.	Exact event identity is unnecessary for the documented resolution-limited estimator.
Ne–Ne open segment	Ordinary log/exp value error exceeded a stringent global tie tolerance.	Numerical tie lists require a separate high-precision value model.

These cases illustrate genuine shell-structure and floating-point effects. They lie outside the documented resolution-limited semantics: boundary mode is unaffected, and minimum mode deliberately coalesces or rejects sub-resolution structure.

Public interface and scope¶

The public IAS API has three entry points:

estimate_proatomic_boundary(...)
estimate_promolecular_density_minimum(...)
estimate_ias_position(..., mode="boundary" | "minimum")

boundary is the default. It is the best fit for stable geometry and Laguerre/Voronoi calibration. minimum remains available for Bader-oriented research, but it searches only the significant-overlap region, returns one resolved candidate rather than every local minimum, and reports ambiguity or a non-result explicitly.

For identical atoms, both modes return exactly R/2. Extreme short-distance cases are handled through explicit statuses such as one_atom_dominates or no_resolved_interior_minimum, rather than a universal distance exclusion.

The remaining scientific limitation is deliberate: neither neutral-proatom mode has yet been calibrated against a curated set of molecular-density QTAIM interatomic-surface coordinates.