Key Concepts

Key Concepts#

This section describes the core architecture and data formats used by AimsPy. Understanding these concepts is essential for using the API effectively and for extending the package with new functionality.

Overview#

AimsPy drives FHI-aims DFT calculations directly from Python by loading a patched libaims.so via ctypes. There is no subprocess and no file-staged I/O on hot paths — Hamiltonian, overlap, energy, and forces are exchanged as in-memory arrays through a callback framework exposed by the bundled FHI-aims patch.

The package is organised in three layers:

Layer

Purpose

Public?

aimspy.calculator, aimspy.matrix, aimspy.structure, aimspy.interface.deeph

User-facing API

Yes

aimspy._callbacks

Callback framework (CallbackManager, hook-point specs)

No

aimspy._binding, aimspy._patches

ctypes binding to libaims.so; versioned FHI-aims diffs

No

The central user-facing class is Calculator, whose lifecycle is governed by the CalcState state machine.

In-Memory Architecture#

The ctypes binding layer#

aimspy._binding.libloader.load_aims_lib(lib_path) is the only place that calls ctypes.CDLL on libaims.so. Two details matter:

  1. MPICH symbol visibility. mpi4py’s own shared object must be loaded with RTLD_GLOBAL and anchored at module level (_mpi_cdll_anchor). Without the module-level anchor, the CDLL would be garbage-collected (and dlclose’d) when the function returns, removing the global symbols that libaims.so needs for lazy MPI symbol resolution. AimsPy loads libaims.so itself with RTLD_GLOBAL as well.

  2. Forward-compatible symbol probing. BindingLib (in aimspy._binding.prototypes) remembers which C symbols were detected and exposes a has(name) predicate. setup_prototypes silently skips symbols missing in older libaims builds, so a single AimsPy release can drive multiple patch versions.

The Calculator lifecycle#

CalcState is a six-state enum with a directed lifecycle:

UNINIT ──init()──> INITED ──calc()──> [RUNNING] ──> DONE
                     │                                │
                     └──close()───────────────────────┘
                                                │
              close()/force_close() ──────────> FINALIZED
   (any state, on error) ──> FAILED ──force_close()──> FINALIZED

State

Meaning

Allowed next

UNINIT

Freshly constructed

init / do / force_close

INITED

aimspy_init done; callbacks registered; info/structure readable

calc / close

RUNNING

Transient inside calc() while aimspy_run is executing

DONE / FAILED

DONE

SCF converged; hamiltonian/forces readable

close

FAILED

Operation aborted; Fortran runtime in unknown state

force_close only

FINALIZED

aimspy_finalize called; terminal

—

State transitions are guarded — calling calc() from UNINIT, or close() from RUNNING, raises AimspyStateError. Use force_close() for safe recovery from any state.

Note: energy is also accessible in the RUNNING state (e.g. from inside a callback during pre-SCF), but the value may be uninitialized before the first SCF iteration completes.

Thread safety: Calculator is not thread-safe. FHI-aims uses chdir(2) (process-global) internally, so concurrent Calculator operations in the same process will race. The typical MPI usage is one rank = one process.

Callback Framework#

The bundled FHI-aims patch inserts trigger points across the FHI-aims driver files — mainly inside src/initialize_scf.f90 (after reshape_matrices, before the initial diagonalisation), but also in src/prepare_scf.f90 (pre-SCF basis export), src/scf_solver.f90 (post-SCF grid export), and src/DFPT_main/DFPT_module.f90 (pre/post-CPSCF first-order Hamiltonian export and injection). The available callbacks are:

Callback

Purpose

get_descr

Capture the CSR sparse-storage layout

export_ovlp

Export the overlap matrix

export_h0

Export the free-atom initial Hamiltonian (H_init)

python_func

Generic Python hook (deferred source generation)

modify_h0

Inject the modified H_init back into FHI-aims

export_dHde

Export DFPT first-order Hamiltonian (post-CPSCF)

modify_dHde

Inject modified first-order Hamiltonian (pre-CPSCF)

export_grid_data

Export real-space grid data (post-SCF convergence)

export_basis_data

Export NAO radial basis splines (inside prepare_scf, during aimspy_init — must be registered before aimspy_init)

When modify_h0 is registered, the patch short-circuits the standard initial diagonalisation: it calls advance_KS_solution directly on the injected Hamiltonian and sets restart_zero_iteration=.true., which is what enables warmstart in several iterations.

