VISION.md

Vision

Purpose

Imago is an all-electron electronic structure code that implements the OLCAO (Orthogonalized Linear Combination of Atomic Orbitals) method. It uses periodic boundary conditions and an LCAO basis to compute electronic properties of a broad range of material systems: crystals, amorphous solids, nanoparticles, molecules, interfaces, and grain boundaries. It is an academic code used for research and student training.

The current development focus comprises five concurrent prongs. Two of them address Brillouin-zone integration: implementing the Linear Analytic Tetrahedral (LAT) method, and correcting the treatment of eigenvector-dependent quantities under IBZ symmetry reduction. The third addresses SCF startup: augmenting today's isolated-atom initial potential database with potentials drawn from converged solid-state calculations, opening a path toward future environment-aware prediction of the starting potential. The fourth is infrastructural: a shared system for submitting, tracking, and harvesting batches of Imago calculations on an HPC cluster, supporting both near-term needs (the database build, convergence sweeps, validation harnesses) and longer-term goals (ab-initio molecular dynamics, high-throughput screening). The fifth is an accumulation prong that knits the previous two together: a curated historical-guidance database that records which convergence settings worked on which families of systems, so that future calculations can predict a converged operating point from past experience and verify it with a small targeted flight instead of an exhaustive grid search.

Goals

1. Implement LAT k-point integration. Offer the Linear Analytic Tetrahedral method (Bloechl, Jepsen, & Andersen, PRB 49, 16223, 1994) for Brillouin-zone integration alongside the existing Gaussian-broadening approach, selected per calculation rather than replacing it. LAT provides exact analytic integration within each tetrahedron, eliminating the arbitrary broadening parameter and improving accuracy at lower k-point densities.

   Both methods are kept, deliberately. Every quantity that has gained a LAT path retains its Gaussian one, and the choice is a user input. Two reasons, and the second is the one that is easy to lose sight of: existing work stays reproducible, and the two schemes fail differently, so a disagreement between them at the same mesh is diagnostic. Retiring the Gaussian path would discard that check. Where a quantity has no LAT path yet, that is unfinished work rather than a decision against it.

2. Correct partial properties under IBZ reduction. Ensure that eigenvector-dependent quantities (PDOS, bond order, effective charge, and the partial optical properties) are computed correctly when the k-point mesh is reduced to the irreducible Brillouin zone. The current approach of multiplying IBZ contributions by star multiplicity is incorrect for these quantities.

3. Augment the initial SCF potential database. Extend today's isolated-atom potential database with potentials taken from converged Imago calculations on a curated set of real solids (drawn from the Crystallography Open Database and/or the Materials Project). The original single-atom potentials are preserved -- both as a fallback when no improved potential is available, and as the baseline against which improvements are measured. The first deliverable is one improved potential per element. The database grows progressively, not exhaustively: new (element, environment) entries are added as evidence warrants, with no attempt to cover every combination. The schema and lookup interface must be designed from the outset to permit future keying by additional attributes -- local environment, expected charge state, coordination -- without disturbing callers. Two payoffs: faster SCF convergence on typical systems, and the ability to run useful non-SCF (single-pass) calculations whose results remain informative for many purposes. A further payoff, realized once the seed exists: the database -- converged reference potentials plus the run records behind them -- becomes the substrate for systematically testing exchange-correlation functionals. Wigner LDA is the most thoroughly tested functional in OLCAO and serves as the trusted baseline; the same flight machinery can later re-run the curated solids under other functionals and weigh their convergence and quality against that baseline. A future integration of an external exchange-correlation library (libxc / xclib), should it prove feasible, would widen the functional set with higher confidence in each functional's correctness.

