Comparing periodic stockholder methods

This notebook applies Hirshfeld, Hirshfeld-I, MBIS, LISA, and AVH-B to the same periodic density and uniform grid. The example illustrates method-specific basis requirements and common numerical checks; it is not intended as a physical model of crystalline hydrogen.

import numpy as np

from horton_part import ProAtomDB
from horton_part.periodic import partition_periodic
from periodic_setup import make_demo_spline_basis, make_demo_system

coordinates, numbers, grid, density, lisa_basis = make_demo_system()
spline_basis = make_demo_spline_basis()
print(f"integrated electrons: {grid.integrate(density):.8f}")

integrated electrons: 1.99999972

Basis requirements

MBIS constructs its own Slater basis. LISA uses normalized exponential functions with fixed exponents and optimized nonnegative populations. Hirshfeld, Hirshfeld-I, and AVH use electron-normalized radial-spline densities of isolated charged atoms. AVH-B selects the physically bound states from that library. The neutral spline schema is shared with finite-system calculations; no basis conversion is required.

finite_proatoms = ProAtomDB.from_spline_file(spline_basis)
print("finite-system elements:", finite_proatoms.get_numbers())
print("available H charge states:", finite_proatoms.get_charges(1))

finite-system elements: [1]
available H charge states: [1, 0, -1]

calculations = {
    "Hirshfeld": ("hirshfeld", {"basis": spline_basis}),
    "Hirshfeld-I": ("hirshfeld-i", {"basis": spline_basis}),
    "MBIS": ("mbis", {}),
    "LISA": ("lisa", {"basis": lisa_basis}),
    "AVH-B": ("avh", {"basis": spline_basis, "avh_variant": "B"}),
}
results = {}
for label, (method, options) in calculations.items():
    results[label] = partition_periodic(
        method,
        coordinates,
        numbers,
        grid,
        density,
        density_cutoff=1.0e-14,
        threshold=1.0e-7,
        maxiter=1000,
        **options,
    )

print(f"{'method':12s} {'q(H1)':>10s} {'q(H2)':>10s} {'iterations':>11s} {'solver':>10s}")
for label, result in results.items():
    print(
        f"{label:12s} {result.charges[0]:+10.6f} {result.charges[1]:+10.6f} "
        f"{result.iterations:11d} {result.solver:>10s}"
    )

method            q(H1)      q(H2)  iterations     solver
Hirshfeld     -0.202526  +0.202532           0      fixed
Hirshfeld-I   -0.254761  +0.254762           8         sc
MBIS          -0.193535  +0.193536          30  optimizer
LISA          -0.200000  +0.200000          27  optimizer
AVH-B         -0.196296  +0.196297          18  optimizer

Check every partition

A converged calculation should conserve the total charge, reconstruct the supplied density from the AIM densities, and satisfy the pointwise partition of unity wherever the promolecular density is nonzero.

for label, result in results.items():
    active = result.promolecule > 1.0e-15
    unity_error = np.max(
        np.abs(result.aim_weights[:, active].sum(axis=0) - 1.0)
    )
    assert result.converged
    assert abs(result.charges.sum()) < 1.0e-5
    assert unity_error < 1.0e-12
    assert result.reconstruction_error < 1.0e-12
print("All five partitions passed the conservation checks.")

All five partitions passed the conservation checks.

Optional self-consistent solver

LISA, AVH, and MBIS also accept solver="sc". The optimizer remains the default because relative performance and convergence depend on the material and basis. The two LISA routes agree for this well-conditioned example.

lisa_sc = partition_periodic(
    "lisa",
    coordinates,
    numbers,
    grid,
    density,
    basis=lisa_basis,
    solver="sc",
    threshold=1.0e-7,
    maxiter=1000,
)
print("optimizer charges:", results["LISA"].charges)
print("SC charges:       ", lisa_sc.charges)
print("maximum difference:", np.max(np.abs(results["LISA"].charges - lisa_sc.charges)))

optimizer charges: [-0.19999981  0.20000009]
SC charges:        [-0.19999983  0.20000011]
maximum difference: 1.6057171592365194e-08

The charge magnitude varies slightly because each method defines a different pro-atom model. This dependence is expected: AIM charges are analysis quantities rather than unique observables. For a scientific comparison, keep the density, quadrature, atomic-density libraries, convergence thresholds, and charge reference fixed across methods.