Note: Matrix extraction and injection (overlap, Hamiltonian, H_init) require a periodic system with use_local_index = .false.. Forward SCF calculations (without matrix extraction/injection) work with any system type. For isolated molecules, use a sufficiently large periodic cell with vacuum.

How callbacks are wired#

CalculatorConfig flags control which default callbacks are auto-registered by Calculator._wire_callbacks:

Flag

Enables

Effect

(always)

get_descr

calc.csr_descr populated

capture_overlap=True

export_ovlp

calc.overlap returns live overlap on all ranks

capture_initial_hamiltonian=True

export_h0

calc.initial_hamiltonian populated

modify_init_ham(...) called

python_func + modify_h0

warmstart / scaling / custom modification

capture_first_order_hamiltonian=True

export_dHde

calc.first_order_hamiltonian = [x, y, z] matrices (DFPT)

modify_init_first_order_ham(...) called

modify_dHde

dH/de warmstart injection (pre-CPSCF)

capture_grid_data=True

export_grid_data

calc.grid_data populated with this rank’s grid subset

capture_basis_data=True

export_basis_data

calc.basis_data populated (registered pre-aimspy_init because it fires during init)

Users can also register custom callbacks via Calculator.register_callback(name, fn, aux, extra_ptr) for advanced use cases.

Error handling#

A Python exception raised inside a callback never crashes Fortran. The CallbackManager records (name, exception, traceback_str) tuples and lets calc() complete. After aimspy_run returns, Calculator._check_callback_errors raises a single AimspyCallbackError aggregating all failures, with the per-callback details preserved on exc.callback_errors. Notably, forces are captured before the callback error check, so they survive even when a callback raises.

Hamiltonian Modification Strategies#

Calculator.modify_init_ham(source=..., strategy=...) configures how the live H_init buffer is mutated before SCF starts. The built-in strategies, dispatched by the pure function _apply_strategy:

Strategy

Behaviour

Required argument

Typical use

REPLACE

Clear and copy external blocks into the live H_init

source=

Warmstart with a DeepH prediction

ADD

Add external blocks on top of the live H_init

source=

Correction (Delta-prediction): add predicted H − H₀ to H₀ to recover H

SCALE

Multiply the live H_init by a constant factor

factor= (float)

Scaling experiments

CUSTOM

Call custom_fn(live, external, structure, aux) to mutate the live matrix in place

custom_fn= (callable)

Arbitrary transforms

Direct vs. deferred source#

modify_init_ham supports two modes:

  • Direct mode (source= is passed): the source object is stored immediately. The python_func callback then converts it via source.to_aimspy(structure) and stores the result in _runtime_aux["external_aimspy"], which modify_h0 reads.

  • Deferred mode (used as a decorator): the user function is called during the python_func callback, after export_h0 and export_ovlp have fired, so it has live access to calculator.initial_hamiltonian and calculator.overlap. This is essential when the external source needs the runtime structure to be built.

# Direct
calc.modify_init_ham(source=data, strategy=Strategy.REPLACE)

# Deferred
@calc.modify_init_ham(strategy=Strategy.REPLACE, option={"path": "deeph_out/"})
def gen_source(calculator, option):
    return DeepHData.from_directory(option["path"])

ExternalMatrixSource Protocol#

The ExternalMatrixSource Protocol (aimspy.interface) is the contract for any external matrix provider:

from typing import Protocol, runtime_checkable

@runtime_checkable
class ExternalMatrixSource(Protocol):
    def to_aimspy(self, structure: AimspyStructure) -> AimspyMatrix: ...

The reference implementation is DeepHData (importable from aimspy directly), which reads the DeepH on-disk format (POSCAR + info.json + .h5) and converts to AimspyMatrix. Adding a new external format is a single subpackage under aimspy/interface/<format>/ with a class satisfying this protocol — no other code changes are required.

AimspyMatrix Block-Sparse Format#

AimspyMatrix is AimsPy’s canonical in-memory representation of block-sparse real-space matrices. It holds a dict mapping atom-pair keys to dense numpy blocks:

from aimspy import AimspyMatrix
import numpy as np

# Direct construction
matrix = AimspyMatrix(
    blocks={
        (0, 0, 0, 0, 0): np.array([[1.0, 0.5], [0.5, 1.0]]),  # R=0, atom 0→0
        (1, 0, 0, 0, 1): np.array([[0.1]]),                     # R=(1,0,0), atom 0→1
    },
    n_spin=1,
)