4. Enable high-throughput Imago calculation flights. Provide shared infrastructure for submitting, tracking, and harvesting batches of Imago calculations on an HPC cluster, so that scripts whose purpose is what to compute and what to learn -- initial-potential database build, convergence sweeps, validation harnesses, future ab-initio molecular dynamics, future high-throughput screening -- do not each reinvent how to submit and watch jobs. The infrastructure rests on three named layers. imago.py continues to run one Imago calculation per invocation, refactored to expose a callable Python API alongside its existing CLI entry point. An ASE Calculator wraps imago.py and adapts imago to the broader computational-materials community; committing to ASE yields LAMMPS interoperability, ab-initio molecular dynamics through ASE's MD integrators, trivial convergence-test loops, and acceptance by every downstream workflow tool that consumes ASE -- essentially for free. Parsl orchestrates calculations on SLURM as the v1 dispatch backend; its Python-as-workflow idiom handles both embarrassingly parallel sweeps (thousands of independent jobs) and tightly iterative inner loops (AIMD, adaptive convergence) under one mental model. Snakemake may join later as an outer orchestration layer if DAG-of-codes workflows -- e.g., SCF -> bands -> optics chains spanning many structures -- become a regular need; in that scenario the inner Parsl logic is unchanged and Snakemake is additive rather than a replacement. This dispatch must be portable and low-ceremony. Imago will run on clusters operated by people outside the original group, so a user must never hand-write SLURM batch scripts: they describe their site once -- its queues, account, per-node cores and accelerators, and how a worker brings up the imago environment -- and choose a few per-run knobs, and the flight layer generates the scheduler submission, filling in defaults and surfacing suggestions for the handful of genuinely site- or run-specific decisions. Two dispatch shapes must both be available, because they suit opposite regimes of the inter-problem axis (Goal 7): a single shared allocation that many small, similar calculations stream through (the convergence sweeps and database seeds), and one scheduler job per calculation for large or heterogeneous runs (the future MPI/GPU solves of differing size). Neither needs nested batch files -- the nesting is of resources within an allocation, not of submissions -- and the per- calculation resource request (cores, ranks, threads, devices, walltime) is exactly what the resource-and-cost dataspace (Goal 6) predicts, so that sizing too becomes automatic rather than hand-guessed.

5. Accumulate historical convergence guidance as a chemistry-and-physics dataspace. Build a curated dataspace that records, for each converged calculation imago has run, the chemistry-and-structure signature of the system (element-group composition, lattice family, system type), the resulting electronic-structure character (band gap, spin polarization), and the convergence settings that worked (initially: the k-point density in points per unit reciprocal-cell volume; later: cell size, basis size, and other knobs). A small regression / nearest-neighbor predictor learns from this dataspace: given a new system's signature, it predicts the converged operating point and an uncertainty. New calculations then run a verification grid around the predicted point -- narrow when the predictor is confident, wide when it is not -- rather than scanning a convergence surface from scratch. Each successful flight appends its result back into the dataspace (via a staging step that a curator promotes), so the artifact grows monotonically with use and the predictor gets better as it accumulates evidence. The dataspace makes the initial-potential-database build (Goal 3) faster because the producer becomes a predict-then-verify client of this prong rather than a guesser of settings, and it makes routine convergence work on new systems vastly cheaper than today's full-grid sweeps.

6. Accumulate resource-and-cost guidance as a hardware-aware dataspace. Build a second curated dataspace -- a sibling of Goal 5, sharing its library scaffolding but kept deliberately separate -- that records, for every imago run, how much it cost: the problem-size signature (atom count, electron count with the core/valence split that an orthogonalized method cares about, basis-function count, the wave-function representation -- 1-component Schrodinger versus 4-component Dirac, which multiplies the secular dimension -- k-point count, and spin channels), the execution configuration it ran under (cores, nodes, cores-per-node, MPI ranks, OpenMP threads per rank, core/socket pinning -- captured as an extensible, registry-validated set so new knobs like GPUs or NUMA policies append rather than break the schema), the build and toolchain it was compiled with (compiler family and version, optimization level, and how the key libraries -- HDF5, ScaLAPACK, BLAS/LAPACK, MPI -- were built, e.g. parallel versus serial or threaded versus sequential), recorded at a deliberately coarse granularity so build choices become a comparable feature rather than a source of fragmentation, while the full compile string is retained as provenance so a targeted study of one specific flag can still recover it; and the measured resources (peak memory, disk footprint, walltime, and optionally per-phase timings). The atomic record is a single execution observation -- one run under one configuration -- never collapsed to a per-system summary, so the same artifact answers several questions. Crucially, unlike the convergence dataspace, this data is hardware-specific and not portable across machines: a walltime measured on one cluster is misinformation on another. The dataspace is therefore partitioned by a hardware fingerprint, and a physics-informed regressor (memory grows like the square of the secular dimension, the eigensolve like its cube) learns the cost surface within each fingerprint. The near-term consumer is provisioning: given a chosen parallel configuration, predict the memory, disk, and walltime a run will need so the flight layer fills SLURM resource requests that neither overflow memory nor hit the walltime limit -- replacing today's hand-guessed requests and advancing the cost discipline of Principle 6. Because every observation records its full configuration, the same artifact later supports recommending the most efficient parallel configuration and pure scaling studies, with no schema change. Like Goal 5, it grows monotonically through a staging step that a curator promotes.

