Periodic real-space partitioning: quick start

This notebook constructs a small periodic density and partitions it with MBIS. It uses only NumPy, qc-grid, and HORTON-Part, so GPAW and external calculation files are not required. Coordinates, cell vectors, quadrature weights, and densities use atomic units.

from pathlib import Path
from tempfile import TemporaryDirectory

import numpy as np

from horton_part.periodic import partition_periodic
from periodic_setup import make_demo_system, print_diagnostics

Construct a periodic grid and density

The helper places two hydrogen atoms in a cubic cell. One atom contributes 1.2 electrons and the other contributes 0.8 electrons, giving a neutral two-electron cell. Periodic images are accumulated through PeriodicGrid.get_localgrid.

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

atoms: 2
grid points: 4,096
integrated electrons: 1.99999972

Run MBIS

MBIS builds its minimal Slater pro-atoms internally and therefore needs no external atomic-density library. The default optimizer is recommended for production calculations.

result = partition_periodic(
    "mbis",
    coordinates,
    numbers,
    grid,
    density,
    threshold=1.0e-7,
    return_weights=True,
)
print_diagnostics(result)

method: mbis; solver: optimizer
charges: [-0.193535  0.193536]
sum of charges: +2.834e-07 e
maximum partition-of-unity error: 2.220e-16
density reconstruction error: 2.220e-16

Inspect and save the result

charges are obtained by integrating the final AIM densities on the supplied grid. model_charges are the analytic populations of the fitted pro-atoms. Their difference is a useful finite-grid diagnostic. to_dict() returns all arrays needed for an NPZ result archive.

print("integrated charges:", result.charges)
print("model charges:     ", result.model_charges)
with TemporaryDirectory() as directory:
    output = Path(directory) / "mbis.npz"
    np.savez(output, **result.to_dict())
    with np.load(output) as archive:
        print("stored fields:", ", ".join(sorted(archive.files)))

integrated charges: [-0.19353544  0.19353573]
model charges:      [-0.19353563  0.19353541]
stored fields: aim_weights, charges, converged, grid_populations, history, iterations, method, model_charges, model_populations, parameter_counts, parameter_labels, parameters, populations, promolecule, reconstruction_error, solver

Production inputs

For a GPAW calculation, run part-from-gpaw calculation.gpw density.npz in the GPAW environment and then part-periodic density.npz mbis.npz --method mbis. The converter preserves the uniform and PAW augmentation blocks required for an all-electron partition. See the periodic user guide for the complete archive schema and basis requirements.