# Access
block = matrix.blocks[(0, 0, 0, 0, 0)]  # ndarray, shape (2, 2)
print(matrix.n_pairs)                    # 2

Each key is a 5-tuple (R1, R2, R3, i_atom, j_atom):

  • R1, R2, R3 — lattice vector components (integers), following R_aimspy = -R_aims = R_deeph

  • i_atom, j_atom — 0-based atom indices in aims native order

Each value is an np.ndarray of shape (n_orb_i, n_orb_j), dtype float64, where n_orb_i / n_orb_j are the number of basis functions on atom i / j.

Conventions#

Property

Convention

R vector

R_aimspy = -R_aims (same sign as DeepH)

Atom indices

aims native order (no reordering)

Orbital order

aims native basis order (no reordering)

Parity

wiki/DeepH convention (phase_i * phase_j already applied)

Units

Hartree (Hamiltonian), dimensionless (overlap)

Hermitian partners

both (R,i,j) and (-R,j,i) are stored

Phase factor#

The parity convention is implemented in AimspyStructure.phase_factor:

phase_factor = np.where((basis_m > 0) & (basis_m % 2 == 1), -1, 1).astype(np.int32)

It is self-inverse (phase² = 1), so applying it once converts aims ↔ aimspy and applying it again undoes the conversion.

Construction#

Three ways to obtain an AimspyMatrix:

# 1. Direct (offline, for testing or custom data)
matrix = AimspyMatrix(blocks={...}, n_spin=1)

# 2. From FHI-aims CSR (after calc.do(), rank 0 only)
H = calc.hamiltonian  # AimspyMatrix

# 3. From DeepH on-disk format
from aimspy import DeepHData
data = DeepHData.from_directory("deeph_out/")
H = data.to_aimspy(calc.structure)  # AimspyMatrix

Conversions require an AimspyStructure (provides atom/orbital info and derived properties like phase_factor, orbit_per_atom, atom_permutation) and a CsrMatrixDescriptor (FHI-aims’ CSR sparse layout — see API Reference for field details).

Conversion#

  • To aims CSR: matrix.to_aims_csr(csr_descr, structure) returns a (n_spin, n_ham_size) C-contiguous array, ready for ctypes.memmove into the Fortran buffer. Hermitian fallback: if (R,i,j) is missing, (-R,j,i) is used with transposition.

  • To DeepH format: DeepHData.from_aimspy(structure, hamiltonian=matrix, ...) handles atom reordering (aims → POSCAR) and unit conversion (Hartree → eV).

Limitations#

  • Spinless only: converters read/write spin channel 0; n_spin=2 leaves channel 1 as zero. Spin-polarised support is on the roadmap.

  • Periodic only: matrix extraction/injection requires use_local_index = .false. (see Troubleshooting).

DeepH Data Format#

DeepHData reads and writes the standard DeepH on-disk format used throughout the DeepH ecosystem. The format is shared with DeepH-dock — for the full field-level specification, see the DeepH-dock Key Concepts page. A summary:

some_directory/
├── POSCAR                   # Atomic structure (VASP format, element-grouped order)
├── info.json                # System metadata + basis set info
├── overlap.h5               # Overlap matrix S (sparse)
├── hamiltonian.h5           # Hamiltonian H (sparse, eV)
├── hamiltonian_init.h5      # Free-atom initial Hamiltonian (sparse, eV)
├── force.h5                 # (optional) Forces + energy (MD-style: cell/energy/force/stress)
└── electric_response.h5     # (optional) DFPT first-order Hamiltonian dH/de (sparse, eV)

Each matrix .h5 file stores four datasets: atom_pairs (N,5), chunk_boundaries (N+1,), chunk_shapes (N,2), and entries (M,). The atom order in POSCAR is element-grouped (different from aims native order); DeepHData handles the reordering via AimspyStructure.atom_permutation.

force.h5 uses a different MD-style layout: cell (3,3), energy (scalar), force (n_atoms,3), stress (6,) (zeros placeholder), with formula and natoms root attributes. Forces are in eV/Å (matching calc.forces), energy is in eV (converted from calc.energy Hartree).

electric_response.h5 stores the DFPT first-order Hamiltonian (dH/de) — the response of the Hamiltonian to an electric field perturbation. It uses the same atom_pairs as hamiltonian.h5, but chunk_shapes rows are 3× larger (one block per Cartesian direction [y, z, x] = real spherical harmonics m = -1, 0, +1), and entries is 3× longer. Units are eV (converted from Hartree). Requires electric_field_response DFPT + electric_field_serial .false. in control.in.

