Document hierarchy: VISION -> ARCHITECTURE -> DESIGN -> PSEUDOCODE -> Code. For goals and principles, see
VISION.md.
src/
imago/ Primary production program
imago.F90 Top-level dispatcher
kpoints.f90 O_KPoints module (mesh, weights, phases)
dos.F90 O_DOS module (TDOS, PDOS)
bond.F90 Bond order (Mulliken overlap population)
populate.F90 Electron population / occupancy
makeKPoints/
makekpoints.F90 Standalone k-point mesh generator
scripts/ Python tools: flat .py files are CLI
entry points; subdirectories are Python
packages
imago.py Imago driver: CLI and callable API (9.2)
makeinput.py Input file orchestrator
makegroups.py Bispectrum type-grouping helper (8.9,
DESIGN 5.10): CLI + importable
matchers.py Matcher protocol library: Reduce/Bispec
matchers + MATCHERS registry (8.9)
ase_imago.py ASE Calculator, ImagoCalculator (9.3)
cod_fish.py COD acquisition front-end -> CIF (9.5)
cif2skl.py CIF -> imago.skl converter (9.5)
kaleidoscope/ Parsl flight dispatcher, "kaleidoscope" (9.4)
kinds.f90 Shared precision kinds
constants.f90 Shared physical/mathematical constants
dev/
VISION.md Goals and principles
ARCHITECTURE.md This document
DESIGN.md Algorithmic design
PSEUDOCODE.md Algorithm specifications
TODO.md Task list by level
LICENSE ECL-2.0, verbatim (1.1)
LICENSES/ Same text as LICENSES/ECL-2.0.txt, the
location the REUSE tooling reads (1.1)
NOTICE Attribution notices that travel with any
redistribution (1.1)
CITATION.cff Machine-readable citation metadata (1.1)
REUSE.toml Blanket licensing for files that cannot
carry a comment header (1.1)
CONTRIBUTING.md Licensing, credit, and chain discipline
README.md Landing page: what Imago is, license,
citation
BUILD.md Production and debugging builds
Imago is licensed under the Educational Community License,
Version 2.0 (ECL-2.0), the same license used by the OLCAO
codebase it succeeds. This section implements VISION principle
15. Six root-level files carry the licensing and attribution
machinery, and they divide the work along a distinction that is
easy to blur: LICENSE states the terms, NOTICE states who
must be credited, and CITATION.cff states how to cite the work
in a paper. REUSE.toml supplies that same information in bulk
for files that have nowhere to put a comment, LICENSES/ holds
a second copy of the license text at the path the REUSE tooling
requires, and CONTRIBUTING.md tells a contributor how to keep
all of them correct.
NOTICE is the load-bearing one. Section 4(d) of the license
requires that any derivative work reproduce its contents, so it
is the mechanism by which credit -- to Imago itself and to
upstream sources such as the Birkbeck College space group tables
used by makeSGDB.py -- survives into work built on top of this
code. Third-party material added to the tree must gain an entry
there.
File header convention. Every shipped source file carries a
two-line header identifying the license and the copyright
holder. This covers Fortran, Python, CMake, shell, and XML, and
it applies wherever the file lives -- not only under src/, but
also the nanoHUB integration in rappture/ and middleware/,
which is distributed just as the engine is:
!! SPDX-License-Identifier: ECL-2.0
!! Copyright (c) 2026 Paul RulisPython, CMake, and shell files use ## in place of !!. The
doubled comment markers mark these lines as structured content
so that rewrap_prose.py will not reflow them. XML takes a
plain <!-- --> comment instead, because the doubling convention
exists solely for the reflow tool and that tool does not process
XML.
Anything the language requires to come first stays first, and
the header slots in beneath it: a #!/usr/bin/env python3
shebang in Python and shell, and the <?xml version="1.0"?>
declaration in XML.
The header is uniform across the tree and does not vary per
file. SPDX-License-Identifier is a machine-readable tag that
license scanners recognize, which makes the tree auditable
without reading it. Per-file authorship is recorded separately,
in an Author: line inside the file's own header block or
docstring, and cumulatively in the authors: list of
CITATION.cff. Copyright ownership and authorship credit are
distinct, and conflating them in the copyright line is the
failure mode this convention exists to prevent.
Verification. The tree is compliant with version 3.3 of the
REUSE specification, which is checked by running reuse lint at
the repository root. It reports every file lacking recognizable
copyright and licensing information and exits non-zero if any
are found, so adding an untagged file breaks a check that
currently passes. A new file either gains a header or gains
coverage from a REUSE.toml annotation.
Why this level. File headers are governed here rather than
in PSEUDOCODE because they carry no algorithmic content: there
are no inputs, outputs, or control flow for a pseudocode section
to specify. The gate requiring pseudocode before editing src/
therefore does not apply to them, and this section is the
specification the header pass was written against.
The runtime half. Everything above is addressed at someone redistributing or packaging Imago. It never reaches the person who ran a calculation and is writing up the results, who will cite whatever the output puts in front of them. That reader is served by the runtime citation banner specified in DESIGN 10, governed by section 12, and the two halves together are what VISION principle 15 asks for.
Modules directly affected by the current development:
-
O_KPoints (
kpoints.f90): Stores the k-point mesh, weights, and phase factors. DeclareskPointIntgCode(0=Gaussian, 1=LAT). LAT data already implemented:numTetrahedra,tetraVol,tetrahedra(:,:),fullKPToIBZKPMap(:),generateTetrahedra,computeTetraVol,initializeKPointMesh(with IBZ folding). Phase F addition for correct IBZ unfolding of eigenvector-dependent quantities:fullKPToIBZOpMap(:), recording which point group operation mapped each full-mesh k-point to its IBZ representative. -
O_AtomicSites (
atomicSites.f90): Holds the atomic site list and, as the other half of the Phase F IBZ unfolding infrastructure, the atom permutation tables:atomPerm(:,:)(the image atom of each site under each point group operation),invAtomPerm(:,:)(its inverse), and thebuildAtomPerm/buildInvAtomPermroutines that construct them from the point operations and the fractional atom positions. The tables live here rather than with the k-points because they are a property of the structure; the k-point side supplies only the operation index.buildAtomPermadditionally enforces the closure that every type-level sum depends on -- it accepts an image atom only within the source atom's type and stops fatally otherwise -- so that closure is verified at startup rather than assumed (DESIGN 2.3). -
O_Optc / O_OptcSpectra (
optc.F90,optcSpectra.F90): Compute and write the optical properties, both the total spectra and the partial (POPTC) decomposition of the momentum matrix element between a pair of groups.computePOPTCPairsbuilds the decomposition index that assigns every basis function to a partial, accumulates the pair matrix per transition per k-point, and applies the IBZ star average to it;printSpectrumPOPTCwrites the resulting spectra in the same partial ordering. The star average consumesatomPermfrom O_AtomicSites andfullKPToIBZOpMap/fullKPToIBZKPMapfrom O_KPoints -- the same pair of inputsbond.F90uses, applied to both indices of the pair matrix instead of one. The offered decompositions and their numbering are DESIGN 11; the unfolding is DESIGN 2.5 and PSEUDOCODE 7a.Brillouin-zone integration is presently Gaussian unconditionally:
getOptcCondandgetOptcCondPOPTCin O_OptcSpectra smear each transition and weight it bykPointWeight, and this is the one consumer that does not branch onkPointIntgCode. It will gain a LAT pathway alongside, as a separate accumulation routine rather than an internal branch, because the LAT loop inverts the nesting (DESIGN 12.2).Taking it on restructures this pair of modules to match the DOS path rather than imitate it.
optcCondandoptcCondPOPTCmove to module scope in O_OptcSpectra, the house pattern that O_OptcTransitions already follows withtransitionProb;printOptcResultssplits intocomputeOptcSpectra, which sets up and selects the pathway, andprintOptcSpectra, which writes the files; andsubroutine optcinimago.F90calls the two in turn, exactly wheresubroutine dosmakes the same choice. The accumulators are renamed toaccumulateOptcCondandaccumulateOptcCondPOPTCwith_LATcounterparts, since they accumulate into an array handed to them rather than getting anything. That pathway adds two dependencies, and they land on O_OptcSpectra rather than on O_Optc, because the accumulation routines live there: the tetrahedra andtetraVolfrom O_KPoints, andbloechlCornerDOSWtfrom O_MathSubs, which it feeds the four corner values of the band-energy DIFFERENCE rather than of a band.The two pathways share the momentum matrix elements and the decomposition index, and diverge in how the resulting transition strengths are indexed, filtered, ranged and occupied. The shared part -- the
conjWaveMomSumconstruction, which is the expensive one -- is extracted into a routine both producers call, so that neither duplicates it. DESIGN 12. -
O_DOS (
dos.F90): Two DOS paths:computeIterationTDOS(in-SCF convergence monitoring) andcomputeDOS(full TDOS/PDOS post-processing). Will gain a LAT branch incomputeDOSdispatched onkPointIntgCode. -
bond.F90: Computes bond order and effective charge (Q*). Already uses
electronPopulation_LATwhenkPointIntgCode == 1. Phase F additions: the accumulation loop will buffer per-IBZ-kpoint projections and distribute them across the star of each IBZ k-point usingatomPermandfullKPToIBZOpMapfrom O_KPoints, giving correct per-atom Q* and per-pair bond order under IBZ reduction. -
populate.F90 (
O_Populate): ComputeselectronPopulation(Gaussian/Fermi-filling path) which folds inkPointWeight.computeElectronPopulation_LAT(already implemented) produceselectronPopulation_LATusing tetrahedra data from O_KPoints. Both occupation arrays serve the same downstream consumers (bond order, effective charge); the integration method (kPointIntgCode) selects which one is used. Upstream of bond order and effective charge. -
makeKPoints (
makekpoints.F90): Legacy standalone program that generates explicit k-point lists with IBZ reduction (foldMesh). Separate executable, not linked into imago. Being retired -- its mesh-building and IBZ reduction capabilities have been ported into imago'sinitializeKPointMesh. Retained only for backward compatibility; new workflows should not depend on it. -
makeinput.py (
src/scripts/makeinput.py): Top-level orchestrator that prepares all input files for Imago. For k-points it supports three pipelines:- Mesh mode (
-kp,-scfkp,-pscfkp): writes a style-code-1 k-point file with axial counts and shift. Imago builds the full mesh internally and reduces to the IBZ, giving it full access to the mesh topology and symmetry maps needed for correct decomposition properties. - Density mode (
-kpd,-scfkpd,-pscfkpd): writes a style-code-2 k-point file with density and shift. Imago computes axial counts from the density and reciprocal cell geometry, then builds and reduces the mesh identically to mesh mode. - Explicit mode (style code 0 in the k-point file): a pre-built list of k-points with weights, read directly by Imago with no internal mesh construction. Supported for special cases (e.g., hand-crafted k-point sets), but Imago emits a prominent warning that decomposition properties (effective charge, bond order, partial DOS, partial optical properties) will not be correct unless the user has taken extreme care to provide a symmetric mesh. Not produced by makeinput.
- Mesh mode (
-
makegroups.py (
src/scripts/makegroups.py): the bispectrum type-grouping helper (8.9, DESIGN 5.10). A Fortran-side descriptor can only come from a completed Imago run, so grouping atoms by bispectrum is a sequence -- makeinput (no grouping) ->imago -loen -scf no-> read the enrichedfort.21-> bucket withBispecMatcher-> rewrite the skeleton with explicit per-element species tags -- orchestrated outsidemakeinput.py. This script owns that sequence and is dual-mode: an importablegroup_by_bispectrum(...)the producer (build_initial_potentials.py) calls, plus a__main__CLI for manual use.makeinput.pystays a plain input-writer that reads whatever explicit types the skeleton already carries; the bispectrum reasoning never enters it. Reduce grouping remains inmakeinput.pyfor now (the two environment schemes are deliberately asymmetric); the long-term direction is formakegroups.pyto absorb all type assignment (reduce, target, block), leavingmakeinput.pypurely an input-writer. -
symmetry.py (
src/scripts/symmetry.py): the first focused module split out ofstructure_control.pyat the symmetry seam (Section 7). It owns reading a cell's conventional-cell-abc space-group point operations (and fractional translations) fromshare/spaceDB/<sg>-- the singleread_conv_abc_point_ops(...)that bothmakeinput.py(the k-point-file writer, forbuildAtomPerm's fractional translations) andbuild_initial_potentials.py(the initial-potential producer, for the mesh climb's axis classes, DESIGN 3.12) read through, so one parser of the spaceDB operation file serves both and cannot drift.StructureControlexposes a thin accessor delegating here. Future symmetry concerns (supercell construction, atom permutation tables) migrate into this module as Section 7's split proceeds. -
The SCF HDF5 module family (
hdf5SCF.F90and the six section fileshdf5SCFIntg.F90,hdf5SCFElec.f90,hdf5SCFExco.f90,hdf5SCFEigVal.f90,hdf5SCFEigVec.F90,hdf5SCFPotRho.f90; the PSCF family mirrors it): the layout follows the HDF5 file itself.O_SCFHDF5is the sole owner of the FILE -- its id, property list, k-point group and the shared attribute dataspace -- and every open or close of the file in the program happens in this one module; the section modules own the groups, datasets and attributes of their subtree and receive shared context DOWNWARD as arguments, never reaching back up (theusearrows run fromO_SCFHDF5into the sections only). The parallelization work (Section 6.5, DESIGN 9.5) extends the pattern rather than bending it: the term-stage file lifecycle (suspend/resume, the shared file-name builder, and the lock-disciplined per-term write) lives inhdf5SCF.F90with the other file-level operations, while each section exports its object NAMES as public parameters beside the code that creates them, so a by-name reopen cannot drift from the creation spelling. -
O_HDF5Utils (
hdf5Utils.f90, new): the small shared layer beneath both file owners, holding only operations about HDF5 files in general -- the lock-retry open (the file lock is the write mutex of DESIGN 9.5) and the close-every-remaining-object sweep that makes "the file is truly closed" a guarantee.O_SCFHDF5uses it now; the PSCF owner uses it when post-SCF parallelization arrives, which is also the natural moment to unify the two owners' duplicated file-name construction here.
makeinput.py (top-level orchestrator)
+-- symmetry.py (read_conv_abc_point_ops: the single
| spaceDB operation reader; kp writer uses the
| fractional translations for buildAtomPerm)
+-- (mesh mode: -kp)
| writes kp-scf.dat / kp-pscf.dat (style 1)
+-- (density mode: -kpd)
| writes kp-scf.dat / kp-pscf.dat (style 2)
+-- makeKPoints (legacy, no longer called by makeinput)
build_initial_potentials.py (initial-potential producer)
+-- symmetry.py (same read_conv_abc_point_ops; the mesh
| climb uses the rotations for axis classes, DESIGN 3.12)
+-- kaleidoscope (dispatch), guidance_db (dataspace)
imago.F90 (top-level dispatcher)
+-- O_KPoints (kpoints.f90)
| +-- O_Lattice (recipCellVolume, invRealVectors)
| +-- O_AtomicSites (atom positions for buildAtomPerm)
+-- O_DOS (dos.F90)
| +-- O_KPoints (eigenvalues, tetrahedra, weights)
| +-- O_PSCFIntg / O_SCFIntg (eigenvectors via HDF5)
+-- populate.F90 (O_Populate)
| +-- O_KPoints (kPointWeight, tetrahedra,
| | eigenvalues, fullKPToIBZKPMap for LAT)
+-- bond.F90
| +-- O_Populate (electronPopulation or
| | electronPopulation_LAT)
| +-- O_KPoints (fullKPToIBZKPMap,
| | fullKPToIBZOpMap, numFullMeshKP:
| | which star each IBZ point owns)
| +-- O_AtomicSites (atomPerm, invAtomPerm:
| | where each atom lands under an operation)
| +-- O_PSCFIntg (eigenvectors via HDF5)
+-- O_Optc (optc.F90)
| +-- O_KPoints (fullKPToIBZKPMap,
| | fullKPToIBZOpMap, numFullMeshKP,
| | kPointWeight, kPointIntgCode)
| +-- O_AtomicSites (atomPerm; also valeDim and
| | atomSites for the decomposition index)
| +-- O_AtomicTypes (radial function counts per
| | type, for the QN_nl resolved partials)
| +-- O_SecularEquation (eigenvectors, momentum
| | matrix elements, eigenvalues)
| +-- O_OptcSpectra (optcSpectra.F90: broadening and
| spectrum output, total and per partial
| pair; carries both accumulation pathways)
| +-- O_KPoints (tetrahedra, tetraVol,
| | numFullMeshKP, fullKPToIBZKPMap,
| | fullKPToIBZOpMap: the LAT
| | accumulation walks tetrahedra over
| | the full mesh)
| +-- O_MathSubs (bloechlCornerDOSWt, fed
| the corner values of the band
| energy DIFFERENCE: LAT only)
- Fortran 90; CMake out-of-source build in
build/ - Compiler set via
$FC(must match HDF5-compiled compiler) - Install prefix via
$IMAGO_DIR - HDF5 required for primary I/O in imago
- Supported compilers: gfortran, ifort
The top-level CMakeLists.txt carries opt-in instrumentation
options -- IMAGO_CHECKS, IMAGO_FPE_TRAP, IMAGO_INIT_SNAN,
IMAGO_WARN_EXTRA, and IMAGO_SANITIZE -- all default OFF. With
every option off the compile and link flags are byte-for-byte
identical to the historical build, so the production default is
never disturbed; the instrumentation is only ever appended.
CMakePresets.json bundles the common combinations as named
presets (gfortran-serial / -mpi / -debug / -audit / -asan
/ -profile / -gprof), each
building into its own tree and installing to a per-preset sandbox,
so a preset can never overwrite the production install.
A build flavor is an alternative build of the engine kept under
$IMAGO_DIR/envs/<flavor>/bin -- e.g. asan for leak hunting,
audit for checks + traps + warnings, and later profiling and
parallel-debug builds. Because only the compiled Fortran
executables differ between flavors (the Python helper scripts and
the share/ database are compiler-agnostic and identical), a
flavor bin is just symlinks to the production bin with the engine
executables overlaid. The envs.sh switcher (tracked in .imago/
from the clone and sourced into the shell, alongside imagorc,
never installed as a program) provides imago_env <flavor>, which
activates a flavor by repointing only IMAGO_BIN (and
PATH / PYTHONPATH) at it -- share/ and .imago/ are reused,
so no data is duplicated and the production install is untouched.
This is the durable replacement for ad-hoc overlay directories.
Activation is per-shell, so different terminals can run different
flavors at once without interference. Operational details live in
BUILD.md; the bug-squashing campaign that motivates the
instrumented flavors is in dev/DEBUG.md.
This section records a decision, not its execution -- none of the
changes below are in the code yet. The intended default build is
the parallel one: MPI with the ELPA eigensolver (the gfortran-mpi
preset -- FC=h5pfc, IMAGO_MPI and IMAGO_ELPA on), rather than
the serial build. This follows VISION Goal 7: intra-problem
parallelism is the path a large secular problem needs, and the
teaching and research environments Imago targets (a provisioned
Ubuntu VM or a teaching cluster) carry the parallel stack, so the
parallel toolchain is the common case. A serial build is kept as an
explicit, named fallback for a host that lacks the parallel
libraries; the nanoHUB path is deferred until its own
parallelization is investigated and does not constrain this.
Naming follows from that. The serial preset is renamed from the
vague gfortran-release to gfortran-serial -- "release" named an
optimization level, not the axis that now matters, which is serial
versus parallel. Both gfortran-serial and gfortran-mpi are
optimized builds and differ only in whether the MPI/ELPA paths are
compiled in. The production install ($IMAGO_DIR/bin) is then built
with the parallel settings -- the gfortran-mpi configuration
installed to $IMAGO_DIR rather than to a per-preset sandbox -- and
gfortran-serial is reached as a flavor (imago_env serial).
Two consequences a reader must not miss. First, the default
eigensolver changes. The serial build calls LAPACK's ZHEGV, while
the MPI build uses ELPA even at a single rank, so making MPI the
default changes the default solver from ZHEGV to ELPA. These are
different algorithms, so the acceptance gate for the change is
physical agreement within the convergence noise floor across the
reference cells -- not bit-for-bit identity. (The electrostatic
setup is already known bit-identical from the performance campaign;
it is the solve that differs.) A one-rank MPI run needs no mpirun
and reproduces the serial user experience -- the launcher appears
only when a job asks for ranks -- but "one rank" is one ELPA rank,
not the serial ZHEGV path.
Second, the switcher's reserved names were decoupled from the build
name ahead of the flip (done). envs.sh reaches the installed
default only through no-argument / prod / production -- whatever
build is installed -- and every other name, serial and mpi
included, is an ordinary flavor. No build name is conflated with the
default, so the switcher stays correct whichever build production
becomes. The retired release alias had no callers.
Retiring gfortran-debug. Its flags come entirely from
CMAKE_BUILD_TYPE=Debug (-Og -g -fcheck=all -fimplicit-none -Wall -fbacktrace -DDEBUG), and both gfortran-audit and
gfortran-asan build with that same Debug type and then add
instrumentation on top -- so every check gfortran-debug provides
is already contained in them. Its checks are fully superseded. The
only property it uniquely offers is mildness: audit adds
floating-point traps, signaling-NaN initialization, and an extra
warning set that can halt or add noise on benign conditions, and
asan is heavyweight. The direction is to retire the generic
gfortran-debug in favor of the explicit builds -- gfortran-audit
for checks and traps, gfortran-asan for memory errors -- with
audit as the everyday interactive-debug build. The one thing a
clean deletion loses is a plain gdb build with no traps; a reader
who wants that can instead keep a minimal build under a name that
says so.
Implementation order, when this is taken up: (1) rename the serial
preset and build the production install from the gfortran-mpi
configuration; (2) prove the noise-floor agreement gate above;
(3) update the envs.sh reserved-name mapping (done -- the switcher
no longer conflates a build name with the installed default);
(4) update BUILD.md and the preset list in 4.1; (5) retire or
rename gfortran-debug.
These elements exist in the code and are leveraged by the current design:
kPointStyleCode: integer controlling how k-points are specified (kpoints.f90)- 0 = explicit list (legacy; not produced by makeinput; Imago emits a warning that decomposition properties may be incorrect)
- 1 = axial counts + shift (mesh built internally;
primary mode for
-kp/-scfkp/-pscfkp) - 2 = minimum density + shift (mesh built internally
via
computeAxialKPoints; used by-kpdetc.)
kPointIntgCode: integer read from input (0=Gaussian, 1=LAT) (kpoints.f90)readKPoints: reads the k-point input file; branches onkPointStyleCode. For style codes 1 and 2, reads point group operations intoconvAbcPointOpsandconvAbcFracTrans(conventional-cell-abc fractional form, written straight fromshare/spaceDB/<sg>), plus theCONV_LATTICEblock (the conventional cell in Bohr) and theCELL_MODEflag (fullorprim). See DESIGN 2.7 for the basis-invariance design. (kpoints.f90)computeAxialKPoints: convertsminKPointDensityintonumAxialKPoints(:)usingrecipMagandrecipCellVolume(kpoints.f90)initializeKPointMesh: builds the uniform mesh fromnumAxialKPointsandkPointShift; whenapplySymmetry=1, folds the mesh to the IBZ usingabcRecipPointOpsand savesfullKPToIBZKPMap(kpoints.f90)computeRecipPointOps: conjugates the conv-abc point group operations inconvAbcPointOpsinto the basis of whichever reciprocal lattice O_Lattice currently holds, using the composed change-of-basis matrixC = invRealVectors^T * M_conv(withM_convread from the kp file'sCONV_LATTICEblock). Infullmode theCELL_MODEflag triggers an identity shortcut so the loop becomes a copy. Output isabcRecipPointOps(kpoints.f90). See DESIGN 2.7.computeRealPointOps: real-space sibling ofcomputeRecipPointOps. ConjugatesconvAbcPointOpsandconvAbcFracTransinto the basis of the loaded real lattice using the sameC, producingabcRealPointOpsandabcRealFracTransforbuildAtomPerm(kpoints.f90, DESIGN 2.7). Same full-mode identity shortcut as the reciprocal-space routine.convAbcPointOps,convAbcFracTrans: on-disk conventional-cell-abc fractional form of the symmetry operations read from the kp file (kpoints.f90). Values matchshare/spaceDB/<sg>entry-for-entry; the consumer-side conjugation rebases them into the loaded cell at runtime.convLattice,cellMode: conventional-cell matrix (Bohr) andfull/primflag read from the kp file'sCONV_LATTICEandCELL_MODEblocks. Together they give the consumer the information needed to form the change-of-basis matrixCand to decide between the identity shortcut and full conjugation (kpoints.f90, DESIGN 2.7).abcRealPointOps,abcRealFracTrans: symmetry operations expressed in the basis of the real lattice currently in O_Lattice (full conventional cell or primitive reduction). Consumed bybuildAtomPerm(kpoints.f90 / atomicSites.f90, DESIGN 2.7)xyzRealPointOps: the same operations in CARTESIAN form,L R_abc L^-1with the real lattice vectors as columns. Needed because the momentum operator is a vector whose components mix under an operation, so unfolding an optical matrix element to a star member requires rotating it and not only permuting atoms. Built besideabcRealPointOpson the same branches, and therefore absent on the SYBD path and for k-point style code 0 (kpoints.f90, DESIGN 13.4)sc.full_cell_real_lattice(structure_control.py): snapshot of the conventional lattice captured beforeapply_space_group()may overwritesc.real_latticewith a primitive reduction; sibling offull_cell_mag/full_cell_angle; same 1-indexed layout asreal_lattice; written into the kp file'sCONV_LATTICEblock by the kp-file writer inmakeinput.py(DESIGN 2.7).- Note: an earlier iteration of this design used a
producer-side
_to_cartesian_opshelper to convert spaceDB operations to a Cartesian xyz intermediate on disk. Under the current design the producer is a near-passthrough: spaceDB operations are written unchanged; the only added responsibilities are emittingCONV_LATTICEfromsc.full_cell_real_- latticeandCELL_MODEfrom the skeleton flag. See DESIGN 2.7 for the diagnostic history. generateTetrahedra: tiles the full uniform mesh with tetrahedra, six for each of the box's long diagonals that are used, so 24 per sub-cube by default. Cutting along all four is what makes the crystal point group carry the decomposition onto itself;NUM_TETRA_DIAGONALSin the k-point file selects 1 instead for runs that need only totals (DESIGN 1.2). UsesgetIndexFromIndicesfor periodic wrapping. Called frominitializeKPointswhenkPointIntgCode == 1(kpoints.f90)- Atom-resolved tetrahedron results are averaged over
the point group before being written, using the
permutation tables the unfolding already builds.
This reaches the lattices the decomposition alone
cannot (hexagonal, rhombohedral) and is on by
default, controlled by
SYMMETRIZE_LAT_PARTIALSin the k-point file alongside the two settings above, since averaging the result is equivalent to averaging the integration weights over the star (DESIGN 1.7) getIndexFromIndices(a,b,c): converts mesh indices to the linear k-point index (kpoints.f90)energyEigenValues(n, i, spin): fully in memory after SCFrecipCellVolume,invRealVectors: available from O_Lattice afterinitializeLatticeatomSites(:)%cartPos,atomSites(:)%atomTypeAssn: atom positions and type assignments from O_AtomicSites, available after input parsing
Measured on imago 2026-08-18/19 (
dev/PERFORMANCE.md, "Baseline" and "Coarse time map"; the record of the numbers lives there, this section carries only what they decide). Two operations account for 72-86 % of every run above toy size: the ONE-TIME three-centre electronic-potential integrals and the PER-ITERATION secular solve, the solve's share growing with system size (9 % -> 62 % along the complex column at 8 -> 1296 atoms, the endpoint at ONE k-point; 27-37 % at 243-1296 atoms real). The real-space grid work that an earlier draft of 6.5 called "the cheap, proven win and the right first target" is 1-3 % of the medium cells. Section 6.8 holds the table; 6.5 and 6.6 are written from it. The programmer's ruling on that evidence (2026-08-18): go straight for distributed-memory MPI with the ELPA eigensolver -- the dominant use case is large cells at one k-point -- with in-node threading as a companion inside a rank, not a separate stage. The intra-node directions 6.1-6.4 (precision, SIMD, GPU) stand but follow it.
The program will support compilation in both double
precision (64-bit) and single precision (32-bit) via a
compile-time kind parameter (e.g., wp for "working
precision") defined in kinds.f90. All floating-point
declarations, literal constants, and MPI/HDF5 type tags
must use this parameter so that a single preprocessor
flag or CMake option switches the entire build.
Motivation. Single precision doubles SIMD throughput on CPU (8 floats vs. 4 doubles in AVX-256) and is essential for GPU performance, where double-precision throughput is 2x slower on data-center GPUs and up to 32x slower on consumer hardware.
Numerical stability. Single precision provides only
~7 significant digits, which raises concerns in
accumulation loops where small contributions are added
to large sums. The lattice loop in the integral engine
(gaussOverlapOL and siblings in integrals.F90) sums
contributions from nearest to farthest replicated cells;
far-away terms are individually small but collectively
significant. At double precision this is safe (~15
digits of headroom). At single precision, compensated
(Kahan) summation or smallest-first accumulation order
may be required for specific accumulation sites.
Stability analysis must accompany the precision switch,
targeting at minimum the lattice-sum accumulators and
the eigenvalue solver interface.
The alpha-pair iteration inside the lattice loop of
gaussOverlapOL (and the analogous nuclear and
three-center routines) is currently control-flow-
dominated: a do while (.true.) loop with mode
switching, early exits, and per-pair negligibility
tests. This structure prevents SIMD vectorization.
Target restructure. Separate the selection phase
(which alpha pairs survive the alphaDist threshold for
a given atom pair and lattice distance) from the compute
phase (evaluate the Gaussian overlap integral). The
selection phase produces a packed list of surviving alpha
pair indices; the compute phase processes that list in a
tight, branchless loop amenable to SIMD.
This "gather, filter, compute in bulk" pattern is the same restructure needed for GPU offload, so the two goals reinforce each other.
Long-term, the integral engine and eigenvalue solver are candidates for GPU acceleration. The restructured inner loops from 6.2 translate almost directly to GPU kernels: the packed alpha-pair list becomes a kernel launch over surviving pairs. Single precision from 6.1 is required to achieve full GPU throughput.
Practical path:
- Introduce
wpkind parameter (6.1) and validate numerical stability at single precision. - Restructure alpha-pair loops into gather/compute phases (6.2); verify SIMD vectorization on CPU.
- Offload the compute phase to GPU (OpenACC, CUDA Fortran, or OpenMP target) as a later step.
Each phase is independently useful: phase 1 halves memory footprint and improves cache behavior; phase 2 speeds up CPU execution; phase 3 adds GPU capability.
An earlier vectorized implementation exists outside this repository at:
/home/rulisp/lewis/CPG/cpg-repo/v34/src/olcao/
integralsSCF.vec.F90 (SCF integrals)
integralsPSCF.vec.F90 (post-SCF integrals)
gaussIntegrals.vec.f90 (vectorized primitives)
The vectorized gaussIntegrals.vec.f90 primitives are
also tracked in this repo under src/imago/
src/imago/.
These files demonstrate the gather/compute separation
described in 6.2. The alpha-pair selection loop
(branchy, scalar) collects surviving pairs into
lmAlphaPairs organized by angular momentum type. A
packing step linearizes them into orderedAlphaPairs
with segIndices marking contiguous same-type regions.
A single batched call (overlap2CIntgVec) then
processes all pairs, enabling SIMD over the contiguous
segments. This pattern is the starting point for the
restructure in A4.
The directions above (6.1-6.4) attack intra-node performance. This section addresses intra-problem scaling across nodes: distributing one large secular problem over many MPI ranks so a system too large for a single node's memory still runs. This is the first of the two parallel axes named in VISION Goal 7; the second, inter-problem throughput, is the flight layer of section 9.
A sibling branch of the common-ancestor OLCAO code (see 6.7) located where MPI must reach into the algorithm; the 2026-08-18 measurement (6.8) says which of those seams carry the time. In that order:
- The three-centre electronic-potential integrals (the
Electronic Potential Integralsstage,elecPotGaussOverlapinintegrals.F90) -- 37-45 % of every medium run and the largest single stage on every deck. Computed ONCE per run. Its structure is a loop over the potential terms (potDimof them: every alpha of every potential type; 32 for a two-type glass, more for crystals with more types), and for each term one full valence-by-valence matrix is assembled over all atom pairs and lattice cells (gaussOverlapEP), orthogonalized to the core (ortho) and written to its OWN HDF5 dataset carrying a completion attribute. Terms are independent of one another, and within a term the atom pairs are independent of one another, so the stage offers TWO orthogonal decompositions, and both have been tried in the ancestor codes (6.7):- By potential term. Each rank takes a subset of the terms,
computes complete matrices for them, orthogonalizes and writes
its own datasets: no reduction, no exchange of partial
matrices, and the per-term completion attribute already makes
the stage restartable term by term. Width
potDim. Its weakness: terms are NOT equal in cost -- a diffuse alpha (small exponent) reaches many more (atom pair, lattice cell) entries than a tight one (theanyElecPotInteractionmask that prunes pairs per term is exactly the record of that inequality) -- so an equal COUNT of terms per rank is not an equal share of work. - By atom pair within a term. Each rank takes a range of the
packed pair index (
numAtomSites*(numAtomSites+1)/2pairs) and accumulates its pairs' contributions into a DISTRIBUTED matrix (the polcao effort of 6.7 did this with Global Arrays; the upolcao design of DESIGN 9.3-9.4 does it with a block-cyclic layout and a few redundantly computed pairs). Width is the pair count, so it scales with the cell, and it produces the matrix already distributed for the solve. Its weakness is the mirror image: pairs are not equal in cost either (basis size, neighbours within the cutoff), and the accumulation is a communication pattern rather than a private write. The choice between the two (or their two-dimensional product: terms across one axis of the process grid, pair blocks across the other) was a DESIGN 9 decision resting on one measurement, the per-term and per-pair cost distributions on the medium decks (PERFORMANCE.md "PA1 cost distributions"). DECIDED 2026-08-18, BY TERM (DESIGN 9.5): per-term costs are mild and monotone in the alpha exponent (max/mean 1.5-1.6), and on multi-k-point cells most of the stage is the per-term core-orthogonalization and dataset write OUTSIDE the pair loop -- work a term distribution carries to the ranks for free and a pair distribution would leave serial. The load balance is decided with it: a STATIC deal of terms, most-diffuse (costliest) first; a dynamic work queue is unnecessary at the measured imbalance, and the pair-level split under 9.4 (with the two-dimensional product it enables) is the later refinement for rank counts abovepotDim. Output through serial HDF5 works for the first version (ranks write disjoint datasets in turn); collective parallel HDF5 (DESIGN 9.7) is a later optimization, not a prerequisite.
- By potential term. Each rank takes a subset of the terms,
computes complete matrices for them, orthogonalizes and writes
its own datasets: no reduction, no exchange of partial
matrices, and the per-term completion attribute already makes
the stage restartable term by term. Width
- The secular solve -- the generalized eigenproblem
H c = e S c, per iteration, per k-point, per spin. 27-46 % of the medium runs and the share that GROWS with size, as N^3 must overtake the integrals' pair count. This is the scaling wall and the subject of 6.6: ELPA over a block-cyclic layout (DESIGN 9.3), with k-point and spin as an outer, communication-free distribution when there are several of either. Note that at the sizes measured (valence dimension 12528 for the 1296-atom glass: 1.26 GB real, 2.5 GB complex -- corrected 2026-08-23 from the run logs, which report valeDim directly; the earlier figure of 5184 was that deck's numSTATES, the first of the three NUM_STATES_TO_USE values, not its matrix dimension) the matrix still fits one node many times over -- the driver for distributing the solve is TIME, not memory, until the dimension reaches tens of thousands. - Valence charge density (
Valence Charge Density, per iteration) -- the third cost, 5-14 %; it consumes the eigenvectors the solve produces and parallelizes over the same distribution. - Real-space grid work -- the site and grid loops in
elecStatandexchCorr(Make SCF Potential,Electrostatic Matrices). Embarrassingly parallel by the one-dimensional load balance plusMPI_REDUCEof DESIGN 9.2 -- and, on imago, 1-3 % of a medium run (the electrostatic setup is 14 % of the 243-atom glass, once). It is the LAST target, not the first; the earlier draft of this list had it first on inherited authority, and the measurement reversed that. - The atom-pair loop family (added 2026-08-25; DESIGN 9.9) --
the two-centre integral stages (overlap, kinetic, mass velocity,
nuclear potential) of the SCF setup and EVERY integral stage of
the post-SCF setup (
intgPSCF: the overlap again, the full HamiltoniangaussOverlapHamPSCF, and the dipole, momentum and k-overlap property integrals). All share one loop over atom pairs accumulating into full in-memory matrices that root then orthogonalizes and writes. A few percent of an SCF run at the medium sizes, N-squared in growth; but on the post-SCF path the Hamiltonian stage carries the ENTIRE three-centre cost with no term axis, so a band-structure or optical run's setup is root-serial in full today. The deal is the 9.2 pattern on the packed pair index with a sum-onto-root reduce, taken in two steps (programmer's ruling, 2026-08-25): the SCF-path stages first, with no change to the run shape; then the post-SCF path, which changes the RUN SHAPE -- every rank entersintgPSCFfor the integral stages and then servesbandPSCF's solves, which the k-point deal of DESIGN 9.6 already distributes. The sum-onto-root form holds a full copy of every matrix on every rank and is accepted for TODAY's cells only: the stated target of 10,000 atoms and more (6.9) makes a full copy impossible per rank, so the block-cyclic assembly of DESIGN 9.4 -- nothing assembled whole -- is the destination, and the first form is built so that its pair partition carries over into it. - The post-SCF analysis family (added 2026-08-26; DESIGN
9.10) -- the density of states, bond order, optical, dipole and
field programs that run after the eigenpairs exist. One shape
for all of them: per spin and k-point, read that k-point's
eigenvectors and integral matrices from root's file, contract
them state by state, add into accumulators that are plain sums
over k-points and states. Root-serial in full today, on both
entry points. Designed on the PDOS first (programmer's choice):
the projection is a matrix product in disguise (
T = S C, one GEMM per k-point in place of a strided intrinsic sum per state-and-basis-function pair -- the largest gain, single-core, taken first and measured on its own), then a deal chosen by DESIGN 9.6's rule -- whole k-points to ranks when k-points are plentiful, contiguous blocks of STATES within every k-point when they are not, which is what gives a one-k-point cell of 10,000 atoms an axis with width -- and a sum-onto-root combine. Run shape ruled 2026-08-26: no server; root shuts the solve servers right after the last solve and every rank calls the analysis routines in lockstep, root alone reading and writing the file. The state deal holds the full overlap on every rank, the same stepping stone as the pair family with the same expiry: at the target sizes the product becomes a PBLAS operation on the block-cyclic layout with the eigenvectors left where ELPA put them.
Two patterns to name so they are not mistaken for progress: the
replicate-and-broadcast form (rank 0 computes, MPI_BCAST makes
every rank consistent) is correct and gives NO speedup; and a
serial LAPACK solve on a broadcast copy of the full matrix on every
rank is the same thing for the eigenproblem. Both are what the
sibling branch stopped at (6.7).
Build and runtime shape. The MPI lifecycle, the one-dimensional
load balancer and the process-grid helper belong in one module
(O_MPI -- mpi.f90 in the sibling branch, mpi.F90 here, the
capital F carrying the preprocessor guard; use mpi_f08); imago
does not yet carry it, and introducing it is the first code-level
step. The toolchain exists and is verified (BUILD.md, "Parallel
toolchain"): the cpgp conda environment (dev/env/cpgp.yml)
supplies openmpi, MPI-enabled HDF5 (h5pfc), ScaLAPACK and ELPA on
the same gfortran/OpenBLAS family as the serial build; the
gfortran-mpi preset (IMAGO_MPI, IMAGO_ELPA) builds the
unchanged serial source against it and reproduces the serial
baselines, so the parallel work starts from a binary that is
already known to be right at one rank. A parallel run is launched
by imago.py through the existing IMAGO_EXE hook (mpirun -np
inside a batch script is the group's precedent).
Ranks only: one thread per rank (programmer's ruling,
2026-08-28). The parallelism model is pure MPI. Threads are
NOT a parallel axis: every rank runs one thread
(OMP_NUM_THREADS=1, OPENBLAS_NUM_THREADS=1 by policy), and
every phase that must go faster is distributed over ranks,
never over threads within a rank. This reverses the
companion-threads position held here since 2026-08-18, and the
measurements are why. A monolithic MPI job holds a FIXED rank
count for its whole life -- the world size is set at mpirun
launch and the allocation cannot be reshaped mid-run -- so it
cannot trade ranks for threads between phases. Yet the phases
pull opposite ways under one fixed split: the setup integrals,
the electrostatic setup and the eigensolve all want ranks (the
solve is thread-proof -- 16 OpenBLAS threads bought ZHEGV four
percent, PERFORMANCE.md "PA4-A"), while the one phase that
wanted threads was the valence-density build, whose rank-k
recast (DESIGN 9.6, PSEUDOCODE 31) speeds up under threaded
BLAS. Serving both at once is impossible with a fixed split,
so the ruling removes the conflict at its root: make the
valence build want ranks too, and the whole run collapses to a
single resource number. The payoffs are a one-axis resource
request (DESIGN 14, the resource specification), scaling that
crosses nodes where threads cannot, and one parallelism story
to document and hand on rather than a per-problem
ranks-by-threads tuning.
Two consequences follow, each its own work item. First, the valence-density build must be DISTRIBUTED over ranks rather than left to threaded BLAS: the rank-k recast (PSEUDOCODE 31) stays as the per-rank kernel, and its state-axis deal -- optional under the retired companion-threads position -- becomes REQUIRED (DESIGN 9.6, amended). Second, thread- oriented code this ruling strands -- the OpenMP-over-atom-pairs companion named just above, any thread-count knob, any unused OpenMP directive -- is to be surveyed and REMOVED, so a dead option cannot sit in the source and confuse a later reader. The build still links a threaded BLAS; this ruling is a runtime policy and a design commitment, not a rebuild.
The secular solve must sit behind one interface with swappable backends, chosen by problem size and available hardware:
- Serial LAPACK (
ZHEGV/DSYGV) -- today's path, correct for systems that fit one node, and the one-rank backend of the parallel binary (with OpenBLAS threads inside the rank). - Distributed: ELPA -- DECIDED 2026-08-18 (programmer's ruling;
the toolchain carries it,
elpa_init/elpa_allocateverified on the cluster). ELPA consumes the block-cyclic BLACS layout of DESIGN 9.3, solves the GENERALIZED problem by a Cholesky transformation followed by its two-stage tridiagonalization, and outperforms ScaLAPACK'sPZHEGVXon the dense generalized problem typical of an LCAO method; one API spans CPU and GPU. ScaLAPACK stays linked beneath it (ELPA requires it) and remains available for the layout utilities (numroc, descriptors) and as a fallback, not as the solver. - GPU -- via ELPA's GPU backend, or a vendor solver (cuSOLVER / MAGMA) for the single-node case; later.
This is the single most important architectural boundary for
parallelization. The sibling branch declared a ScaLAPACK
generalized-eigensolver interface (PZHEGVX) but never called it;
imago defines the boundary cleanly -- one call site that takes H
and S in whatever layout the backend needs and returns eigenpairs
-- and completes the distributed arm with ELPA rather than
inheriting an unfinished stub. What the measurement adds (6.8): at
today's sizes the solve is distributed for TIME (its share is 46 %
at 135 atoms and rising), so the interface must be cheap to enter
and leave per iteration; the redistribution of H and S from
whatever layout the integral stage produced into the block-cyclic
one, and of the eigenvectors back out, is part of the solve's cost
and belongs inside the boundary where it can be measured. Because
the abstraction names CPU and GPU backends behind one call site, it
is also where VISION's "device placement is per-kernel" principle
(Principle 14) is realized for the most expensive kernel in the
program.
Two earlier parallel efforts on the common-ancestor OLCAO code live outside this repository, and they took the two decompositions of 6.5's first bullet:
~/olcao/src/polcao/ the earlier effort (J. Currie): atom-pair
distribution of the integral routines
~/olcao/src/upolcao/ the later effort: MPI lifecycle, grid
parallelism, block-cyclic design
polcao split every SCF integral routine (overlap, kinetic,
nuclear, and the electronic-potential contributions) over ranks BY
ATOM PAIR: loadBalMPI deals a contiguous range of the packed pair
index numAtomSites*(numAtomSites+1)/2 to each rank, each rank
computes only its pairs, and the contributions accumulate into
distributed matrices held in Global Arrays, with block-cyclic
descriptors for the core-valence orthogonalization; the
electronic-potential term loop stays serial with all ranks working
each term's pairs (ga_sync per term). Its parallelSubs.F90
carries the load balancer, the packed/unpacked pair-index maps and
the distributed-write helpers, and its own header says what it is:
a summer's blitz, results-first, with several implementations of
the same balancing idea. What it proves: the pair decomposition
works end to end and yields the matrix already distributed. What
it leaves open: the pair-cost imbalance. Its Global Arrays
dependency is NOT carried forward (programmer's ruling,
2026-08-18): the later upolcao work re-implemented, in plain MPI
and the block-cyclic layout, the parts of Global Arrays polcao
actually used -- the distributed matrix, the local/global index
mapping, the accumulation into it -- so the dependency had already
been retired before imago; imago's distributed matrices live in
the block-cyclic BLACS layout of DESIGN 9.3, which ScaLAPACK and
ELPA consume directly and the cpgp toolchain already provides.
polcao is mined for its decomposition and its index maps, not for
its communication layer.
upolcao is the later branch:
It is not a parent of imago and not a merge source. upolcao branched from an older serial OLCAO to add MPI, while imago rebranded and further developed that same ancestor on a separate path. The two have since diverged file by file, so upolcao is mined for design, not cherry-picked.
What is real and reusable:
mpi.f90-- MPI lifecycle (initMPI/closeMPI/stopMPI), a one-dimensional load balancer (loadBalMPI), a most-square process-grid helper, and a parallel-HDF5 attribute check (checkAttributeHDF5).- Genuine grid parallelism in
elecStatandexchCorr(theloadBalMPI+MPI_REDUCEpattern of 6.5). - ScaLAPACK interface declarations (
PZHEGVX,pzgemm,pdgemm,numroc, a BLACS descriptor type) that are tedious to write correctly and reusable as-is. - The block-cyclic distribution reasoning in
mpi.f90's header, including the deliberate decision to recompute a few atom pairs redundantly rather than communicate partial sub-matrices (transcribed in DESIGN section 9).
What is design-only and must not be mistaken for working code:
createProcessGridandcreateBlockCyclicDistributionare empty stubs.- The ScaLAPACK eigensolver interface has zero call sites; the SCF solve still runs serial LAPACK on a broadcast copy of the full matrix, on every rank.
- The integral and secular paths use replicate-and- broadcast, which is correct but does not scale.
The honest summary: upolcao finished the grid-parallel half and left the distributed-linear-algebra half as an elaborate design with stubs. imago can take the finished half directly and treat the unfinished half as a vetted design study rather than a starting implementation. On the measurement in 6.8 the finished half is also the half that matters least: what upolcao parallelized is 1-3 % of an imago run, and what it left as stubs is where the time is.
The record -- decks, environment, repeats, the per-operation
tables -- is dev/PERFORMANCE.md ("Benchmark deck set", "Baseline",
"Coarse time map"); it is produced by dev/tools/timemap.py from
the operation stamps every run already writes. This section keeps
only the shape of the result, because the shape is what the
architecture rests on. Share of the run's wall time, serial, one
core, plain SCF from the initial potential:
deck (atoms, variant) 3-centre secular valence grid
integrals solve density work
BN 8, complex, 4 k 66 % 9 % 5 % 8 %
BN 8, real, Gamma 77 % 0.6 % 0.3 % 8 %
KNbO3 40, complex, 4 k 37 % 37 % 14 % 8 %
KNbO3 135, complex, 4 k 40 % 46 % 9 % 2 %
SiO2 243, real, Gamma 45 % 27 % 5 % 17 %*
SiO2 1296, real, Gamma 35 % 37 % 3 % 22 %*
SiO2 1296, complex, 1 k 22 % 62 % 3 % 11 %*
(* dominated by the once-per-run electrostatic-matrix setup, not a repeated loop: 14.5 % at 243 atoms, and ~12050 s at 1296 atoms on BOTH variants -- real-space work, independent of k -- which is 21 % of the real run and 11 % of the complex one.) The stamps cover more than 99.8 % of every run, so the table is complete, and every repeat reproduced its wall time to under 1 % on the medium cells.
What the shape decides:
- Two costs, of two kinds. The three-centre integrals are paid ONCE per run; the solve is paid every iteration and every k-point. Parallelizing the first removes a fixed cost; parallelizing the second removes a cost that scales with the iteration count and, through N^3, with size. Both must be done; 6.5 orders them by their share and 6.6 owns the second.
- The solve overtakes the integrals as the cell grows, and at 1296 atoms it has: 9 -> 14 -> 37 -> 46 -> 62 % along the complex column (the 62 % at ONE k-point), 37 % real at Gamma. The large pair, being the same cell real and at one shifted k-point, also gives the real-vs-complex cost ratio directly: whole run 1.93x, solve 3.27x, three-centre integrals 1.17x.
- The grid work is not where the time is. The inherited claim in
an earlier draft of 6.5 was reasonable for the sibling code and
is wrong for imago; DESIGN 9.2's load balancer stays valid as a
design and moves to the end of the queue. One caveat the large
pair adds: the once-per-run electrostatic-matrix SETUP grows
with size (21 % of the 1296-atom real run) and becomes the
serial Amdahl bound once the integrals and the solve are
distributed -- it is the
elecStatwork upolcao already parallelized (6.7), so the last target is last, not never. - What A3-A6 target -- the alpha-pair inner loops of the integral engine, precision, SIMD, GPU -- is inside the largest stage, so those directions are not displaced by this measurement; they are intra-node work on the same stage that MPI distributes between nodes, and they follow it in time (banner above).
6.5-6.8 order the parallel work by where the time is TODAY, at one to two thousand atoms. This subsection answers a different question the programmer asked: what stands between imago and a ten-thousand-atom calculation held in distributed memory across many nodes. The answer is worth recording because it is not "more of 6.5" -- the arithmetic permits the calculation, and the program's structure is what does not.
The arithmetic permits it. Scaling from the measured 1296-atom glass (valence dimension 12528, core dimension 3024, 5184 states, 32 potential terms), a 10,000-atom cell of the same composition gives a valence dimension near 97,000, core dimension near 23,000, and about 40,000 states. The potential term count does NOT grow: it is set by the number of potential TYPES times their alphas, so it stays at 32 for a two-type glass however large the cell.
quantity at ~10,000 atoms size (real / gamma)
one full square matrix 75 GB
one packed matrix 37.6 GB
the 35 stored integral matrices 1.32 TB (~230 GB on
disk at the
measured 5.7x)
the eigenvector block 31 GB
the core-valence block 18 GB
Spread over a thousand ranks the integrals are 1.3 GB each and the solve's working set is a few hundred MB each. MEMORY IS NOT THE OBSTACLE, provided nothing is ever assembled whole. Extrapolating the MEASURED ELPA time (155 s per solve at dimension 12528 on eight cores) by the cubic law and allowing for imperfect scaling puts one diagonalization near a quarter hour on a thousand cores and a twelve-iteration run near three hours of solving. That is an extrapolation, not a measurement, and it will be optimistic -- but it describes a real calculation on an ordinary allocation, not a fantasy.
The structure does not permit it. Four facts, counted in the source rather than estimated:
- The whole-matrix assumption is pervasive. 48 sites allocate a full square valence matrix, and the full valence matrices are referenced more than 350 times across 14 files (integrals, secular equation, valence density, DOS, bond, optical, field). At dimension 97,000 each such allocation wants 75 GB in ONE process. Some would still fit on a large node, which is worse than failing: the run would proceed with one process holding everything and the rest waiting.
- The parallel increments built so far inherit that assumption. PA3 through PA5a distribute WORK while keeping root the custodian of DATA: root assembles each Hamiltonian, ships pieces out, collects results, writes. Even the ELPA arm has root unpack a complete matrix and scatter it, and the valence-density design of DESIGN 9.6 ships wave functions the same way. This is correct at today's sizes and does not extend. The hard wall is near 15,000-20,000 atoms, where one matrix exceeds any node; the practical wall -- where root's serial custody dominates the run -- arrives well before it.
- The FILE is the interface between stages, not merely storage. The integral stage writes matrices, the solve reads them back to assemble H, the valence density reads them again. No stage hands a matrix to the next in memory. An in-memory calculation therefore replaces the layer through which every stage talks to every other; it is not an optimization applied to the present flow.
- The integral stage's width is capped at the term count. PA3 divides by potential term, and that count is fixed by the type count -- 32 for a two-type glass. On a thousand ranks, 968 would idle. DESIGN 9.4's pair decomposition stops being a refinement at this scale and becomes mandatory.
Three further consequences follow. RESTART is currently free (completion attributes in the file) and would have to be rebuilt on collective checkpoint writes, which reopens DESIGN 9.7's chunk-alignment question from the write side. EVERYTHING AFTER THE SCF LOOP -- DOS, bond, optical -- is serial and reads the eigenvector block through one process. And MULTI-NODE OPERATION IS UNTESTED: we have verified that a job launches across two nodes and that ELPA's handshake completes, and nothing else. No calculation has ever run across nodes.
What is already right. The distributed eigensolver is in and working, and it is the piece that matters most. The electrostatic setup already divides on an axis that grows with the atom count. The messaging layer, the control protocol for waking ranks into a distributed region, and the file suspend/resume discipline all exist and are exercised. DESIGN 9.4 already commits to the born-distributed strategy; nothing implements it.
The judgment. Ten thousand atoms is not a continuation of the present campaign but a second one that would replace parts of the first: the k-point deal, the scatter-based ELPA path and the valence-density deal would each be reworked rather than extended, because all three rest on a custodian that the target forbids. That is not an argument against doing it; it is an argument for deciding it deliberately and sizing it honestly. The prerequisite measurement is one we have never taken -- a real multi-node run at TODAY's sizes -- since every number above assumes the interconnect and the parallel filesystem behave across nodes on the evidence of a single handshake.
Finally, a fork this subsection names rather than assumes past: at ten thousand atoms and beyond it is a genuine question whether DENSE matrices are the right representation at all. The basis is local orbitals, so the matrices are increasingly sparse with cell size, and 6.10 measures how sparse they already are.
Because the basis is LOCAL orbitals, a pair of atoms beyond some separation contributes nothing above the negligibility threshold and its block of the matrix is exactly zero. The consequence is a crossover rather than a gradient: while the cell is SMALLER than the interaction range every atom reaches every other and the matrix is dense; once the cell exceeds it, the nonzero count PER ROW stops growing, and from there the fill fraction falls as 1/N. 6.9 asks whether dense matrices are the right representation at ten thousand atoms, and this subsection measures where that crossover actually sits. Counted directly from stored matrices of two glass decks -- exact zeros, and elements above 1e-8 of each matrix's own peak.
deck / matrix fill nonzeros atoms in
per row range
243 atoms overlap 95.8 % 2,251 233
1296 atoms overlap 31.8 % 3,981 412
1296 atoms nuclear 33.2 % 4,164 431
1296 atoms kinetic 37.0 % 4,639 480
The 243-atom cell is too small to have any sparsity to exploit. Its 233 atoms in range out of 243 say that every atom sees every other. That is not a property of the material but of the cell: the cell vectors (which are in BOHR) give half-widths of 7.5, 7.5 and 8.2 Angstroms, and the interaction radius implied by the 1296-atom count -- 412 atoms at that deck's 0.0787 atoms per cubic Angstrom -- is about 10.8 Angstroms. The 1296-atom cell's smallest half-width, 11.9 Angstroms, only just exceeds it. So THE CROSSOVER IS AT ROUGHLY 900 ATOMS FOR THIS MATERIAL, and the largest cell we run sits just past it.
Projecting with the per-row count held fixed at about 4,000 gives 16 % fill near 2,600 atoms, 4 % near 10,000, and 1 % near 40,000. That projection rests on TWO measured sizes, one of which is range-limited by its own cell, so the saturated per-row count is really established at a single size; a third deck near 2,500-5,000 atoms, where the count should be unchanged, is what would confirm it. The radius is also material-dependent -- it follows from the basis functions' extent and the threshold, so a material with more diffuse valence orbitals reaches further.
The matrices are not one population. The stored potential terms span a thirty-fold range, monotone in the alpha exponent and visible both in fill and in how well each compresses on disk: the most diffuse term of a type is 30.7 % full (compressing 3.3x), the middle of the series 16.2 %, and the tightest 1.1 % (compressing 54x). The long-ranged terms are a small minority of the set.
Three architectural consequences, none of them decided here:
- Sparsity would pay first where the measured bottleneck actually is. At 4 % fill the ten-thousand-atom integral set falls from the 1.3 TB of 6.9 to under 60 GB, and the 44 GB per iteration the SCF loop currently inflates (DESIGN 9.7) falls with it. That is a storage and I/O win available WITHOUT changing any algorithm.
- It does not by itself help the solve. ELPA consumes dense block-cyclic matrices, so a sparse Hamiltonian handed to it gains nothing in arithmetic. Sparse storage that expands to dense blocks only at the solve is coherent and is probably the first step; making the SOLVE exploit sparsity is a change of method, not of data structure.
- The program already computes the sparsity pattern and
throws it away. The
anyElecPotInteraction(pair, cell) negligibility mask of DESIGN 9.5 is exactly the record of which blocks are nonzero. Storing in a sparse format would be using information imago already has.
One thing that will NOT become sparse: the Coulomb interaction is genuinely long-ranged and no threshold makes it local. imago already handles that separately, by the reciprocal-space treatment in the electrostatic setup -- which is WHY the stored potential terms can be short-ranged at all. The long-range problem is already factored out of these matrices.
Open directions, named so they are weighed rather than drifted into: numerical thresholding beyond the present mask (measured here to buy a further 30 % on the overlap at 1e-8 of peak); hierarchical low-rank treatment of the far field, which approximates rather than discards distant blocks; iterative subspace solvers that need only sparse products; and density-matrix methods that avoid eigenvectors altogether and are the only route to genuine linear scaling -- with the caveat that their nearsightedness argument requires a gap, and imago deliberately targets metals as well as insulators.
The Python side of the toolchain is currently anchored
by src/scripts/structure_control.py, a single module
holding the StructureControl class plus a large set
of related utility routines (9,335 lines, 171 tests).
It was produced by the Perl-to-Python port of
StructureControl.pm and inherits that file's
monolithic shape: everything related to an atomic
structure -- coordinate manipulation, fractional-to-
Cartesian conversion, neighbor finding, bond analysis,
PDB/skeleton I/O, symmetry operations, supercell
construction, reaction-template surgery, and more --
lives in one class and one file.
This shape has carried the codebase well enough to reach feature completeness on the port, but it is approaching the limits of what a single file should be asked to hold:
- Navigation cost. A file of 9,000+ lines is
hard to search, hard to diff meaningfully in code
review, and hard for students to study as an
example. The "Prefer concise, self-documenting
names" rule in
CLAUDE.mdaims at readability; file size is the same problem at a higher scope. - Domain bleeding. New force-field concerns
(UFF bond parameters, geometry-derived angle
clustering, LAMMPS type-ID unification) that
arise in
condense.pyandmake_reactions.pywant a natural home. Adding them tostructure_control.pywould couple its atomic-structure responsibility to force-field data prep -- a domain mismatch. Forcing everything that touches a structure intostructure_control.pyturns it into a miscellaneous bin. - Test concentration. 171 tests in one suite already means slow runs and noisy coverage reports; adding more erodes the value of running the full suite on small changes.
- Change risk. Every edit to
structure_control.pynow touches a file imported by almost every script. Small localized changes carry outsized blast radius on import-time errors and accidental namespace collisions.
Proposed direction (future work; no timeline).
Split structure_control.py along natural domain
seams, with the StructureControl class kept as
the main public entry point but its implementation
distributed across focused modules. Candidate
seams, in rough order of lowest churn first:
element_data.pyalready exists; the element-lookup code that currently lives instructure_control.pyshould migrate there.geometry.py-- coordinate math, fractional/Cartesian conversion, distance and angle kernels, minimum-image wrapping. Pure functions, heavily reused, no state.symmetry.py-- space group operations, supercell construction, atom permutation tables.neighbors.py-- neighbor search, bond analysis, coordination shells.io_skeleton.py,io_pdb.py,io_lammps.py-- format-specific I/O, each file small and focused on one external schema.structure_control.py(remaining) -- theStructureControlclass itself, delegating concrete operations to the modules above. Imports stay stable from external callers' perspective; internal implementation becomes discoverable.
Force-field concerns (UFF bond parameters, angle
clustering) live in their own separate modules
(bond_utils.py, angle_utils.py, or a shared
ff_utils.py) and never mix into
structure_control.py's domain.
Criteria for starting the split. This is a refactor, not a bug fix, so it should be scheduled deliberately rather than bolted onto feature work. Reasonable triggers:
- The file crosses 12,000 lines or the test suite crosses 250 tests.
- Two independent pieces of feature work stall on
each other due to merge conflicts inside
structure_control.py. - A new concern (e.g., machine-learning potential support, a new I/O format) wants to live near the structure code but clearly does not belong in the monolith.
Until one of those triggers hits, the conservative
choice is to keep structure_control.py as-is and
route new force-field-adjacent logic into small
dedicated modules (the path taken for
angle_utils.py in the angle-handling rework).
Starting the split at the symmetry seam (2026-07-14).
The initial-potential producer's adaptive mesh climb
(DESIGN 3.12; ARCHITECTURE 9.7) is the first concern to hit
trigger #3. To seed and step the climb the producer must
compute a cell's k-point axis classes before it dispatches
any run, which means reading the cell's conventional-cell-abc
space-group point operations from share/spaceDB/<sg> in
Python. That reader currently lives only inside
makeinput.py (the private _extract_point_ops, which the
k-point-file writer uses); duplicating it in the producer
would leave two parsers of one file free to drift, and the
climb's later self-check -- do the producer's Python axis
classes match the operations imago actually used (DESIGN 2.7)
-- is only meaningful if the climb and the kp writer start
from the same operations. So the reader is promoted into a
new focused module, symmetry.py, as the first tenant of
the module this section already earmarks: a single
read_conv_abc_point_ops(...) that both the kp writer (which
keeps consuming the fractional translations for
buildAtomPerm) and the producer (which uses the rotations
for axis classes, DESIGN 3.12) read through, so one parser of
the spaceDB operation file serves both and cannot drift.
StructureControl gains a thin accessor delegating to it, so
a caller can still ask a loaded cell for its operations
without the reading logic entering the monolith. The
remaining seams (geometry, neighbors, I/O, and the rest of
symmetry -- supercell construction, atom permutation tables)
stay future work; this starts the split in the small, driven
by a real need rather than a scheduled refactor.
The third active prong (VISION Goal 3) augments today's
isolated-atom potential database with potentials extracted
from converged Imago runs on curated reference solids. This
section covers on-disk layout, file format, the lookup path
through makeinput.py, the build pipeline, and the
validation harness. Algorithmic details (interpolation across
nearby labels, environment-descriptor computation) are
deferred to DESIGN.
The augmented database lives inside the existing per-element
directories under share/atomicPDB/. Each element keeps its
own unique store, into which additional files may be added
over time without disturbing the rest of the tree.
share/
atomicPDB/ Per-element database directory
au/
pot1 Existing: nuclear/electronic metadata
coeff1 Existing: Gaussian coefficients + alphas
s_gaussian_pot.toml New: augmented database, multiple
labeled potentials per element
o/
pot1
coeff1
s_gaussian_pot.toml
ti/
pot1
coeff1
s_gaussian_pot.toml
...
pot1 and coeff1 remain the canonical atomSCF output
for the isolated-atom potential. The new per-element file is
a derived artifact regenerated from atomSCF output plus
converged Imago runs on curated reference solids (8.5). It
contains the same isolated-atom data as one labeled entry,
plus one or more additional labeled entries from the
solid-state runs.
The producer's own working artifacts live apart from the
databases, under a share/curation/ tree:
share/
curation/ Producer working artifacts. Every
file here is reconstructible by
re-running the producer.
structures/ Materialized reference structures:
the CIF fetched from COD and the
skeleton converted from it, keyed
by reference_id (DESIGN 5.7)
workspace/ The kaleidoscope flight workspace
the producer dispatches into (9.6);
its run-reuse cache persists across
producer runs
run_log.toml One record per reference solid of
the last run -- converged mesh,
k-point density, SCF iteration
count -- read by the validation
harness (8.6)
The split is along reconstructibility. A harvested entry
in s_gaussian_pot.toml is the product of cluster time and
cannot be recovered by re-running a script; everything under
curation/ can. Keeping the two apart means the campaign
footprint can be cleared wholesale, in one gesture, without
anyone having to reason about which files were expensive.
Why the structure cache is not under atomicBDB/. It
once was, at share/atomicBDB/cache/structures/, placed
beside a per-solid SCF cache that DESIGN 5.7 later dropped
in favour of kaleidoscope's run-reuse cache (9.6). Once
that sibling was gone, the surviving directory left
downloaded crystal structures sitting inside the atomic
basis database, where they are neither per-element nor
basis data and where no reader would look for them.
Grouping them with the producer's other working artifacts
keeps atomicBDB/ purely per-element.
The augmented database uses TOML. Files are small (well
under 1 MB per element), human inspection of checked-in
files is a hard requirement, and TOML reads via the Python
stdlib (tomllib). Writes are hand-formatted from inside
the helper module (initial_potential_db.py) -- the subset
of TOML emitted is small enough that a focused emitter is
simpler than pulling in a third-party writer dependency.
Architectural invariants the format enforces:
schema_versionon every file. Lets the reader refuse unknown versions and lets future schema changes be additive without breaking older files. Bumped from v1 (Phase 1) to v2 (Phase 2) to carry the per-entry default tag and the per-entry fingerprint sub-blocks.- Provenance is required, not optional. Every numerical-potential entry must carry enough information to retrace its origin -- source structure, atom site, Imago commit, convergence parameters, SCF iteration count. Non-negotiable per VISION Principle 5 and feeds the validation harness (8.6).
- Geometric-alpha layout preserved but explicit.
alpha_min,alpha_max, andnum_gaussiansdefine an implicit geometric series;alphas[]is written explicitly to allow future non-geometric layouts without a format change. - Numerical content is the same as the legacy
pot1+coeff1pair (radial Gaussian fit of the nuclear and electronic potentials); the TOML file wraps it in a metadata envelope.
Full schema (top-level keys, per-entry keys, fingerprint sub-blocks, provenance, validation rules), the deterministic emitter contract, and a worked sketch live in DESIGN §5.2, §5.3, and §5.5. This section names the architectural choice; DESIGN owns the bytes.
A future arrays-grow-long-enough-to-hurt-readability contingency (hybrid TOML metadata + sidecar columnar block referenced by relative path) is the format-level escape hatch. The reader/writer abstraction in §8.7 is designed so that swap is local.
share/atomicPDB/<elem>/{pot1,coeff1} (atomSCF output;
canonical
isolated-atom data)
|
| ingested at regeneration time (8.5)
v
share/atomicPDB/<elem>/s_gaussian_pot.toml
(augmented database;
"isolated" entry
mirrors atomSCF;
extra labeled
entries from
curated Imago runs)
|
| makeinput.py reads, picks a label per atom
| site, emits the chosen coefficients into the
| Imago input file in the existing on-the-wire
| format
v
Imago input file (unchanged format)
|
v
Imago SCF (Fortran consumes input identically regardless of
database source)
The Fortran side does not change. All format awareness lives
in makeinput.py. The augmentation is a Python-side feature
end to end.
A "label" is an arbitrary string key under which one
numerical potential entry lives within an element's
database file. The format imposes no naming convention
on labels beyond one reserved value: "isolated", which
is the atomSCF-derived baseline and must be present per
DESIGN §5.2 rule 6. Every other label is curator-chosen
("default_solid", "fcc_bulk_metal", "tetrahedral_O",
etc.) and conveys meaning only to humans reading the file.
Selection -- the question of which labeled entry to use for each atom -- happens through two cooperating mechanisms:
- The
defaulttag (DESIGN §5.2 rule 7). Exactly one entry per file carriesdefault = true. This is the entry picked for any atom that no scheme grouped by local environment. - Fingerprint records (DESIGN §5.2 fingerprint
sub-blocks). Each entry optionally carries one or
more
[[potential.fingerprint]]records keyed by(method, sub_spec). When the user requests an environment-based grouping (-reduce,-bispec), the matcher protocol (§8.9) computes a per-atom fingerprint, buckets atoms into species by descriptor similarity, and picks the manifest entry whose recorded fingerprint best matches each species.
Phase 1 ships only the literal-label override
(-pot LABEL applied uniformly across the structure)
plus the default tag fallback. No fingerprint records
are required in Phase-1 files; entries that lack them
simply don't participate in fingerprint matching.
Phase 2 layers the matcher-driven path on top of
Phase 1 without changing the legacy override: -pot LABEL
still wins when given. The matcher protocol is the
extension point for new descriptor families (element-
aware bispectrum, SOAP, future schemes) -- adding one is
a new matcher class in makeinput.py plus a new
[[potential.fingerprint]] shape, with no schema
rewrite.
Full algorithm (CLI surface, mutual-exclusion rules, spatial scoping, per-species pick, type inheritance) is in DESIGN §5.6. The Phase-3 interpolation question -- what to do when the best fingerprint match exceeds the matcher's similarity floor -- is parked in DESIGN §5.9.
The augmented database is a build product, not a hand-edited artifact (VISION Principle 5).
src/scripts/
build_initial_potentials.py Curates and regenerates the
augmented per-element files
under share/atomicPDB/<elem>/
The producer is a kaleidoscope client. It does not run reference SCFs itself. For each curated reference solid it builds a small predict-then-verify flight, hands that flight to kaleidoscope to dispatch and track, and then harvests the converged potential from the run directories. This is the relationship set out in 9.7: the producer is a what-to-compute script, while kaleidoscope owns running, caching, and tracking the batch.
The producer keeps no cache of its own. Kaleidoscope's run-reuse cache
(9.6) serves it, keyed on structure.dat -- makeinput's resolved
output, which bakes in the type/species assignment, basis, functional,
and potential -- plus the SCF convergence limit. Because an entry
customization touches none of those, editing one re-triggers no SCF:
only the cheap harvest step re-runs. The engine build is recorded
alongside each run but is not part of the key (9.6, VISION 16), so
rebuilding the engine costs nothing until a curator asks for the re-run.
Prepare runs in the driver, not the worker. Keying
on structure.dat means makeinput must run before the
hit-test can read it (DESIGN 6.2.5). So the producer
prepares each unit itself, in the driver: it runs
makeinput into a per-unit staging area -- separate from
the run directory, so the prior run's structure.dat
survives untouched for the byte-compare -- and points the
unit's key-file source at the staged copy. Kaleidoscope's
cache core is unchanged: it still byte-compares the unit's
key-file source against the run directory's copy, knowing
nothing about makeinput. The unit carries the staging
reference as opaque domain data the core round-trips but
never interprets (Principle 9), exactly as it does the
prediction record and the makeinput options. On a cache
hit the decision is made in the driver from local files,
so the unit never reaches the scheduler; on a miss the
wingbeat commits the staged inputs into the run directory
and runs the prepared directory (no rebuild). Where the
driver itself runs -- a login node or its own batch job --
is settled in DESIGN 6.2.11.
Open item: decoupling the potential harvest from the
producer. Running the calculations is already a
separate concern from the producer (kaleidoscope owns
it, above). Reading the results back is only partly
separated. The convergence-detection helper and the
guidance-staging step already live in a standalone
script (guidance_harvest.py, DESIGN 7.8), which points
at a finished workspace and produces staged output
without re-running anything. The potential harvest --
walking each converged run directory to extract the
per-environment potential and its fingerprints into the
per-element database -- is not separated: it runs
in-process inside build_initial_potentials.py (DESIGN
5.7 step 5). Giving that harvest the same standalone
shape (a tool that reads a finished workspace and writes
or updates the database) would let a finished result set
be re-read by different or newer harvest tools without
re-invoking the producer. This is a genuine
architecture change, deliberately not undertaken
yet; it is recorded here so the trade-off is visible
before anyone acts on it.
- For. Re-harvest a finished result set with a new or
different tool without touching the producer. Symmetry
with
guidance_harvest.py, which already demonstrates the standalone-post-processor shape works. A cleaner split of concerns -- producer = "make the calculations exist and converge," harvester = "read results into a database" -- each testable on its own. The[harvest]manifest block (DESIGN 5.7) is exactly the interface such a standalone harvester would read, so it is already in place. - Against (the stronger case, hence no action yet).
The potential harvest reads far more context than the
guidance harvest does -- entry customizations, derived
labels, the dedup rule, the
datSkl.mapsite mapping, the resolved per-solid settings -- all of which the producer already holds in memory while it dispatches. A standalone harvester must reconstruct every one of those from the workspace plus the manifest, which is real work and a real surface for drift. And the urgent need it would serve is already met: because of the run-reuse cache and the regenerate model (9.6), a producer re-run over a warm workspace skips every calculation and goes straight to harvest, so "re-harvest with a new recipe" works today -- it is simply expressed as re-running the producer, not as a separate tool. Decoupling is therefore an ergonomics and separation-of-concerns improvement, not a missing capability, and it is not worth the reconstruction cost or the drift risk until a concrete need pushes past the warm-cache path.
Inputs:
- A curation manifest in TOML, schema v2 (full spec in DESIGN
5.7). It declares the reference solids, each solid's
system_type, a database-wide[characterization]fingerprint recipe, and optional per-entry customizations. The harvest is automatic -- one representative per distinct environment (DESIGN 5.2.2/5.2.3) -- so the manifest does not enumerate atom sites; a customization, when present, only annotates an auto-discovered environment (itsdefaulttag, description, label, or a pinningatom_site). - The historical guidance dataspace (section 10), consulted per reference solid to predict the converged k-point density and the width of the verification grid around it. When the dataspace is too sparse to predict, the producer falls back to a wide default grid.
- The existing
pot1/coeff1files (for"isolated"entries). - An Imago build, used through kaleidoscope to run the
verification flights and, for Fortran-side fingerprint
matchers, the follow-on
imago.py -loen -scf noruns per Fortran-side characterization fingerprint.
Structure materialization. A reference solid is named
either by a local structure_path or by a cod_id pinned to
a cod_revision. A thin materialization step turns either
form into one local structure file: a structure_path is read
from disk; a cod_id is fetched once from the Crystallography
Open Database at its pinned revision and written to a plain
local location. This is the producer's only network access,
and it is deliberately decoupled from any cache -- its sole
job is to hand kaleidoscope a local structure file for its
run-reuse cache to key on. Fetch failures are strict (named
error, no silent fallback to another revision).
Manifest authoring. A complete manifest is not
hand-written from scratch. cod_fish.py (9.5) discovers and
pins structures from COD and, by default, prints a complete,
runnable manifest -- a schema_version header, the shared
[characterization] recipe and [defaults] run settings, and
one [[reference_solid]] stub per structure whose reference_id
is auto-derived from each CIF's metadata -- so cod_fish.py pin <ids> > manifest.toml is the whole authoring step for the common
COD case; the recipe and default values come from the shared
library, not hardcoded in cod_fish. Its --sketch-only mode
prints the bare stubs instead, for the cases that still want
expand_manifest.py: local (non-COD) structures, a non-default
recipe, or interactive per-structure customization.
expand_manifest.py reads such a sketch (from a file, or from
standard input when none is named) and fills in the
[characterization] recipe and [defaults] run settings -- the
same shared-library values -- writing the finished manifest
(per-structure customizations optional, added interactively or by
hand). The producer reads only a finished manifest. The schema
-- the dataclasses, the strict and relaxed readers, and the
writer plus its default-recipe and default-run-setting helpers --
lives in the shared leaf library curation_manifest.py, imported
by both authoring tools (to write a manifest) and the producer
(to read one), so they cannot drift. The relaxed
reader (load_structure_sources) also backs the producer's
--materialize-only pre-flight, which fetches and converts every
reference structure before the run and harvest fields are filled
in.
Outputs:
- Updated augmented database file in each affected
share/atomicPDB/<element>/directory (grown in place; see the Incremental property below). - A guidance contribution. The same converged grid point that feeds a potential entry is also harvested into the historical guidance dataspace's staging area (section 10), so every reference solid the producer converges sharpens the predictor for the next one.
- A run log capturing SCF iteration counts and convergence metrics per reference run (input to 8.6).
Properties:
- Reproducible at supported precision. The TOML
emitter is deterministic at the bit level: given a
fixed in-memory database, it always writes the same
bytes. The pipeline's numerical outputs are
reproducible only to within SCF / Gaussian-fit
precision -- not bit-identical -- because
floating-point accumulation order, threading, and
external library versions can shift the last few
bits, and development changes to the SCF or fit
code that improve converged potentials are expected
to perturb the numbers. Provenance metadata (e.g.,
generated_attimestamps) refreshes on each run and is exempt from any reproducibility guarantee. - Incremental, on two levels. The database grows in place: the producer loads each existing per-element file, appends the environments new to it, and skips those already present (DESIGN 5.2.3), so several manifests accrete into one database and re-running an unchanged manifest changes nothing. The compute is incremental too: reusing kaleidoscope's run-reuse cache (9.6), a run re-executes SCF only for reference solids whose structure, options, or build identity changed, and adding harvest declarations to an unchanged solid re-runs only the harvest.
- Scriptable. Invokable from CI or a developer workstation without manual intervention.
VISION Principle 7 sets the headline metric: improved- initial runs must average at least 20% fewer iterations than isolated-initial runs across a representative benchmark.
src/scripts/
bench_initial_potential.py Runs the benchmark set
under both "isolated" and
"default_solid" labels,
records iteration counts,
emits a comparison report.
The benchmark set should partially overlap with the curation manifest (validating that the database performs on the solids that produced it) and partially extend beyond it (guarding against overfitting). Both sets are checked into the repo as part of the deliverable. The harness report is the artifact that decides whether the 20% target has been met; it is regenerated whenever the augmented database is regenerated.
Python:
makeinput.py: gains a database reader (the format helper from below), a CLI argument that selects the global default label, and the augmented-vs-legacy lookup path. Legacypot1/coeff1lookup remains for backward compatibility but is no longer the default once the augmented database is populated.build_initial_potentials.py: new (8.5).bench_initial_potential.py: new (8.6).
The database reader/writer is encapsulated in a small
helper module (proposed src/scripts/initial_potential_db.py)
so that any future format swap or schema-version bump is a
one-file change rather than a spread of edits across the
three scripts above. The version table and the per-bump
field derivations that DESIGN 5.2.5 specifies live there
too, for the same reason.
src/scripts/potential_migrate.py: future bulk schema-migration tool, rewriting an installedshare/atomicPDB/tree in place. Mirrorsguidance_migrate.py(10) andresource_migrate.py(11), and like them is not in day-1 scope: the reader's in-memory migration plus the producer's next save already carries forward every file a run touches (DESIGN 5.2.5), so this tool is only needed for files no producer run will revisit.
The curation-manifest schema -- its dataclasses, the strict
load_manifest_v2 and relaxed load_structure_sources readers,
and the format_manifest writer -- is likewise encapsulated in a
small leaf library, src/scripts/curation_manifest.py, imported
by the producer and both authoring tools; it depends only on
the lower libraries it validates against (initial_potential_db,
guidance_db) and carries the shared default-recipe and
default-run-setting helpers the authoring tools emit. New
scripts: curation_manifest.py (the schema library) and
expand_manifest.py (the sketch-to-manifest authoring tool,
8.5); cod_fish.py gains complete-manifest emission (8.5).
Fortran: no changes for the Phase-1 chain or the Phase-2
base chain. Imago consumes the same input file format;
the choice of initial potential is opaque to it. The
Phase-2 follow-up (element-aware bispectrum, TODO C62)
requires extensions to O_Input (a new bispecByElement
parameter) and loen.f90 (per-neighbor-element
accumulation in computeBispectrumComponent plus an
extended fort.21 output format). Those extensions are
scoped under DESIGN §5.9 and TODO D10 / C62 -- they are
not part of the Phase-2 base deliverable.
Shared modules: no changes. The new files live inside the
existing share/atomicPDB/<elem>/ directories alongside
pot1 and coeff1; no new top-level directory is
introduced.
Flagged here so DESIGN can pick them up rather than litigating them inline:
- Curation manifest format. Resolved (2026-05-11). TOML,
schema v1 (Phase 1), bumped to v2 in Phase 2 to carry the
default tag, the database-wide
[characterization]fingerprint recipe, and optional per-entry customizations. Full schema, validation rules, cache layout, and COD-fetch contract are specified in DESIGN 5.7. TOML was chosen so the manifest reader uses the sametomllibstdlib machinery as the per-element database file. - Phase 2 selection / interpolation method. Selection
resolved (2026-05-19; decoupled from grouping 2026-06-17,
C93); interpolation parked for Phase 3. The per-species
database pick (DESIGN 5.6.5) runs for every species,
independent of how it was grouped -- crystallographic,
position-based (
-target/-block), or environment-based (-reduce/-bispec) -- and is default-on, with-nofingerprintthe opt-out. The pick forms one order-independent representative per species and matches it against the element database: when the user chose an environment scheme it matches that family at the user'ssub_spec; otherwise it reads the database'spreferredrecord per family (5.2.2, the[characterization]convention). The interpolation question (what to do when the best match exceeds the similarity floor) is parked for Phase 3; Phase 2 falls back to the default-tagged entry with a warning. Full algorithm in DESIGN 5.6, parameter mapping in 5.10.
The Phase-2 selection algorithm dispatches on a small abstraction called a matcher. Each matcher knows exactly one descriptor family (e.g., reduce shell-codes, bispectrum components) and exposes a uniform interface so the species pass and the producer can call it without caring which family it is.
Location. The matcher protocol and its concrete
implementations live in src/scripts/matchers.py, a
neutral library module (not a CLI script, so it does not
add to script proliferation). The protocol first lived
inside makeinput.py, but moved here once makegroups.py
needed BispecMatcher: keeping it in makeinput would
have forced makeinput to import upward from makegroups
(a cycle). Every caller now imports it downward --
makeinput.py (ReduceMatcher / ReduceStructureView
for the reduce grouping), makegroups.py (BispecMatcher
for the sequential loen flow), and
build_initial_potentials.py (the MATCHERS registry).
This also stages the eventual split in which makegroups
owns all environment-based grouping.
Protocol surface. Each matcher class exposes:
name Matcher identifier
(the string written
into manifest
method fields).
E.g., "reduce",
"bispectrum".
needs_loen_run (bool) True for matchers
whose descriptor is
computed by Imago
(loen path), false
for matchers that
compute in Python.
True means the
descriptor is obtained
by the sequential loen
flow that makegroups
orchestrates (DESIGN
5.10), not in-process.
default_similarity_floor (float) Per-matcher default
distance threshold
used by the
fingerprint-match
step (DESIGN 5.6.5).
Atoms whose nearest
manifest fingerprint
is farther than this
fall back to the
default tag. Users
can override per
scheme on the CLI.
to_loen_input(sub_spec) Translates a
sub_spec inline
table into the
parameter dict the
LOEN block of
imago.dat
expects. Only
meaningful when
needs_loen_run is
true.
parse_loen_output(path, Reads the enriched,
sub_spec) self-describing
fort.21 of a loen run
(DESIGN 5.10.3) into
per-site records --
each carrying the row's
identity (element,
species, type) plus its
fingerprint vector.
Only meaningful when
needs_loen_run is
true.
compute_query(structure, Python-side matchers
sub_spec) only: computes
per-atom fingerprint
vectors from
StructureControl
in-process. Loen-side
matchers do not
implement it -- their
vectors come from the
sequential loen flow
(DESIGN 5.10), read off
fort.21 by
parse_loen_output.
distance(vec_a, vec_b) Symmetric scalar
distance in the
matcher's descriptor
space, used for both
species bucketing
and manifest-entry
selection.
representative(members) Reduces a list of
member-atom
fingerprints into one
representative
fingerprint per
species (DESIGN
5.6.5 step 2).
BispecMatcher
returns the
element-wise mean;
ReduceMatcher
returns the first
member's shell-code;
future matchers may
use a medoid or any
scheme appropriate
to their descriptor
space. The protocol
does not pin the
choice, only the
shape (members in,
one fingerprint
out).
build_payload(vector) Wraps one fingerprint
into the dict stored
as a
FingerprintRecord
payload (DESIGN 5.2,
5.4): BispecMatcher
emits {values: ...},
ReduceMatcher emits
{shell_code: ...}
(element-only, so the
descriptor transfers
across structures).
The producer (the
fingerprint harvest)
calls it.
extract_query_vector(payload) The inverse: reads a
stored payload back
into the form this
matcher's distance
expects, so producer
and consumer agree on
field naming.
Registry. matchers.py defines a module-level dict
mapping matcher names to matcher classes (imported by every
caller -- makeinput.py, makegroups.py, and
build_initial_potentials.py):
MATCHERS = {
"reduce": ReduceMatcher,
"bispectrum": BispecMatcher,
}Adding a new descriptor family (e.g., element-aware
bispectrum, SOAP) is a new class plus a new dict
entry; no other code needs to change. The
initial_potential_db.load() validator consults
MATCHERS.keys() to enforce per-element-database
rule 9 (unknown method is a hard error).
Concrete Phase-2 matchers.
ReduceMatcher(needs_loen_run = false). Wraps the existing reduce algorithm (DESIGN 5.6.4); the current in-place implementation ingroup_reduceis refactored to live behind this class so producer and consumer reach it through the same surface.BispecMatcher(needs_loen_run = true). Mapssub_spec = {twoj1, twoj2}to the LOEN input block; parsesfort.21rows into per-site records whose vector hastwoj2 + 1components (the coupling channelsjin|j1 - j2| <= j <= j1 + j2, withtwoj1 >= twoj2). Element-aware mode is gated by an optionalby_elementkey insub_spec, currently ignored (Phase-2 follow-up; DESIGN 5.9 and TODO).
The pick is a per-species step, decoupled from
grouping. The fingerprint-based potential pick is not a
side effect of an environment grouping flag: it runs for
every species, however that species was grouped --
crystallographic, position-based, or environment-based
(the C93 decoupling; DESIGN 5.6.4/5.6.5). It is on by
default and disabled per run by -nofingerprint. Which
matcher and parameters the pick uses split by regime:
-
User grouped with an environment scheme (reduce in makeinput; bispectrum is grouped upstream by
makegroups). That family and the user'ssub_specdrive the match, reusing the per-atom descriptors the grouping pass already computed. The database never overrules the user: a database lacking thatsub_specis a silent best-effort miss to the default entry, never an error. -
Species are file-dictated (crystalline or pre-assigned). The matcher carries no user
sub_spec, so the database supplies it: each element's database marks onepreferredfingerprint record per family, and the preferredsub_specfor a family is uniform across the whole database (the curation convention, DESIGN 5.7). The pick reads the preferred record -- bispectrum if present, else reduce -- computes exactly one query, and matches. No search across families or sub_specs, and at most one loen run per structure.
This adds one protocol obligation beyond grouping: a
needs_loen_run matcher must be able to produce a query
for the pick's file-dictated branch. It does so through
the same loen seam used for grouping (to_loen_input ->
loen run -> parse_loen_output), reusing a makegroups
run when one exists and otherwise paying one fast loen run.
The bispectrum compute is cheap enough that the consumer
simply pays it; optimizing it away is not warranted unless
it ever becomes a measured problem.
Why this is architectural, not just design. The
matcher abstraction is what isolates the Imago Fortran
side from the manifest schema's growth. Adding a
descriptor family does not touch imago.f90 unless
the family needs new loen capability; conversely,
adding loen capability (element-aware bispectrum, a
SOAP path) does not touch the manifest schema or the
selection algorithm. Listing the protocol here pins
that contract down before code lands.
This section specifies the infrastructure of VISION Goal 4: shared machinery for submitting, tracking, and harvesting batches of Imago calculations. Its purpose is that the scripts which decide what to compute -- the initial-potential database build (Section 8.5), convergence sweeps, validation harnesses, and future ab-initio molecular dynamics and high-throughput screening -- do not each reinvent how to submit and watch jobs on a cluster. The flight dispatcher is named kaleidoscope.
The layers and the four scripts/packages that realize them:
imago.py-- runs one calculation per invocation; both a CLI and a callable Python API (9.2).ase_imago.py-- an ASE Calculator,ImagoCalculator, adapting Imago to the materials-simulation community (9.3).kaleidoscope/-- a Parsl-based flight dispatcher that drives many calculations in parallel and/or series batches (9.4).cod_fish.py+cif2skl.py-- the structure- acquisition front-end: pull CIFs from the Crystallography Open Database, convert them toimago.skl(9.5).
Three VISION principles are the load-bearing constraints. Principle 8 (decouple adapter from orchestrator) keeps 9.3 and 9.4 independent. Principle 9 (domain-specific machinery lives at the adapter layer; the flight layer is ordinary scientific Python) keeps kaleidoscope free of materials-specific coupling. Principle 10 (complete-and-report at the flight level) means one failed calculation never fails the whole flight by default.
The dependency arrows point in one direction only, so each layer can be replaced without disturbing the others:
cod_fish.py --> cif2skl.py --> imago.skl
|
kaleidoscope/ --drives--> makeinput.py --> inputs
| |
+----------------- calls ---------> imago.py (API)
^
|
ase_imago.py (ImagoCalculator)
imago.py is the foundation: it has no dependency on
ASE, on Parsl, or on the flight layer. The ASE
adapter and kaleidoscope both sit above it and call
down into it. This is what lets the cluster jobs and
kaleidoscope run without ASE installed, and lets a
future orchestrator (Snakemake on top, a different
dispatch backend) slot in without touching imago.py
or the adapter.
Today imago.py is a command-line driver: it stages
prepared input files, selects the gamma vs. non-gamma
executable, runs the Fortran binaries, collects and
renames outputs, checkpoints completed work, and
manages a lock file. imago.py on its own performs a
plain SCF.
The refactor exposes that same orchestration as a
callable Python API, with the existing CLI reduced to a
thin wrapper over it. The name stays imago.py (not a
separate run.py): the module is new, so there is no
namespace collision to avoid, and one obvious name for
"run an Imago calculation" is better than two.
The API offers two entry granularities, so a caller can join at whichever level it already has inputs for:
- Prepared-directory mode -- given a run directory
that already holds the Imago inputs (
imago.dat,structure.dat,scfV.dat, kp files), run it as-is. No makeinput call; the caller (or a prior step) produced the inputs. - Structure-and-options mode -- given a structure and
a set of makeinput options, drive
makeinput.pyto build the run directory first, then run it.
Both return the same small result object
(success/failure, SCF iteration count, paths to the
output files such as the converged
gs_scfV-<basis>.dat). Input preparation still lives
in makeinput.py; structure-and-options mode simply
calls it on the caller's behalf.
This API is the single seam every higher layer reaches through, so its contract is deliberately Imago-native and dependency-free.
ImagoCalculator subclasses ASE's Calculator. It
translates between ASE conventions (an Atoms object,
a results dict keyed by ASE property names, eV/Angstrom
units) and Imago's native world (an imago.skl
structure, makeinput options, Hartree/Bohr), and runs
the calculation by calling the 9.2 API.
It is a separate module, not folded into imago.py,
for three reasons. (1) Dependency isolation: ase is
a third-party package; keeping it out of imago.py
means the callable API, the cluster jobs, and
kaleidoscope do not require ASE to be installed. (2)
VISION Principles 8 and 9 treat ASE as a swappable
adapter, with domain-specific machinery confined to
this layer. (3) Contract mismatch: the native API is
Imago-shaped while a Calculator must be ASE-shaped;
separating them lets each stay clean.
The module is named ase_imago.py, not ase.py: a
flat module named ase.py on the scripts path would
shadow the real ase package on import. It stays a
single flat module: ImagoCalculator is one class that
exposes many ASE properties (energy, forces, stress,
charges, dipole, and Imago specialties such as DOS,
bands, and bond order as custom result keys) and plugs
into ASE's optimizer, MD, and analysis ecosystem -- all
of which is one calculator, not many entry points.
It would be promoted to an imago_ase/ package only if
auxiliary adapter modules later accrete (a dedicated
DOS/bands adapter, an ase.db writer), the same
trigger-based split policy as Section 7.
Committing to ASE buys the materials-community interoperability VISION Goal 4 describes -- LAMMPS, ASE's MD integrators for future AIMD, and acceptance by ASE-consuming workflow tools.
The Atoms -> Imago structure translation splits along
the ASE boundary, so that no ASE dependency leaks into
the core:
structure_control.pygains an ASE-free factory that builds aStructureControlfrom plain arrays (lattice vectors, fractional coordinates, element symbols). This is the genuinely shared, broadly reused primitive; it has no knowledge of ASE.- The small ASE-specific step -- reading those
arrays off an
Atomsobject -- stays here at the adapter layer.
This matters because structure_control.py is imported
by nearly every script (Section 7); putting ASE-domain
machinery there would risk making the whole toolchain
ASE-dependent the first time someone added an import ase for an isinstance check or to construct an
Atoms. Keeping the factory ASE-free and the
Atoms-reading glue at the adapter layer honors
Principle 9 and keeps import ase confined to the
modules that genuinely need it. cif2skl.py (9.5)
reuses the same factory.
kaleidoscope/ is a Parsl-based package that drives a
set of calculations. Given a flight specification
(which structures, with which makeinput options), it
dispatches the per-structure work -- makeinput to build
inputs, then the 9.2 imago.py API to run -- across
SLURM via Parsl, handling both embarrassingly parallel
sweeps (thousands of independent SCFs) and tightly
iterative inner loops (adaptive convergence, future
AIMD) under one model.
Kaleidoscope dispatches a unit of work through a
pluggable wingbeat seam. The default wingbeat is the
imago.py API directly -- the flight layer is
ordinary scientific Python (Principle 9), and most
flights (the database build, convergence sweeps) need
nothing more. But a unit may also be dispatched
through the ASE adapter (for ASE-MD or ASE relaxation
semantics), or through a future adapter, and a single
flight may blend wingbeats -- some units plain Imago
SCF, others adapter-wrapped -- so that Imago
calculations can be mixed with other ASE-compatible
calculations and dispatch activities under one flight.
Keeping the wingbeat pluggable is what lets new adapters
and new blends slot in without changing kaleidoscope's
dispatch core (Principle 8).
Per Principle 10, kaleidoscope records and surfaces per-job outcomes (converged, non-converged, cluster-side loss, post-processing error) but does not abort the flight on a single failure; the client script decides whether the aggregate result is acceptable for its scientific purpose. It is a package (not a flat module) because it carries real substance: dispatch, status tracking, harvest hooks, and workspace management (9.6).
Cluster dispatch is configuration, not code. The
mapping onto SLURM lives entirely in the flight's Parsl
Config (DESIGN 6.2.3): the same dispatch core serves a
laptop, an interactive node, and a batch allocation, and
only the Config changes. Two cluster topologies are
both supported because they suit opposite regimes (VISION
Goals 4 and 7):
- One shared allocation (pooled). Parsl requests one (optionally auto-scaled) block of nodes and streams many units through its workers. Best for many small, similar units -- the convergence sweeps and the database seed.
- One scheduler job per unit. Each unit is its own SLURM submission. Best for large or heterogeneous units -- future MPI/GPU solves of differing size -- where a single uniform worker slice would waste resources.
No client hand-writes a SLURM script. A Config is
assembled from three layers: a per-site resource-control
file (queues, account, per-node cores and accelerators,
and the commands a worker runs to bring up the imago
environment), a few per-run choices (topology, partition,
nodes, walltime), and the per-unit resource request that
the resource-and-cost dataspace (section 11) predicts. A
generator fills the defaults and surfaces only the
decisions that are genuinely the user's. Right-sizing
heterogeneous parallel units -- routing each to a block
sized for it -- needs per-size executors keyed on a cost
hint and is staged with section 11, not built up front.
How those layers are settled -- the rc-file home, the CLI choices,
the right-sizing deferral, and the generator's home -- is recorded
in DESIGN 6.2.11, resolving the questions once gathered in 9.8.
The resource-control file is about the cluster, plus a bounded
set of job-class defaults. Most of the per-site file is cluster
fact -- the queues, account, per-node cores and memory capacity,
and the environment bring-up -- true no matter what runs. The
remaining settings size a job, and there are only three job
CLASSES on the cluster, not one per builder: a worker (a
single calculation), an orchestrator (a driver that
prepares units and fans them out to workers, DESIGN 6.2.11), and
an md job (an external molecular-dynamics program run under
MPI as many ranks on one node). A future producer that runs as a
driver is an instance of the orchestrator class -- it reuses the
one orchestrator shape, overriding on the command line if it
differs -- so the file grows by job class (bounded), never by
orchestrator instance. Each class carries a single site-default
resource shape (the worker sizing, and small orchestrator and
md groups), overridable per run. These hand-set shapes are
transitional: as the resource-and-cost dataspace (section 11)
learns to predict a job's memory and walltime from its problem
size, the per-job requests migrate out of the file toward
prediction, leaving it closer to pure cluster fact -- which is
why per-node memory is a capacity ceiling, not a request.
The file is the site's scheduler truth, not kaleidoscope's
alone. The md class exists because condense.py writes a
LAMMPS submission file, and a script that emits a batch job needs
the same queues, account and per-node facts a flight does.
Recording those separately in each script's own *rc.py would
leave one fact written down in several places with nothing to say
which copy had gone stale, so every Imago script that submits
work reads the cluster facts from this one file, and keeps in its
own *rc.py only what is not about the cluster.
The reader that performs that load sits in kaleidoscope beside the dispatch builders, and imports no Parsl: the flight layer defers that import to the point of use. A script may therefore share the reader without taking on the dispatcher's requirements, and a LAMMPS input builder reads the site file without depending on a flight layer it never calls. That is a property to keep deliberately rather than an accident to rely on -- it is what lets the two submission-file generators live together, where a reader looking for one will find the other.
The md class carries its own bring-up rather than reusing
worker_init, because worker_init starts imago and an external
MD program does not need it. That separation is what keeps the
install location of a program Imago merely launches out of
src/ entirely: a site points the md bring-up at wherever it put
LAMMPS, and no path to it appears anywhere in the source.
Structure acquisition is a front-end, separate from the
flight dispatcher, in two small CLI tools. Both lean on
existing services rather than reinventing them: COD's own
HTTP query endpoint for search, ASE for CIF parsing, and
Imago's own apply_space_group for symmetry expansion.
Neither tool adds a dependency beyond urllib (standard
library) and ase (already required to read CIF).
cod_fish.py -- COD acquisition front-end. Four
verbs, all thin wrappers over COD endpoints; it keeps no
local database and writes no SQL of its own:
get <cod_id> [--revision REV]fetches one structure's CIF viaurlliband is strict: with--revisionit verifies the CIF's$Revision:$and refuses on a mismatch; a network/COD failure or a missing pinned revision errors rather than silently substituting a different revision (the reproducible path the producer and a manifestcod_iduse; matches DESIGN 5.7).search --elements E [E ...] [--text AUTHOR] [--max-extra K]queries COD'sresult.phpand prints a numbered candidate table (cod_id, formula / mineral / common name, spacegroup + number, the six cell parameters, volume). Exact composition is the default: the query sets the distinct-element count bounds to the length of the element list (strictmin = strictmax = len), so--elements Sireturns single-element Si phases, not every Si-containing compound;--max-extra Kloosens the upper bound to admit up to K further elements. Search also writes the result set to a small session file in the working directory so later verbs can refer to rows by their index rather than their eight-digit id.pin <index|id> [...]resolves the chosen rows' revisions (fetching only those few) and prints a complete, ready-to-read sketch manifest: aschema_versionheader plus one[[reference_solid]]stub per structure, each withcod_id/cod_revisionfilled in and areference_idderived from that CIF's own metadata --<formula>_<H-M symbol>_<IT number>_<year>, e.g.si_fd-3m_227_2010. The year both dates the entry and separates phases that share a space group (4H and hcp silicon are both P 6_3/m m c); a residual name clash gets a trailing counter so every id stays unique (manifest rule 5). Each stub also carries two discovery hints read from the same CIF -- the composition (elements) and asource_description(chemical name + space group + year) -- which the authoring tool uses to auto-fill each customization's element and description so the curator invents neither. The two part ways in the finished manifest:elementsis a transient, non-schema hint the producer ignores and the finished manifest omits, whilesource_descriptionis persisted as areference_solidfield (DESIGN 5.7), the structure-level description the harvest qualifies per environment. Socod_fish.py pin <ids> > sketch.tomlwrites a sketch the authoring tool reads directly. This is the bridge from browsing to a reproducible pull; indices come from the saved search, so a student never retypes long ids, and the auto-named id makes a wrong pick obvious.rank(orsearch --rank) is advisory triage, not validation: when a composition has many COD entries it annotates and orders them by interpretable signals drawn only from the returned metadata -- ambient vs high-P/T conditions, presence of a mineral/common name, spacegroup consensus across the candidates, and cell-volume-per-formula-unit outliers -- and prints the reasons so the curator decides. It narrows "which of these is the real phase," it does not answer it.
Discovery (search/rank/the index session) is kept
strictly separate from acquisition (get/pin by pinned
id+revision): a manifest only ever pins an id and a
revision, so reproducibility never depends on a search.
cif2skl.py -- CIF to imago.skl, preserving the
space group. An imago.skl is not a flat P1 list of
atoms: it stores the asymmetric unit plus a space-group
token, and read_imago_skl's apply_space_group
regenerates the full cell. Preserving the space group
matters for correctness, not just compactness -- the
Brillouin-zone integration samples the irreducible wedge
using that space group, so flattening a crystal to P1
would silently corrupt every k-point-dependent result
(the same hazard makegroups.py guards against). So the
converter recovers the CIF's space group rather than
discarding it:
- ASE parses the CIF (
ase.io.cif), which yields the authored asymmetric unit (the raw_atom_site_*list), the cell, the IT number and origin setting, and -- read the ordinary way -- the fully symmetry-expanded cell. This is the only step that touchesase, and it spares us a hand-written CIF parser (CIF'sloop_blocks and symmetry expansion are exactly what ASE already handles). Genuine partial/mixed occupancy is refused for now; a refined occupancy within a small tolerance of 1 (e.g. 0.9999, a rounding on a physically full site) counts as full. - The space group is resolved by verification, not by
parsing operation tables. The
spaceDBentries for the CIF's IT number are the candidate settings (the<IT#>_a,<IT#>_b, ... variants, ordered by ASE's reported origin setting so the likely one is tried first). For each candidate the authored asymmetric unit is written with thatspacetoken and run throughapply_space_group; the variant whose regenerated cell matches ASE's full expansion (atom count and positions, up to lattice translation and permutation, within a tolerance) is the answer. This reuses Imago's own symmetry engine as the judge, so nospaceDBoperation set is ever parsed or canonicalized, and the origin choice falls out of the match (translations are origin-dependent). It is the automated form of the established manual practice -- default to_a, check the resulting structure, try the others if it is wrong. - On a unique match the skeleton is written as the
asymmetric unit plus that token. No candidate
verifying is a hard error that names the IT number
and the tried variants; a
--space <token>override forces a specific setting (still verified). P1 fallback is deliberately not offered for a crystal, since it would reintroduce the corruption above.
No symmetry library (e.g. spglib) is added: the CIF's
symmetry is taken as authored (ASE reports it) and
checked against Imago's expansion, rather than
re-derived. The ASE-free StructureControl factory
(9.3) remains the shared primitive for genuinely P1
inputs and the ASE adapter; cif2skl's symmetric path
writes the asymmetric unit and token directly and
verifies through apply_space_group.
Splitting acquisition from kaleidoscope keeps the flight
layer agnostic about where a structure came from:
kaleidoscope consumes imago.skl files, whether they
came from COD via this front-end, from a hand-authored
structure_path, or from any other source.
A flight that touches hundreds or thousands of
structures needs a directory and naming scheme so the
inputs and outputs do not become an unnavigable mess.
A flight owns a workspace rooted at a single
directory, keyed throughout by a stable per-structure
id (e.g., a COD id or a curation reference_id).
The layout is pinned in DESIGN 6.2.4 (the id
charset/uniqueness rule, the <calc> tag format, and
the status.toml schema); the shape is:
<flight_root>/
flight.toml What to run + global options.
structures/<id>/ Acquired inputs.
<id>.cif
<id>.skl
wingbeats/<id>[/<calc>]/ One working dir per calculation:
<makeinput inputs: imago.dat, structure.dat,
scfV.dat, kp-*>
<run outputs: gs_scfV-<basis>.dat, imago.out>
cache_key.toml Cache identity snapshot.
result.toml Wingbeat-persisted native result.
status.toml queued / running / done / failed /
lost, plus detail and timings.
results/ Harvested / aggregated outputs.
logs/
A single structure may host more than one calculation
(e.g., different bases or property runs), hence the
optional <calc> level under wingbeats/<id>/.
This layout also serves as the producer's content-keyed SCF
cache, so the producer needs none of its own (DESIGN 5.7): a
kaleidoscope wingbeats/<id>/ directory is a cached run,
so cache-hit/miss logic becomes "is there a completed
run directory for this id whose inputs still match?"
The cache is a general kaleidoscope feature, split into mechanism and policy so that generality does not cost correctness:
- Mechanism (kaleidoscope). Write a key snapshot into the run directory, compare it on the hit-test, skip a run whose snapshot still matches, and resume a flight over already-completed runs.
- Policy (client). The client supplies the key fields -- only it
knows which inputs define identity for its calculations. The
database producer declares its key as the SCF convergence limit plus
the resolved
structure.datbytes (DESIGN 5.7): the inputs that settle which calculation this is. Each swept k-density is its own unit with its own run directory, so the mesh never enters the key.
This keeps kaleidoscope from guessing input identity (a too-broad key
risks false hits and wrong results; a too-narrow key risks needless
re-runs) while still giving every flight one cache implementation. The
boundary with imago.py's existing checkpointing is clean: imago.py
resumes within a run directory; kaleidoscope decides whether to
launch the run directory at all.
The key asks "the same calculation?", not "still good?" (VISION 16).
The engine build stays out of it. A rebuilt engine does not make a
stored potential wrong -- that potential is a starting point every later
SCF re-converges -- while keying on the build would discard every stored
result each time the code was touched, which is an expensive answer to a
question nobody asked. So the build identity is recorded -- a field
in the status.toml each run directory already carries, where a curator
can read later what produced a number -- and it is never compared.
It is not a run option either: it reaches neither makeinput nor imago,
so it rides on the unit rather than in the options dictionary those
tools are handed (DESIGN 6.2.10). What stands in for the automatic
guard is a
report: the driver prints its reuse plan before it spends anything --
which units it will reuse, which it will run, and why -- and re-running
is requested outright rather than inferred from a file that changed.
Two tiers, and reclaiming the larger one. A run
directory spans two filesystems. The kept tier is the
directory itself -- staged inputs, result.toml,
status.toml, cache_key.toml, the SCF potential, the
log. The scratch tier sits behind the intermediate
symlink that imago.py creates, pointing at a temporary
area, and holds the engine's working files. Scratch
dominates: measured at seed scale it was 99.7% of the
bytes, almost entirely the HDF5 that carries the
wavefunctions, against a couple of hundred kilobytes of
kept files per calculation.
Reclaiming it is split exactly as the cache is, and for the same reason -- only the client knows when a finished run is finished with:
- Mechanism (the reclamation tool). Walk a root, resolve each run's scratch, preview it, and remove what a policy marks reclaimable.
- Policy (client). Decide when a unit's scratch is
spent. The k-point producer is done with a rung once
its
result.tomllands; a client that post-processes wavefunctions is not done until that step has run.
Unlike the cache, the mechanism sits beside kaleidoscope
rather than inside it. Reclamation reads imago's own names
-- the intermediate link, the lock file, the completion
marker -- out of imago.py, and within kaleidoscope the
one place engine knowledge belongs is the wingbeat (9.4),
the pluggable piece that knows how to run a unit
(Principle 9: domain-specific machinery lives at the
adapter layer). Reclamation is not a way of running a
unit, and the dispatch core beneath the wingbeat names no
imago file at all. So the tool stands on its own, and all
three layers reach it from the client side.
Scratch is not only a flight's problem, though. An
ordinary imago.py run plants the same intermediate
link and leaves the same tens of megabytes behind, and
those runs -- a student's job directory, a hand-driven
convergence test -- are the common case rather than the
exception. Reclamation therefore recognizes two kinds of
root and handles exactly one per call: a workspace,
whose units prove completion in status.toml, and a job
tree, whose runs prove it by holding no imagoLock and
ending their runtime log with the completion marker.
src/scripts/tidy_scratch.py: the standalone reclamation tool -- the home of the logic, which the producer's--clean-after(sweep once the harvest is written) and--tidy-run(prune each unit as it lands) both call rather than reimplementing. The in-flight prune reaches it through theon_outcomehook a flight already fires, so it is a producer wiring and adds nothing to kaleidoscope. Design in DESIGN 6.2.12, including the five refusals that hold for every root (never the run directory, never an unfinished unit, never a link out of the scratch area, never a symlink descended while walking, and never a tree that holds another run's scratch -- scratch mirrors the run path, so a nested run directory nests its scratch too), the two the job tree adds (never a run that has not declared completion, never descending into a workspace from a job tree), and the dry-run-by-default rule.
Kaleidoscope's clients are the what-to-compute
scripts: build_initial_potentials.py (the
initial-potential producer, C48),
bench_initial_potential.py (the validation harness,
C50), and the future AIMD and screening flights.
This reshapes the producer. As currently written,
DESIGN 5.7 / PSEUDOCODE 11.4 / ARCHITECTURE 8.5 have
build_initial_potentials.py run reference SCFs itself,
with its own COD fetch and its own per-solid cache.
Under this section the producer instead becomes a
kaleidoscope client: it hands kaleidoscope the curated
structures and the makeinput options, lets kaleidoscope
run and track the batch, and then harvests the converged
potentials from the run directories. The converged scfV
output carries every potential type in the material (a
NUM_TYPES header + per-type Gaussian blocks under a
TOTAL__OR__SPIN_UP channel); the harvest selects each named
site's type block and takes its coefficients and alphas
together (the site's type number comes from the datSkl.map
described below; the producer runs non-spin, so the
TOTAL__OR__SPIN_UP channel is the total potential -- spin
handling is deferred; see DESIGN 5.7 / PSEUDOCODE 11.4).
Those three
producer sections will be revised to delegate to
kaleidoscope once this section stabilizes; that
revision is tracked as follow-up work, not performed
here.
The harvest contract reads one companion output beyond the
converged scfV. Because the producer-derived entry label
encodes the OLCAO (species, type) of each harvested site
(DESIGN 5.2.1), and those numbers are assigned by makeinput
during input preparation rather than known at manifest
time, the harvest stage needs them back at storage time.
Rather than add a new output, the producer co-opts the
datSkl.map file makeinput already writes (the sorted-dat
to skeleton atom-number mapping): it gains three columns
carrying each site's element, atom_species_id, and
atom_type_id. The harvester looks up the manifest
entry's atom_site (a skeleton-numbering index) and reads
its (species, type) straight off the map. The assigner
thus records its own verdict where the run lives, and the
harvest stage reads it back to mint the storage label
without re-parsing the run's Imago input. This keeps the
makeinput/harvest interface explicit: the run produces the
numbers, the producer only consumes them. The
producer-side change is tracked as TODO C87.
One seam this section once left open -- how the producer
supplies the dispatch Config (9.4) -- is now closed by
C100. The producer no longer forces a local executor: it
turns its dispatch choice into a flight Config through
kaleidoscope's shared generator and lets the driver
auto-select Local versus Parsl (9.4; DESIGN 6.2.11), so the
seed and database builds reach the scheduler rather than
running one calculation at a time on a login node.
Scheduler dispatch is the default -- a run with no site rc
file present is a configuration error, not a quiet local
fall-back -- and an in-process local run is the deliberate
opt-out that tests, a laptop, and the materialize pre-flight
request explicitly (DESIGN 6.2.11, decision 2). The wiring is
settled in DESIGN 6.2.11 -- the tiered site rc file, the
per-run CLI choices, the Config generator's home in
kaleidoscope, and the uniform-slice deferral of right-sizing --
and the code is C100 (done).
The producer's interaction with kaleidoscope is not always a
single fan-out. The predict-then-verify convergence search
(DESIGN 3.12) runs the k-point rungs as an adaptive climb: the
producer launches every material's seed rung (a small grid for
a confident prediction) at once, then collects rungs as they
land and, for each landing, reads that one material's energies
back, decides its next mesh, and launches it -- so a material
climbs on the instant its own rung lands and no material waits
on another, until each converges or hits a ceiling. This makes
the producer an iterative kaleidoscope client -- a running
stream of sends and collects -- rather than a single-batch one.
Its only addition to the flight layer is a small, additive one:
the dispatch's two phases, send-off and collect, exposed as
separately callable steps (DESIGN 6.2.3). It adds no dispatch-
side domain machinery -- each unit is still a flat, independent
CalcUnit, and the loop that reads energies and chooses the
next mesh lives entirely in the producer's Python. That is
precisely the division Principle 12 draws -- dependent, per-unit
iteration in client code, the dispatcher dumb -- so the climb
needs only the producer driving kaleidoscope one rung at a time.
The verify step of the guidance dataspace (section 10) is
correspondingly a confidence-gated grid-or-climb, not a single
fixed grid.
Most of the early questions are resolved and recorded
in the subsections above: the callable API offers both
a prepared-directory and a structure-and-options entry
(9.2); the run-reuse cache is a general kaleidoscope
mechanism with client-supplied key fields (9.6);
the ASE adapter stays a flat ase_imago.py (9.3); the
Atoms-to-StructureControl translation splits into
an ASE-free factory in structure_control.py plus
adapter-layer glue (9.3, 9.5); and kaleidoscope
dispatches through a pluggable wingbeat seam, defaulting
to imago.py but able to use the ASE adapter or future
adapters (9.4).
The workspace scheme (9.6) was deferred here and is
now resolved: the stable-id convention, the <calc>
tag format, and the status.toml schema are pinned in
DESIGN 6.2.4 (with each wingbeats/<id>[/<calc>]/ cache
directory carrying cache_key.toml, result.toml, and
status.toml).
What remains open:
-
Producer-section revisions (9.7) -- RESOLVED (C69). DESIGN 5.7, PSEUDOCODE 11.4, and ARCHITECTURE 8.5 have been rewritten so the producer delegates SCF to kaleidoscope: it materializes each reference structure (a thin local read, or a one-shot fetch of the pinned COD revision), builds a predict-then-verify flight, dispatches one flat batch, harvests each solid's converged grid point, and contributes that point back to the guidance dataspace (section 10). The bespoke per-solid SCF cache is gone; kaleidoscope's run-reuse cache (9.6) subsumes it. The matching code is TODO C74.
-
Cluster dispatch configuration (9.4) -- RESOLVED (DESIGN 6.2.11). The four decisions are settled:
- Site-config home. A dedicated, tiered
*rc.pyresource-control file (pure data, like every other*rc.py) -- a tiny required core (queues,worker_init, and account where the cluster demands one) with every performance and advanced knob optional and defaulted -- rather than a section ofimagorc. A separate tool,cluster_probe.py, reads the discoverable tiers off the scheduler and node hardware (sinfo,scontrol,lscpu) and emits a starter file, leaving the non-discoverable required core as blanks to complete. - Per-run choices. CLI flags (
--dispatch,--partition,--nodes,--walltime) defaulting from the rc file, with an optional resolved-config file written beside the run for a reproducible record. The command-line default isslurm-per-job: on a cluster the producer and seed reach the scheduler with no flags, and a run with no settings file is a configuration error rather than a quiet local fall-back.localis the deliberate opt-out (no settings file, noConfig), which the test suite and laptop sessions request explicitly. - Per-unit right-sizing. Deferred: both shapes give every unit a uniform slice now; right-sizing waits for a parallel imago and the section-11 cost predictor that would feed it.
- Generator sharing. The
Configgenerator lives inkaleidoscope, the dispatcher every flight already imports, so the producer and future clients share one copy.
Recorded at VISION Goals 4/6/7 and DESIGN 6.2.11; the producer change-over and the generator are TODO C100.
- Site-config home. A dedicated, tiered
-
Reclaiming a retired chain's workers -- DESIGNED, deferred (DESIGN 6.2.11). As the k-point climb retires materials, the pooled block goes on holding cores whose workers now have no work; handing those freed cores to the chains still climbing lets a late, expensive chain finish on several cores instead of one. It is distinct from the per-unit right-sizing above -- it predicts nothing, merely handing out what retirement has demonstrably freed -- and, like right-sizing, waits on a parallel imago. The producer keeps the core count itself (Parsl counts busy workers, not cores), so it needs no dispatch-side machinery; the design, the safety argument, and the three disciplines that keep the climb's wait-for-any code able to grow it are all in DESIGN 6.2.11.
The fifth prong (VISION Goal 5) accumulates a curated dataspace of converged calculations and a small predictor trained on it. Each entry is a datapoint: the chemistry-and-structure signature of a converged system, the resulting electronic-structure character (band gap, spin polarization), and the convergence settings (initially k-point density) that worked. The predictor is a two-stage k-nearest-neighbor regression: first it predicts the electronic character from chemistry, then it predicts the k-density from the electronic character. New calculations query the predictor, get a predicted operating point plus an uncertainty, and run a verification grid whose width is set by the uncertainty. Each successful flight appends back into the dataspace, so the predictor gets better as it accumulates evidence.
What the dataspace deliberately does not hold: metals. Every entry asserts that a k-density converged, and a metal cannot assert it -- its energy keeps moving as the mesh crosses the Fermi surface, which is why the climb stops at the first gapless mesh and settles there as a knowingly rough result (DESIGN 3.12.3). Recording that stopping point as a converged density would feed the predictor a claim nobody made, and in a young collection a lone metal can be the only member of its lattice family and so dominate every prediction for it. The harvest therefore skips metals (DESIGN 7.8). The initial potential database is unaffected and still takes their potential, which is exactly the rough starting point it wants (section 8).
This section covers on-disk layout, file format choice, the data-flow through kaleidoscope's flight builder, the harvest pipeline, and the module boundaries. Algorithmic details (feature-vector definition, k-NN metric, predictor stages, confidence measure, grid widening) live in DESIGN section 7.
The dataspace lives under a new top-level directory in
share/, parallel to share/atomicPDB/ (Goal 3) and
share/atomicBDB/ (curation manifest cache):
share/
historicalGuidanceDB/
entries/ Canonical entries. One TOML
file per entry. Two-level
partition by system_type:
crystalline/, amorphous/,
nanostructure/, molecular/.
File slug:
`<system_type>-<short_sha>.toml`.
staging/ Auto-harvested entries
awaiting curator review.
Mirrors `entries/`s system_type
partition. Once promoted,
files move from staging/ into
entries/; the staging area is
not consumed by the predictor.
superseded/ Entries retired at promotion
because the collection already
holds a record for that
structure under those settings
(DESIGN 7.8). Mirrors the
same system_type partition.
Retired, not deleted: the
record of what a re-run
produced stays recoverable,
and staging/ does not accrete
files that every later
promotion pass re-examines.
Not consumed by the predictor.
SCHEMA_VERSION Single line containing a bare
decimal integer (e.g. `1\n`).
Bumped on schema change.
Readers refuse files whose
`schema_version` mismatches.
elemental_groups.toml Element-group classification
table (alkali / alkali-earth /
halide / chalcogen / ... / H).
Loaded by the library at init
to compute composition
vectors. Checked-in data,
not code, per Principle 11.
Why partition only by system_type, not by gap or chemistry. The predictor operates on a continuous feature space (composition vector + lattice family + measured gap and spin-polarization); the only categorical axis is system_type, because the predictor's stage-2 relationship (gap+spin → k-density) is qualitatively different across system_types (a crystalline insulator's k-density depends on chemistry; an amorphous insulator's is set by cell size via the density convention; a molecular cluster is always Γ-only). Partitioning by system_type lets the predictor switch sub-models cleanly; partitioning further by gap or chemistry would force a continuous variable into bins prematurely.
One file per entry, not one big file. Easy append, easy git diff, no atomic-write contortions during a 500-entry parallel seed run, and a directory of 500-5000 small files stays tractable. The predictor loads all entries into an in-memory dataspace at init time -- TOML parsing of a few thousand small files takes well under a second on commodity hardware.
TOML, same rationale as ARCHITECTURE 8.2: small files,
human inspection a hard requirement, Python stdlib
tomllib for reads, hand-formatted emitter for writes.
The full schema, deterministic emitter contract, and
worked sketch live in DESIGN 7.2 / 7.3 / 7.5.
Architectural invariants:
schema_versionon every file. Lets the reader refuse unknown versions; lets future schema additions remain additive. Day-1 ships v1.- Signature block first. Every entry carries an
[entry.signature]block declaringsystem_type, the 13-dimensionalcomposition_vectorkeyed by element-group, and (for crystalline only) the 6-axis one-hotlattice_family. This is the predictor's feature input. - Measured-quantities block.
[entry.measured]carries the quantities harvested from the converged calculation:gap_ev,gap_kind(direct/indirect/none),spin_polarization,total_magnetization, and the convergedkpoint_density(the target the predictor learns to produce). - Context block.
[entry.context]carries parameters that influenced the measurement:basis(mb/fb/eb),functional, andscf_threshold. Future predictors may use these as additional regression features; v1's predictor conditions onbasisandfunctional(separate sub-models) and ignores the rest. - Provenance is required. Every entry carries the flight id, the source structure id, the imago commit, and a UTC timestamp. Non-negotiable per Principle 11.
new calculation (structure, options, system_type)
|
| The kaleidoscope flight-builder helper
| (DESIGN 6.2.8) asks the predictor:
| "given this structure, predict the converged
| k-density and tell me how confident you are."
v
guidance_db.predict(dataspace, query_signature,
basis, functional, kpoint_integration)
| (free function over the loaded Dataspace, DESIGN
| 7.4; query_signature carries system_type, and
| (basis, functional, kpoint_integration) selects
| the sub-model.)
| 1. Switch on system_type:
| - amorphous, nanostructure, molecular:
| consult the small dedicated sub-model
| (typically returns "use the density
| convention floor", i.e. ~Gamma-only for
| large cells), and return.
| - crystalline: run the two-stage regression
| below.
| 2. Stage 1 (chemistry -> electronic character):
| k-NN regression in composition+lattice
| feature space, predicting (gap, intensive
| magnetization = |M|/atom).
| 3. Stage 2 (electronic character -> k-density):
| k-NN regression in (gap, magnetization) space,
| predicting kpoint_density.
| 4. Confidence: variance over the k nearest
| neighbors at each stage, combined into one
| uncertainty score in [0.0, 1.0].
| 5. Return PredictionResult(
| predicted_kpoint_density, confidence,
| neighbor_entry_ids, predicted_gap,
| predicted_magnetization,
| is_under_trained,
| ).
v
The helper builds the verification grid using the
variance-aware widening function (DESIGN 7.7): narrow
when confidence is high, wide when low, wide-grid
fallback when `is_under_trained` is set.
|
| Kaleidoscope dispatches the grid (DESIGN 6.2).
v
flight harvest hook (10.5)
|
| 1. Pick the converged grid point per structure
| (smallest density where consecutive grid
| points' energy delta < threshold).
| 2. Read measured gap and total magnetization from
| the converged calc's result.toml (the kpd comes
| from the calc tag; cell facts from the structure).
| 3. Build a richly-populated GuidanceEntry and
| emit it to staging/<system_type>/.
v
Curator promotion (manual one-at-a-time, batch-promote
for trusted automated flights, or dry-run preview)
moves staging entries into entries/<system_type>/.
|
v
The promoted entry joins the in-memory dataspace on
the next library load and is visible to all future
predict() calls.
The predictor runs at flight-construction time -- well before any Imago run. The harvest runs at flight- completion time. Imago itself is unaware of the dataspace. All format awareness lives in the new helper module (10.6).
The predictor operates on a feature space that combines a chemistry vector with a lattice-family vector (crystalline only) and a coarse system-type partition:
- Composition vector (13-dim): atom-fraction weight
in each of 13 element-group buckets -- alkali,
alkali-earth, halide, chalcogen, pnictogen, group-IV,
group-III, transition-metal (lumped 3d/4d/5d),
lanthanide, actinide, metalloid, noble-gas, hydrogen.
The group classification table (
elemental_groups.toml, 10.1) is a checked-in data file; Principle 11 requires the chemistry knowledge to be auditable, not buried in code. - Lattice family (6-axis one-hot, crystalline only): cubic / hex / tet / ortho / mono / tri. Cheap-to-extract from the structure file via StructureControl. Mitigates the gap-vs-polymorph smoothness risk (DESIGN 7.6).
- System type partition (4-way): crystalline, amorphous, nanostructure, molecular. Hard switch for which sub-model runs.
The predictor is k-nearest-neighbor with inverse- distance weighting, in two stages for crystalline: stage 1 maps composition+lattice → (gap, intensive magnetization); stage 2 maps (gap, magnetization) → k-density. The split exploits transferability: stage 1 is chemistry-heavy and needs broad chemistry coverage; stage 2 is physics-heavy and is roughly material-independent. Each stage carries its own confidence (variance of the k neighbors); the two combine into one uncertainty for the verification grid.
For non-crystalline system types (amorphous, nanostructure, molecular), the predictor returns a canonical result driven by the density convention plus any system_type-specific corrections; chemistry plays little role. The implementation is correspondingly simpler.
Full algorithm (k value, distance metric, normalization, confidence formula, weighting under sparse data) is in DESIGN 7.6. The exact bootstrap behavior when the dataspace is too thin for stage 1 to be trustable is in DESIGN 7.9.
The dataspace is a build product like the initial- potential database (8.5), but with one key difference: every successful flight is a potential contributor, not just a hand-curated reference set. This makes the harvest hook central machinery:
src/scripts/
guidance_db.py Library: read,
compute composition
vectors, run the
predictor, emit
entries. No
orchestration.
guidance_harvest.py Producer-side
helper: given a
finished flight,
examine each
structure's
verification grid,
read measured
quantities from
each result.toml,
write staged
entries.
guidance_promote.py Curator helper:
review and promote
staging entries
into the canonical
entries directory.
guidance_harvest.pyis invoked at flight-completion time (either as a post-step the flight driver calls, or as a standalone CLI run after the fact). It reads the flight's workspace, identifies each verification grid, picks the converged point per structure, reads the measured electronic-structure quantities from the converged calc's result.toml, and emits one staged TOML per structure into the appropriateshare/historicalGuidanceDB/staging/<system_type>/subdirectory.guidance_promote.pyis the curator's tool. Four modes: interactive one-at-a-time review (default); batch--auto-promotefor trusted automated flights meeting an objective acceptance rule (the converged k-density landed in the middle 60% of the verification grid AND the top three grid points' total-energy spread -- max minus min, read from the staged entry'sgrid_energiesarray -- is below threshold; a spread, not a variance, so the quantity compares like-for-like against the linear threshold, DESIGN 7.8);--allto promote the whole staging directory after manual review;--dry-runpreview. The auto-promotion rule lets a 500-entry seed flight promote ~80% of entries unattended, with the curator reviewing only the ~20% outliers. Promotion is also where a re-run of an already-promoted solid is caught, because it is the only stage that sees both the incoming entry and the promoted corpus (DESIGN 7.8). The harvest cannot: it writes one entry per converged non-metal solid and has no view of what a curator accepted months ago. The rule the check enforces is simply that the collection holds one record per structure per settings combination -- the duplicate-structure guard, and nothing more. Every mode applies it,--allincluded, because keeping one structure out of the predictor's neighbor set twice is a correctness guard rather than a quality judgment, and--allwaives only the latter. What the check does not do is adjudicate: an arriving record for an occupied slot is reported with both meshes and both build identities, the existing record is left exactly as the curator accepted it, and the newcomer is retired tosuperseded/. Interactive review adds a REPLACE choice for the occupied case, and that is the only way an entry ever leavesentries/-- a curator's verb, never an unattended mode's (VISION 16).- Both databases grow incrementally and in place (the
initial-potential DB does too, DESIGN 5.2.3); the
difference is what gates the growth. The
initial-potential DB grows only from a curated manifest
of reference solids, whereas the dataspace grows
monotonically from any successful flight, with no
curation manifest in the loop. Old entries are not
deleted on schema bumps; a migration tool
(
guidance_migrate.py) rewrites them in place.
The seed flight (TODO C75) is what populates the initial dataspace. It is a one-time stratified sweep across element-group pairs and common stoichiometry patterns (~150-250 calculations covering the chemistry surface representatively rather than at random) feeding the auto-promotion rule above. After the seed lands, ongoing flights (the C48.3 producer, future characterization runs, etc.) contribute additional entries as they finish.
Python:
src/scripts/guidance_db.py: new library. TOML reader + validation; element-group classifier (loadselemental_groups.toml); composition-vector and lattice-family computation from a StructureControl; the in-memory dataspace; the two-stage k-NN predictor; the deterministic hand-formatted emitter. Imports onlytomllib,math, and the existingstructure_control.py. Module-level docstring describes its role as the library half of the library/producer/consumer split (10.5).src/scripts/kaleidoscope/flight-builder helper (DESIGN 6.2.8): consumes aPredictionResultand builds aFlightofCalcUnits laid out per the tag convention of DESIGN 6.2.4. This is the option-axis half of the builder split per VISION Principle 12; the structure-axis half remains domain-specific and lives instructure_control/ acquisition.src/scripts/guidance_harvest.py: new producer helper.src/scripts/guidance_promote.py: new curator helper.src/scripts/guidance_migrate.py: future schema- migration tool. Not in day-1 scope.
Fortran:
- A small extension to imago.py's harvest path so that
gap_ev,gap_kind, andtotal_magnetizationare surfaced (alongside the total energy already present). These quantities are computable from existing SCF output (the eigenvalue spectrum); the change exposes them in the iteration data so any plain SCF run yields them (DESIGN 6.1). Tracked as TODO C76 (under Phase K).
The library / producer / consumer split mirrors DESIGN 5: the library knows the format and runs the predictor; the producer (harvest + promote) writes entries; the consumers (kaleidoscope builder and any future client) read entries.
- DESIGN 5 (initial potential database): the
guidance dataspace stores convergence settings +
electronic-structure character; the potential DB
stores converged potentials. They share the
library/producer/consumer discipline (Principle 11)
and the same
share/shape (per-element or per-signature TOML), but they are independent artifacts with independent lifetimes and no cross-references. A guidance entry never names a potential-DB entry, and vice versa; the two share only the curation discipline, not their contents. Considered and rejected in DESIGN 7.10 ("Closed by decision"): the two serve different audiences and update cadences, and entangling their schemas would couple their lifetimes unnecessarily. - DESIGN 6 (kaleidoscope): kaleidoscope is the
dispatch layer that runs the verification grid the
predictor produces. The dependency goes one way:
the kaleidoscope flight-builder helper (6.2.8)
reads the dataspace; the dataspace does not depend
on kaleidoscope. Kaleidoscope still works without
the dataspace -- callers can construct
CalcUnits directly without going through the helper. - C48.3 producer (initial-potential database build): the first major consumer. Its workflow under the new prong: for each curated reference solid in the potential-DB manifest, the producer queries the guidance predictor (system_type = crystalline; structure file already known), receives a predicted k-density and confidence, lets the kaleidoscope helper build a verification sub- grid around the prediction, dispatches the sub-grid, and harvests the converged potential from the run dir whose grid point converged. Without the dataspace seeded, this would run as a wide-grid sweep per solid; with the seed in place, the sub- grid shrinks to 3-5 points per solid (predict-then- verify acceleration). The seed flight (C75) therefore directly accelerates Goal 3.
- File naming uniqueness under parallel harvest. Two flights finishing nearly simultaneously could collide on a slug. Plan: include a short SHA over (flight_id, source_structure, generated_at) in the filename; the flight+structure pair is unique by construction, so collisions only occur if two threads in the same flight harvest the same structure at the same instant. Detailed slug algorithm in DESIGN 7.5.
- Polytype confusion within the predictor. Two polymorphs of the same compound at the same composition vector but different lattice family currently produce two entries with the same chemistry features but different lattice features. The k-NN distance metric weights these features; whether the default weights separate polytypes cleanly enough is an empirical question that won't be settled until the seed flight lands. Open knob: the relative weight of composition vs lattice in the stage-1 distance metric. Reasonable default in DESIGN 7.6; calibration after seed.
- Spin polarization unit and interpretation.
Recorded as
total_magnetization(Bohr magnetons per formula unit); but for antiferromagnets the net magnetization can be zero while the magnetic ordering is non-trivial. A future schema bump may add a separatelocal_moment_per_atomfield. Day-1 we document the limitation and proceed. - Functional / basis as sub-model dimensions vs features. The predictor conditions on basis and functional by running separate sub-models per (basis, functional) combination. Alternative: treat them as additional regression features. The split approach is cleaner (no spurious cross-functional interpolation) but proliferates sub-models. Open; default approach in DESIGN 7.6 is sub-models per combination, with fallback to the most-similar combination when the queried one is empty.
VISION Goal 6 adds a second curated dataspace that records, for every imago run, what it cost to compute -- peak memory, disk footprint, and walltime -- as a function of the problem size and the parallel configuration the run used. It is a deliberate sibling of the historical-guidance dataspace (section 10), not an extension of it: the two share the library / producer / consumer discipline (Principle 11), the staging-then-promote curation flow, the schema-versioning-plus-migration pattern, and the registry-validated key discipline -- but they are independent artifacts with independent schemas, predictors, and lifetimes.
The reason they must stay separate is portability. A
section-10 entry ("MgO converges at k-density 60") is
hardware-independent: it is equally true on a laptop and a
national supercomputer, so that dataspace accumulates
globally and never goes stale when hardware changes. A
resource entry ("this run used 42 GB and 3.1 h") is the
opposite -- it is meaningful only on the machine that
produced it. Folding cost fields into section-10 entries
would contaminate the one artifact whose value depends on
being portable. The resource dataspace is therefore
partitioned by a hardware fingerprint (11.2), the way
section 10 partitions by system_type.
The artifact lives under its own root in share/ (proposed
share/resourceGuidance/, final name in DESIGN 8),
mirroring section 10's entries/ + staging/ split and
its SCHEMA_VERSION marker. Entries are partitioned one
directory level down by hardware fingerprint rather than by
system_type:
<root>/
SCHEMA_VERSION bare-integer marker.
hardware_registry.toml known fingerprints + the
probed attributes that
define each (analogous to
section 10's elemental_groups).
entries/<fingerprint>/*.toml promoted observations.
staging/<fingerprint>/*.toml harvested, awaiting a
curator.
The hardware fingerprint is a short stable key (e.g.
<cpu-model>-<cores-per-node>-<mem-per-node>, exact recipe
in DESIGN 8) naming the node type a run executed on. It is
the coarse partition within which cost is comparable. The
finer parallel-configuration knobs (rank/thread counts,
pinning) and the build configuration (11.3) instead live
inside each observation and feed the predictor as features,
not partitions: hardware is partitioned because cost is never
comparable across machines, but parallel-config and build
choices are kept comparable on purpose -- the whole point of
recording them is to learn which configuration, and which
build, is cheapest.
The atomic unit is one execution observation -- a single run under a single configuration -- never collapsed to a per-system or per-fingerprint summary (VISION Goal 6). Two runs of the same structure under different rank counts are two observations. This granularity is what lets the same artifact answer "how much will this configuration need" (provisioning, near-term) and, later, "which configuration is most efficient" (optimization) and "how does cost scale with size" (scaling studies), with no schema change.
Each observation TOML carries four content blocks plus provenance (field-level schema in DESIGN 8):
- Size signature -- the cost-driving dimensions of the problem: atom count, electron count (with the core/valence split, since the orthogonalized secular dimension is the real driver), basis-function count, wave-function representation (1-component Schrodinger vs 4-component Dirac -- the 4-spinor structure multiplies the secular dimension), k-point count, and spin channels. Derived from the makeinput inputs and the structure; known before the run.
- Execution configuration -- how the run was launched:
node count, cores per node, total cores, MPI rank count,
OpenMP threads per rank, and process/thread binding (core-
or socket-pinning). Captured as an extensible,
registry-validated key-value table: a checked-in
EXECUTION_KNOB_REGISTRYenumerates the recognized knobs, and a new knob (a GPU count, a NUMA policy) is added by extending the registry and bumping the schema version -- never by silently redefining fields. Known at dispatch time from the flight's Parsl provider / launch spec. - Build configuration -- the toolchain the binary was
compiled with, recorded in two layers so it serves both
the accumulating artifact and a targeted experiment. The
coarse layer is a registry-validated set of normalized
knobs that act as predictor features and the axis the
artifact groups by: compiler family and major version,
optimization level and a coarse arch/SIMD tag, and for each
key library (HDF5, ScaLAPACK, BLAS/LAPACK, MPI) its
implementation, major version, and cost-relevant variant
(HDF5 parallel vs serial, BLAS threaded vs sequential).
Like the execution configuration this is an extensible
BUILD_KNOB_REGISTRY, but deliberately bucketed (an optimization level, not a flag string; a major version, not a patch) so a build is a comparable feature, not a fragmenting one. The fidelity layer is the complete compile string and full library build details, stored verbatim as provenance (never a predictor feature, just a string -- cheap to keep), so nothing is lost. A later study of one specific flag then has two routes: mine its state out of the verbatim string post-hoc, or promote the flag to a named knob in the registry so it becomes a first-class comparable feature for that study -- without coarsening the global artifact. Such a study is itself just a flight whose varied axis is the flag (on/off) across a limited grid of systems (the DESIGN 6 sweep). Known at compile time (11.4). - Measured resources -- what the run actually used: peak
resident memory, disk footprint (output plus scratch
high-water mark), and walltime, with optional per-phase
timings (setup, SCF iteration, eigensolve,
post-processing). Like section 10's
METRIC_REGISTRY, the metric set is itself a registry so new measurements append cleanly. Captured after the run (11.4).
The producer is again the flight layer; the consumer is the provisioner. Capture draws from four sources, joined at harvest:
- Dispatch-time (known before the run): the size signature (from makeinput + structure) and the execution configuration (from the kaleidoscope flight's Parsl provider and launch spec). The wingbeat records these into the run directory alongside its inputs.
- Build manifest (compile time): the CMake build emits a
build_info.toml-- compiler and full flag string, plus the detected versions and variants of HDF5, ScaLAPACK, BLAS/LAPACK, and MPI -- installed beside the binary, and imago can echo its compiled-in configuration. The harvest reads this for both layers of the build block (11.3): the coarse knobs and the verbatim compile string. - Scheduler accounting (after the run): SLURM
sacctsuppliesMaxRSS(peak memory), disk read/write high-water marks, andElapsed(walltime); a/usr/bin/time -vwrapper is the fallback off-scheduler. - Imago self-report (optional, after the run): the
Fortran side may write its own peak allocation and
per-phase wall timings into
result.toml, which is more precise than scheduler granularity for the eigensolve phase specifically.
A dedicated harvest helper (resource_harvest.py, parallel
to section 10's harvest) walks the finished flight's
workspace, assembles one observation per run dir from these
sources, and writes it to staging/<fingerprint>/. A
curator promotes with the same staging-then-promote
discipline as section 10. The two harvests are separate
tools reading the same workspace; a single completed run
can feed both a section-10 convergence observation and a
section-11 cost observation.
Unlike section 10, failed runs are not always discarded: an out-of-memory or walltime-exceeded run is positive evidence that a configuration is insufficient at that size, which is exactly the signal a provisioner needs. How much of that negative data to retain -- and how to mark a censored "needed at least X" so the predictor treats it differently from a measured "used exactly X" -- is an open question (11.8).
Within a hardware fingerprint, cost is a smooth, physics-grounded function of the size signature and the execution configuration. This argues for a physics-informed regressor rather than the pure k-NN of section 10: memory grows roughly as the square of the secular dimension and the eigensolve as its cube, so the model can fit scaling exponents (a power law in the secular dimension, scaled by the parallel configuration) instead of interpolating raw neighbors. A k-NN fallback is reasonable when a fingerprint is thinly populated. The predictor choice is deferred to DESIGN 8 and is explicitly not required for the artifact to begin accumulating data.
The near-term consumer is provisioning (VISION Goal 6): given a chosen parallel configuration for a new run, the predictor estimates the memory, disk, and walltime the run will need, and the flight layer turns those into SLURM resource requests with a safety margin. This directly serves Principle 6 (cost discipline at SCF setup) and replaces today's hand-guessed requests, whose failure modes are wasted allocation (over-request) or killed jobs (under-request).
Because every observation stores its full configuration, the same artifact later supports three forward-looking consumers with no schema change: configuration optimization (compare predicted cost across candidate configurations to recommend the most efficient), build comparison (compare cost across compiler flags or library builds on one fingerprint -- the targeted flag study of 11.3), and scaling studies (fit and report cost-vs-size curves). These are not near-term, but the granular-observation design (11.2) is chosen precisely so they need no migration.
src/scripts/resource_db.py(new): the library -- load / validate / hand-formatted emit, the hardware fingerprint plus registries, and the predictor. Built on the same patterns asguidance_db.py(section 10 / C70).src/scripts/resource_harvest.py(new): the harvest producer (11.4).- A provisioning consumer in the flight layer (the kaleidoscope flight-builder helper, DESIGN 6.2.8, or a thin sibling) that reads the predictor and annotates the Parsl provider's resource request. As with section 10, the dependency is one-way: the flight layer reads the dataspace; the dataspace does not depend on the flight layer.
src/scripts/resource_migrate.py(future): the schema migration tool, mirroring the plannedguidance_migrate. Not day-1 scope.- CMake build-system hook: emit
build_info.toml(compiler and full flag string, detected library versions/variants) at configure/install time, so the build block is captured without hand-entry (11.4). - Wingbeat /
imago.pycapture hooks (11.4) to record the dispatch-time config and scrape scheduler accounting into the run directory.
- Section 10 (historical guidance): orthogonal axes of the same run. Section 10 answers "what operating point is accurate" (portable, chemistry-keyed); section 11 answers "what will it cost" (hardware-keyed). One completed run can produce one observation in each. They share library scaffolding and curation discipline only, never contents or cross-references -- the same boundary section 10 draws against DESIGN 5 (10.7).
- DESIGN 6 (kaleidoscope): the flight layer is the producer (it captures config + accounting) and the home of the provisioning consumer. The dataspace does not depend on kaleidoscope; flights run fine without it (callers request resources manually).
- DESIGN 5 (initial-potential build) / C48.3: a heavy batch consumer. Once seeded, the potential-DB build can size its many per-solid SCF jobs from predicted cost rather than a one-size-fits-all request, cutting both queue rejection and wasted allocation across a large flight.
- Hardware fingerprint granularity. Too coarse (just CPU model) lumps differently-provisioned partitions together; too fine (every BIOS/microcode revision) fragments the data so no fingerprint accumulates enough observations to predict from. The right granularity is empirical; DESIGN 8 proposes a default recipe and a fallback when a fingerprint is under-populated.
- Reliable peak-memory capture.
sacct MaxRSSis sampled and can miss short-lived spikes;/usr/bin/timeis per-process, not per-job; the Fortran self-report is precise but only for what imago itself allocates (not MPI buffers or libraries). Which source is authoritative, and how to reconcile them, is open. - Node sharing and contention noise. A run that shares a node with other jobs sees inflated, noisy walltime and memory pressure. Whether to record an exclusivity flag and weight or filter shared-node observations is open.
- Censored (failed-run) data. As in 11.4, an OOM or timeout is a bound, not a point measurement. Whether and how the schema and predictor represent censored observations is open.
- Any portable normalization. Walltime is not portable, but core-hours or a hardware-normalized cost index might transfer partially across fingerprints. Whether a normalized cost is worth recording (to bootstrap a new fingerprint from related ones) is open.
- Which build knobs deserve coarse-layer promotion. The
two-layer build record (11.3) bounds this -- the full
compile string is always kept, so nothing is lost -- but
which knobs are stable and cost-relevant enough to live
in the coarse
BUILD_KNOB_REGISTRY(and thus drive the predictor by default) is a judgment call DESIGN 8 must seed and the curator revisits as evidence accumulates. - Build effects on numerics, not just cost. A build
block primarily conditions cost here, but the same choices
(fast-math, a different BLAS summation order) can perturb
low-order digits of the result. That is a reproducibility
concern straddling section 10, where today only
imago_commitkeys the binary's identity. Whether the build block should also be referenced from the convergence side -- against the no-cross-reference boundary (11.7) -- is an open question to settle in DESIGN, not assume now.
Serves VISION Goal 7 (parallelize across scales and devices) and VISION Principle 15 (attribution). Consumed by DESIGN 10.
The engine currently writes to its log unconditionally. Every
write (20,...) in the Fortran fires on every run, so the only
way to reduce output is to stop calling the code that produces
it. That is workable while the log is small, but it does not
survive high-throughput use: a flight of a hundred units
multiplies every decorative line by a hundred, and the debugging
and parallelization work of Goal 7 will want diagnostics that
are useful when requested and intolerable by default.
The facility has two constituencies, which is why it is established here rather than inside the design that first needs it. Principle 15 requires the program to state its own citation, and that block must be suppressible or it becomes noise a high-throughput user learns to ignore. Goal 7 requires diagnostics that would drown an ordinary run. Both are asking for the same thing: output that is present when wanted and absent when not.
This section establishes a single facility for that control. It is deliberately minimal today -- one category is defined and one call site uses it -- because the interesting decision, which categories exist and what each includes, belongs with the debugging and optimization work rather than ahead of it. What is settled now is the mechanism, chosen so that deciding the categories later costs nothing.
Verboseness is requested through the environment variable
IMAGO_VERBOSENESS, whose value is a comma-separated list of
category names:
IMAGO_VERBOSENESS=normal # the default when unset
IMAGO_VERBOSENESS=none
IMAGO_VERBOSENESS=banner,timing # order is irrelevant
Internally the selection is one integer whose bits correspond to
categories, tested with the btest intrinsic. That
representation is private. Names are the public contract and
bit positions are an implementation detail, and the
distinction is the point of the design rather than a stylistic
preference.
A numeric form (IMAGO_VERBOSENESS=25) was considered and
deliberately not adopted. Accepting a number publishes the bit
assignment the moment anyone writes one into a job script, and
this project records every invocation into a command file, so
such settings persist in the permanent record. Renumbering
afterward silently changes what old scripts do. Names have no
such cost: they can be reordered freely, and an unrecognized one
can be reported where a stray bit cannot. A numeric form can be
added later without breaking anything if it proves wanted; it
cannot be withdrawn.
Names are also the readable choice. banner,integrals states
what it asks for; 25 requires the reader to decompose it into
bits and consult a table, which is the same objection this
project makes to naming a variable np instead of nuc_pot.
An O_Verboseness module holds the parsed mask and the name
table that defines it. The table -- one entry pairing each
category name with its bit position -- is the single source of
truth, which is what makes adding a category a one-line change
rather than a structural one.
initVerboseness reads the environment variable once, splits it
on commas, matches each token against the table, and sets the
corresponding bits. It is called from parseCommandLine
immediately after the log file is opened and before the first
operation timestamp. That slot is forced: the log unit does not
exist earlier, and the first timestamp is already output that a
category should be able to govern.
isVerbose answers whether one category is active. Callers
pass a named parameter and never a literal bit number. If any
call site hardcodes a position, renumbering stops being free and
the property that makes the design safe is lost.
Three behaviours are fixed:
- Unset means
normal, not silent. A fresh installation with no environment configured must produce a useful log; a mute program would look broken. - Unknown names warn and continue. A typo must not kill a long cluster job, but it must not silently do nothing either, so the parser reports what it did not recognize to the log.
noneis explicit. Silence is a thing a caller asks for by name, not a thing it falls into by leaving a variable empty.
One category is defined: banner, governing the identity and
citation block specified in DESIGN 10. normal includes it.
The categories that will matter for the debugging and parallelization work -- SCF iteration detail, integral and setup diagnostics, resource reporting, developer trace -- are deliberately not enumerated here. Naming them before knowing what those campaigns need would fix a vocabulary against guesswork. Each is added later by appending one row to the name table, declaring one parameter, and gating the relevant writes.
Retrofitting the existing unconditional writes across the Fortran is a separate and much larger task. It is not required by DESIGN 10 and is not begun here.
- O_Verboseness (
verboseness.f90): the mask, the category name table,initVerboseness,isVerbose, and the category parameters. Depends on nothing butO_Kinds, so any module may use it without creating a cycle. - O_Banner (
banner.f90): the identity block and the citation text of DESIGN 10. UsesO_Verbosenessand nothing else. - O_MethodCitations (
methodCitations.f90): the method-citation registry of DESIGN 10.5 and its predicates. UsesO_KPoints, and later whatever further engine state a new predicate reads.
The last two were one module in the first draft of this
section, and that arrangement cannot be compiled. The identity
block prints from parseCommandLine, so O_CommandLine must
use it; the registry's predicates read kPointIntgCode, so it
must use O_KPoints; and kpoints.f90 already uses
O_CommandLine. A single module holding both halves closes the
cycle O_CommandLine -> O_Banner -> O_KPoints -> O_CommandLine, which no build order resolves.
Splitting them costs nothing, because the seam is one DESIGN 10.2 has already drawn. That section separates the two blocks by when their information becomes available, and the same fact governs what each may depend on: a block that prints before the work begins can depend on nothing the work produces, and a block that prints after it is finished may depend on all of it. PSEUDOCODE 17.1 records the reasoning and the alternative that was rejected.
Note that O_TimeStamps declares a banner variable
(timeStamps.f90:24) that nothing reads -- it is dead, and
predates this work. It should be removed rather than confused
with O_Banner.
Settled by the security campaign of 2026-09-07 and 2026-09-22,
tracked in dev/SECURITY_TODO.md -- a ledger of findings rather
than a level of the chain. No VISION goal or principle drives
this section: it was forced by an external vulnerability report
rather than chosen. It belongs in ARCHITECTURE anyway, because
what the campaign settled constrains module boundaries and call
shapes, and a ledger entry is read once while this has to govern
code nobody has written yet.
Imago reads files it did not write. A skeleton comes from a
collaborator, a structure comes from a public database, a
manifest is hand-edited by a curator, a scratch directory is
shared on a cluster. The engine also launches other programs:
scripts drive makeinput.py, imago.py and the Grasp2K
binaries as child processes. Those two facts together are the
whole of this section's subject. Everything below exists
because a campaign against an external vulnerability report
(SEC-001 through SEC-010, two rounds, 2026-09-07 and
2026-09-22) found real defects at exactly their intersection.
The rules are stated as architecture rather than as advice because they constrain module boundaries and call shapes, not just style, and because three of them were learned by getting them wrong first.
No value that arrives from a file, a database, a prompt or an
environment variable may be evaluated. eval, exec and
pickle.load are excluded from this codebase for anything with
an external origin, without exception and regardless of how
constrained the namespace looks. Four of the ten findings in
the first round were eval or exec on file content, and one
of them -- viewCell.py executing a BZ.<n> geometry file as
Python -- meant that any such file could run anything.
TOML is the project's data format for this reason and not
merely by taste. tomllib yields dicts, lists, strings,
numbers, booleans and dates, and has no way to express a call.
The contrast worth holding in mind is pickle.load, which
takes the same binary handle and is indistinguishable at the
call site from tomllib.load to anyone skimming, and which
executes arbitrary code by design with no injection required.
A reader checking this codebase for that distinction is
checking the right thing.
Where a value genuinely must be computed from text -- a lattice equation in a SYBD file, a k-point coordinate -- it goes through a purpose-built evaluator that understands only the grammar it is meant to accept. That is a parser, not a sandbox, and the difference is that a parser cannot be escaped from because it has nowhere to escape to.
Every call that starts a child process in this codebase takes the same form, and a new one that does not is wrong:
- The command is a list, never a string. The elements
become
argvthroughexecve, so no shell parses them and metacharacters carry no meaning. shell=Falseis written out at the call even though it is the default. The absence of a shell is the reason the arguments need no quoting, and a reader checking that should not have to recall a default.- The program name comes from a closed set of constants
declared in the module, checked before it is joined to a
directory.
graspElems.py'sGRASP_PROGRAMSis the model. This also closes path traversal: an unchecked name joined to a binary directory can escape it with../. - Data for the child goes over its standard input as a
file, not as an argument to the call that starts it. Writing
the input to a file first and passing the open file as
stdin=keeps the payload out of the call entirely, and it matches how these programs are driven from a shell anyway.
This is the rule that took two rounds to learn and is the one most easily got wrong, because the wrong version looks identical to the right one at a glance.
A test constrains what may proceed:
if not NUMERIC.fullmatch(str(value)):
raise ValueError(...)
command = [..., *values] # the ORIGINAL values
A rebuild changes what proceeds:
def numeric_argument(value):
text = str(value)
if not NUMERIC.fullmatch(text):
raise ValueError(...)
try:
return str(int(text))
except ValueError:
return repr(float(text))
command = [..., *[numeric_argument(v) for v in values]]
Both refuse the same inputs. Only the second guarantees the character content of what reaches the child, because the word it passes on was constructed by Python's own formatter and is not the word that arrived. A string reaching the rebuild cannot contribute a single character to the command.
The property to want is this: the safety of the call should be a local property of the function that makes it, not an inherited property of a caller several files away. A test leaves the guarantee upstream, where a later edit can invalidate it silently. A rebuild brings it home.
Where the value is genuinely free text and cannot be rebuilt
from a type -- a curator's description, an answer script for an
interactive program -- validate it at the single point it
enters the program, hold it to a positive shape where one
exists, and refuse rather than escape. expand_manifest.py's
clean_answer is the model.
The rule with the largest blast radius, learned from being wrong in the ledger.
makegroups.py does not import initial_potential_db. That
was checked, correctly, and then used to conclude that no path
existed between the two -- which was false. makegroups.py
exports loen_input_values, and build_initial_potentials.py
imports both modules and calls that helper with data that
originated in the database. The helper's return value then
flows, in makegroups.py's own other callers, to a
subprocess.run. No single run traverses the join. A
whole-program analysis traverses it easily, because it merges a
function's parameters and return value across every call site.
Three consequences for how modules are laid out and reasoned about:
- Reachability is not an import question. Before asserting that module A cannot reach module B, check who calls INTO A, not only what A imports. The import list answers a narrower question than it appears to.
- A helper exported to another module joins their data
flows. That is usually fine and often the point -- sharing
loen_input_valuesis exactly what keeps the producer's LOEN block and the grouping flow comparable. But the sharing is a seam, and it belongs in the seam inventory that any DESIGN or PSEUDOCODE section touching existing code has to write. - Put the guard between the shared helper and the sink. A check upstream of the helper protects one caller. A rebuild downstream of it, immediately before the call that consumes the value, protects every caller including ones written later.
src/scripts/ leaf utilities are worked inline against these
rules. A finding inside the compiled engine -- anything under
src/imago/, src/atomSCF/ and the other subprograms -- goes
down the document chain normally, because a change there is
engine algorithm work and the gate exists to protect exactly
that.
Two cautions, both earned.
Do not deform the code toward a scanner. These rules are written because each makes the program better to read and harder to misuse. None of them is here to satisfy a tool. Where a construct cannot be removed without harming the program -- and a script whose job is to launch another program must contain a call that starts one -- the right response is to state that plainly in the ledger and stop.
Do not predict what a scanner will report. This project has now made one confident prediction in each direction and been wrong both times: that removing a pattern would produce a clean report, and that an irremovable pattern would be reported forever. Neither prediction was worth making. Judge a change by whether it closes a real weakness and whether it leaves the code clearer, and let the report say what it says.