7. Parallelize across scales and devices, keeping every axis open. Imago must scale along two orthogonal axes that serve different needs and must never be conflated. Intra-problem parallelism distributes a single large secular problem -- thousands or more basis functions -- across many MPI ranks and nodes, so that a system too large for one node's memory still becomes tractable. Inter-problem parallelism runs many independent calculations at once; this is the kaleidoscope flight axis of Goal 4, and the two axes trade against each other on a fixed core budget. A third, cross-cutting concern is device placement: some phases of imago (dense integral assembly, the eigenvalue solve) may run best on GPU while others (branchy setup, I/O orchestration, latency-bound logic) stay on CPU. The architecture must let each phase run where it performs best rather than forcing one global device mode. The near-term deliverable follows the measured cost structure (ARCHITECTURE 6.8, 2026-08-18): distribute the three-centre integral stage and the secular solve -- MPI with the ELPA eigensolver -- while the inter-problem axis of Goal 4 continues alongside. The real-space-grid work measures at a few percent of a run and comes last, and GPU offload is a staged research target, validated one phase at a time, not a day-one commitment. A sibling branch (upolcao) already explored MPI parallelization of the common-ancestor OLCAO code; its design is mined as reference (ARCHITECTURE 6.7), not merged, because the two codebases diverged from that shared ancestor and have drifted file by file. The parallel path is the default, not an option bolted on: the standard build links MPI and the ELPA eigensolver, and a serial build is kept only as an explicit fallback for a host that lacks the parallel libraries. This is affordable because the environments Imago targets for teaching and research -- a provisioned Ubuntu virtual machine or a teaching cluster -- already carry the parallel stack, so a student never has to assemble it by hand and serial is the exception rather than the common case. The nanoHUB deployment is a separate, deferred question: its own parallelization must be investigated before its packaging is settled, and it does not constrain this default.

Design Principles