Unit conventions#

Quantity

AimsPy internal

DeepH on-disk

Hamiltonian

Hartree

eV

Overlap

dimensionless

dimensionless

Coordinates

Å

Å (in POSCAR)

Force

eV/Å (calc.forces)

eV/Å (force.h5)

Total energy

Hartree (calc.energy)

eV (force.h5 energy dataset)

First-order H (dH/de)

Hartree

eV (electric_response.h5)

DeepHData.from_memory converts Hartree → eV on write; DeepHData.to_aimspy converts eV → Hartree on read. Force requires no unit conversion (eV/Å throughout); only atom reordering is applied. Energy in force.h5 is converted from Hartree via from_aimspy(force=, energy=) or set_force(force, structure, energy=). First-order Hamiltonian is converted from Hartree via from_aimspy(first_order_hamiltonian=) or set_first_order_hamiltonian(), and back to Hartree via to_first_order_aimspy().

FHI-aims Patch System#

Three patch versions are bundled (v0.1.0 ~1100 lines, v0.2.0 ~1400 lines, v0.2.1 ~2200 lines — the latest, adding the grid/basis exports and the DFPT dH/de hooks). The patch does three things:

  1. Adds src/aimspy_api/ with Fortran modules:

    • callback.f90 — TAimspyCsrMxDescr (bind(C) struct), TAimspyCallback handle type, abstract callback interfaces.

    • api_bank.f90 — module-level save arrays (c_hamiltonian, c_overlap), aimspy_energy, aimspy_forces accessors.

    • info.f90 — TAimspyInfo bind(C) struct + aimspy_get_info populating a save buffer.

    • register.f90 — the aimspy_register_*_callback bind(C) subroutines.

    • main.f90 — aimspy_init / aimspy_run / aimspy_finalize / aimspy_all lifecycle entry points.

    • export_grid_data.f90 / export_basis_data.f90 — module-level buffer assemblies + triggers for the post-SCF grid export and the pre-SCF NAO basis export (both finalized in aimspy_finalize to release their buffers).

  2. Hooks into src/initialize_scf.f90 — trigger points after reshape_matrices, and the warmstart short-circuit calling advance_KS_solution on the injected Hamiltonian with restart_zero_iteration=.true..

  3. Exposes pbc_lists.f90 arrays — adds target attributes to index_hamiltonian / column_index_hamiltonian so they can be exposed via c_loc.

The patch is versioned (currently v0.2.1) and managed by the aimspy patch CLI, which can apply, uninstall, dry-run, and list bundled versions. Multiple patch versions can ship side-by-side; the CLI auto-detects the currently-applied version by reading a PATCH_VERSION line that the patch itself writes into the source tree’s Makefile.

Grid Data (Real-Space)#

GridData is AimsPy’s in-memory representation of the FHI-aims real-space integration grid and the scalar fields living on it. It is captured after SCF convergence via the export_grid_data callback.

Fields#

Field

Shape

Units

Description

coords

(3, n)

bohr

Grid point coordinates (mapped to center cell for periodic)

partition_tab

(n,)

bohr³

Grid point integration weights

index_atom

(n,)

—

0-based atom index for each point

index_radial

(n,)

—

Radial shell index

index_angular

(n,)

—

Angular grid index

rho

(n_spin, n)

e/bohr³

Converged electron density

vks

(n_spin, n)

Hartree

Kohn-Sham potential (includes vdW if active)

vks0

(n_spin, n)

Hartree

Free-atom reference potential (no vdW)

vh

(n,)

Hartree

Hartree potential (includes nuclear attraction)

vh0

(n,)

Hartree

Free-atom Hartree potential

rho0

(n,)

e/bohr³

Free-atom superposition density (rho0 == rho_free)

atom_coords

(n_atoms, 3)

Å

Atomic coordinates (from in-memory structure)

atom_symbols

(n_atoms,)

—

Element symbols

lattice

(3, 3)

Å

Lattice vectors

Key semantics#

  • rho0 is rho_free: FHI-aims exports free_rho_superpos which carries a factor of 4π; AimsPy normalises at import so rho0 IS the free-atom density.

  • vks includes vdW: When use_vdw_correction_hirshfeld_sc, use_mbd_std, or use_libmbd is active, vks = V_H + V_nuc + v_xc + v_vdw. vks0 does NOT include vdW (free-atom reference has no vdW correction).

  • LDA scalar only: The GGA non-local (vector) term 4*xc_gradient_deriv is NOT exported. vks is exact for LDA, scalar part for GGA. Hybrid functionals are not supported.

Derived quantities#

gd.delta_rho    # rho - rho_free (density difference)
gd.delta_vks    # vks - vks0 (potential difference)
gd.vxc          # vks - vh (exchange-correlation potential)
gd.vxc0         # vks0 - vh0 (free-atom XC potential)
gd.coords_ang   # coords in Å (converted from bohr)
gd.vks_ev       # vks in eV

MPI gather#

GridData.gather(local, comm) collects per-rank subsets to root using mpi4py.MPI.Comm.Gatherv (zero-pickle, memory-efficient). Root peak memory is ~1x the total dataset, compared to ~3x for the default comm.gather on a Python dict.

npz serialization#

gd.save_npz("grid.npz")           # save (includes structure fields if present)
gd2 = GridData.load_npz("grid.npz")  # load

The npz format is self-describing: it stores n_full_points, n_spin, n_atoms, all grid arrays, and optionally atom_coords/atom_symbols/lattice.

NAO Radial Basis (BasisData)#

With capture_basis_data=True, the export_basis_data callback (registered before aimspy_init, because it fires inside prepare_scf during init) captures the complete cubic-spline representation of the NAO radial basis: spline coefficients for u(r), (e−v)·u(r), and du/dr, plus per-species logarithmic grid parameters (r_grid_min, r_grid_inc, n_grid) and per-function outer_radius.

  • Evaluation: basis_data.evaluate_u(i_fn, r) / evaluate_phi / evaluate_du_dr — the per-function 0-based species map is attached automatically at init (basis_data.species_of_fn); u(r) evaluates to zero outside [r_grid_min, outer_radius].

  • Identity metadata (n, l, type, species per radial function) comes from calc.info (basisfn_n / basisfn_l / basisfn_type / basisfn_species).

  • Units: lengths in bohr, energies in Hartree; u(r) normalized so ∫u²dr = 1 (bohr^−1/2).

basis.h5 format#

BasisData.save_h5(path, info) builds an incremental element-per-group library: the file is created if missing; existing element groups are skipped silently (not overwritten), so one file can accumulate basis sets across calculations.

/attrs: format_version, generator, date, units, species_list, n_species
/<El>/attrs: element, z, r_grid_min, r_grid_inc, n_grid,
             n_basis_rad, n_orbitals, l_max
/<El>/r_grid        (n_grid,)                bohr, shared log grid
/<El>/n, l, zeta    (n_basis_rad,) int32     quantum numbers; zeta = index
                                             among same (n,l) duplicates
/<El>/type          (n_basis_rad,) S8        atomic/hydro/...
/<El>/outer_radius  (n_basis_rad,)           bohr, per-function cutoff
/<El>/spline_wave   (n_basis_rad, 4, n_grid) cubic coeffs for u(r)
/<El>/spline_kinetic (same)                  for (e−v)·u(r)
/<El>/spline_deriv  (same)                   for du/dr

Visualize offline with aimspy viz-basis basis.h5 or aimspy.viz_basis.plot_radial_basis.

Data Flow in AimsPy#

A typical warmstart workflow:

  1. Input: FHI-aims control.in + geometry.in in work_dir; an external Hamiltonian source (e.g. a DeepHData directory from a DeepH-trained model).

  2. Processing:

    • Calculator.init loads libaims.so, calls aimspy_init, builds AimspyStructure.

    • _wire_callbacks registers the default callbacks based on CalculatorConfig flags.

    • Calculator.calc calls aimspy_run. Inside FHI-aims, after reshape_matrices:

      • get_descr populates the CSR layout.

      • export_ovlp / export_h0 capture overlap and free-atom H_init (if enabled).

      • python_func converts the external source via to_aimspy(structure).

      • modify_h0 applies the Strategy, writes the result back via memmove, and the patch short-circuits the diagonalisation.

  3. Output: calc.hamiltonian (AimspyMatrix, Hartree), calc.energy (Hartree), calc.forces (eV/Å), and optionally calc.overlap / calc.initial_hamiltonian. These can be exported to DeepH format via DeepHData.from_aimspy(...).save(...).

For more on the API surface, see Basic Usage. For extending AimsPy with new callbacks or matrix sources, see the Development Guide.