1. Eigenvalue symmetry != eigenvector symmetry. Point-group operations guarantee e_n(k) = e_n(Rk), but eigenvectors at Rk are related to those at k by a basis transformation that mixes orbital coefficients. Any integration scheme that weights eigenvector-dependent quantities by IBZ star multiplicity alone is incorrect. This was confirmed empirically: bond orders for symmetry-equivalent O atoms in KNbO3 disagreed under IBZ reduction but agreed with the full mesh.
2. Full mesh for post-processing. Use the IBZ for the expensive SCF diagonalization, but reconstruct the full BZ mesh before computing partial properties. For eigenvalue-only quantities (TDOS), unfolding eigenvalues by symmetry suffices. For eigenvector-dependent quantities, either use the full mesh directly or apply atom-permutation corrections.
3. Phased implementation. Implement LAT TDOS first (eigenvalues only), validate against Gaussian-broadened results at high k-point density, then extend to integrated partial properties via electronPopulation_LAT, then to energy-resolved PDOS with cornerIntgWt_LAT, and then to the optical properties, whose integrand is a joint density of states over band PAIRS rather than a density of states over bands. Each phase must be validated before proceeding to the next.
4. Stable potential representation, flexible container. The per-potential numerical representation (the existing radial / Gaussian-fit form) is preserved. How potentials are stored, labeled, and combined into a starting potential is open to redesign -- one file per (element, label) pair, one file per element containing many labeled potentials, or other arrangements as needed. Downstream consumers see only the assembled starting potential, not how it was retrieved or composed, so the labeling scheme and retrieval logic can grow more sophisticated (e.g., interpolation across nearby labels when no exact match exists) without disturbing callers.
5. Curation and regeneration are both part of the deliverable. The value of the improved-potential prong depends directly on the breadth and quality of the reference solids used to generate the database. A documented curation strategy -- which structures, why, and how to extend the set as new chemistries are encountered -- is required from the start, not added as an afterthought. The database itself must be regeneratable from the curated set, not a hand-edited artifact, so that adding entries or rebuilding after a code change is a scriptable operation rather than a manual exercise.
6. Cost discipline at SCF setup. Any prediction or lookup performed at SCF initialization must remain cheap relative to the SCF iterations it is meant to shorten or replace. This bounds how elaborate the environment-keyed prediction step can grow, and guards against over-engineering the prediction layer for marginal accuracy gains.
7. Iteration count is the headline metric. Every reference calculation that contributes to the database records its SCF iteration count. Across a representative benchmark, runs that start from improved initial potentials should take at least 20% fewer iterations on average than runs starting from the original isolated-atom potentials. The original potentials are kept in the database precisely so this comparison can be made on the same code path, with only the starting potential changed.
8. Decouple adapter from orchestrator. The ASE Calculator adapts imago to the materials-simulation community; Parsl orchestrates calculations on the cluster. These two layers are independent. Replacing either should not require changing the other, and new adapters (other ASE-compatible calculators) or new orchestrators (Snakemake on top, or a different dispatch backend) can slot in along well-defined seams.
9. General-purpose orchestration. Workflow tools are chosen so that the skills students acquire while developing or running Imago transfer beyond computational materials. Domain-specific machinery lives at the adapter layer (ASE); the flight layer (Parsl) is ordinary scientific Python. A student who leaves the field carries portable workflow skills rather than Imago-specific muscle memory.
10. Complete-and-report at the flight level. A single failed calculation never fails an entire flight by default. Per-job failures are recorded and surfaced; the flight driver script decides whether the aggregate outcome is acceptable for its scientific purpose. This applies whether the cause is convergence non-attainment, cluster-side job loss, or post-processing error.
11. Experience as a curated artifact, not tribal knowledge. Hard-won operating knowledge -- which k-point density a family of oxides converges at, which supercell size a defect chemistry needs, which initial potential worked for which coordination environment -- belongs in a checked-in, schema- governed, regeneratable database, not embedded in scripts or held in researchers' heads. Tools consult that database to predict; flights verify the prediction; the verified result feeds back into the database through a deliberate curation step. This applies equally to the initial- potential database (Goal 3) and the historical-guidance database (Goal 5); both share the same curate/regenerate/ consult discipline. The implication for code: scripts must never silently encode "experience" as hardcoded constants; anything tunable enough to need empirical justification belongs in a database that future readers can audit and future flights can update.
12. The flight layer stays dumb; description lives in Python. The kaleidoscope dispatcher is thin: it drives a flat list of independent CalcUnit records. It does not interpret per-calculation options, has no notion of "this calc depends on that calc," and does not grow a flight description language (no flight DSL, no YAML/TOML workflow grammar, no DAG engine). Higher-order flight shape -- multi-axis sweeps, dependent phases, per-unit internal iteration -- lives in ordinary client- side Python that builds the flat list before handing it off, or in custom wingbeats that own the iteration for one unit. This is the deliberate design choice that distinguishes kaleidoscope from workflow engines like Snakemake or LAMMPS' input language: by refusing to grow a domain language, kaleidoscope stays open to every flight shape its clients invent, and the skills a student picks up while building flights are ordinary Python rather than a one-off DSL. Principle 9 says "domain-specific machinery lives at the adapter layer"; Principle 12 says specifically what the flight layer refuses to absorb.
13. Two axes of parallelism, never conflated. Intra- problem (one large solve spread across many ranks) and inter-problem (many small solves, one rank-set each) are different problems with different bottlenecks and different ideal hardware shapes. A design that serves one must not silently assume the other. The resource- and-cost dataspace (Goal 6) records both ranks-per-run and runs-per-flight precisely because they trade against each other on a fixed core budget, and the flight layer must be able to request either shape.
14. Device placement is empirical and per-kernel. Whether a phase belongs on CPU or GPU is decided by measurement, not assumption. The cost of moving data across the CPU-GPU boundary can erase a kernel's on-device speedup, so a phase that is dense, regular, and single-precision-tolerant may win on GPU while a branchy, memory-bound, or latency-sensitive phase does not. The architecture exposes a placement boundary so individual kernels move between devices without rewrites, and the decision for each kernel is backed by the cost dataspace (Goal 6) rather than asserted up front.
15. Attribution is a non-negotiable, not a formality. Imago cites the work it derives from, and it is built so that credit survives being passed on. Two obligations follow from that. Where an algorithm or dataset comes from published work, the source is cited in the code beside the thing it explains, so that a student reading cold can find the paper rather than guess at it. And where Imago itself is redistributed or built upon, its own credit travels with it: the license requires that NOTICE accompany any derivative work, which is why upstream sources belong there rather than in a comment a fork can quietly drop. Copyright ownership and authorship credit are distinct, and are recorded separately -- a uniform copyright holder in every file header, and per-author credit in the file itself and in CITATION.cff. Conflating the two is the common failure mode. This is the same obligation that governs a methods section, and the source is where students first learn that it applies to software as well as to prose.
16. Where the answer is a judgment, report it and let the curator decide. A tool that cannot settle a question from the data in front of it prints what it found, together with the facts a person needs in order to choose, and leaves the choosing to them. It does not encode a policy that guesses on their behalf. Principle 10 says a failed calculation is recorded and surfaced rather than failing a whole flight; this says the same of an ambiguity -- a stored result whose engine build has since moved on, a harvested record claiming a structure the database already holds. Each of those is a question with a human answer, and a branch that settles it silently is a decision nobody reviewed. What makes the unhurried answer affordable is what these databases actually hold. An initial potential is a starting point that every downstream SCF re-converges; a predicted k-point density is a starting point that the climb then verifies against its own two-sided test. A stale or imperfect stored value therefore costs iterations, not correctness -- which is exactly the condition under which elaborate automatic safeguards stop paying for themselves (Principle 6). Treating either database as though its numbers were final results buys caution nobody needed and complexity every later reader must still read. The implication for code: prefer a printed line to a branch; make the expensive action explicit, so that a re-run is requested rather than inferred; and record the facts a judgment would want -- which build produced a number, which mesh a record claims -- even where nothing compares them.