diff --git a/AGENTS.md b/AGENTS.md index 16cdb97..2f3b32f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -78,15 +78,16 @@ record keeps the convention it was published with, so a consumer must read the base from the record rather than assume it. The ladder is built from the k-distances at which `ceil(|b_i| / k_distance)` -changes on any axis — that is, from `|b_i| / n`. The enumeration cap on `n` is -applied per axis, so axes with different `|b_i|` exhaust their breakpoints at -different k-distances: the ladder is a complete set of meshes only for -`k_distance >= max(|b_i|) / n_max`, and below that it silently skips reachable -meshes. Any recomputation or published `k_index` column must state the cap it -used. - -This convention is shared with `goldilocks-core`. Changing it invalidates every -`k_index` value already recorded or trained on. +changes on any axis — that is, from `|b_i| / n`. Enumeration stops at a +resolution floor `min_k_distance` (default `0.03` Å⁻¹, on the solid-state 2π +lengths — the AiiDA-QuantumESPRESSO convention). The floor is the same for every +axis, so all axes stop together and the ladder is a complete, gap-free set of +meshes over the whole `[min_k_distance, ∞)` range. Any recomputation or +published `k_index` column must state the floor it used; a `k_index` from a +different floor is comparable only where the two ranges overlap. + +This convention must stay in step with `goldilocks-core` and `goldilocks-ml`. +Changing it invalidates every `k_index` value already recorded or trained on. ## Commands diff --git a/campaigns/qe/kpoints/README.md b/campaigns/qe/kpoints/README.md index e51561d..ab66774 100644 --- a/campaigns/qe/kpoints/README.md +++ b/campaigns/qe/kpoints/README.md @@ -9,7 +9,7 @@ PseudoDojo PBEsol campaign and the SSSP comparison analysis. ```text kpoints/ campaign.yaml human-readable campaign configuration - scripts/ initial submission, extension, and monitoring + scripts/ extension and monitoring notebooks/ curated analysis and visualisation results/ snapshot, summary, and provenance manifest ``` @@ -49,8 +49,6 @@ uv run --extra aiida --extra kmesh python campaigns/qe/kpoints/scripts/monitor.p ``` The monitor calls `extend.py` in a fresh process. Stop it with `Ctrl-C`. -`submit_initial.py` is for the original seed campaign and additionally needs -the historical convergence summary files; run `--help` for its full interface. ## Convergence definition diff --git a/campaigns/qe/kpoints/scripts/extend.py b/campaigns/qe/kpoints/scripts/extend.py index 556760a..17dba5a 100644 --- a/campaigns/qe/kpoints/scripts/extend.py +++ b/campaigns/qe/kpoints/scripts/extend.py @@ -16,7 +16,8 @@ from goldilocks_data.aiida.submit import existing_kindices_by_source, submit_jobs from goldilocks_data.codes import DftCode from goldilocks_data.intents import CalculationIntent -from goldilocks_data.sweeps import AiidaJobSpec, KindexExtension, kindex_points, plan_well_not_ultra_extensions +from goldilocks_data.kmesh import kindex_points +from goldilocks_data.sweeps import AiidaJobSpec, KindexExtension, plan_well_not_ultra_extensions TASK_ROOT = Path(__file__).resolve().parents[1] DEFAULT_SNAPSHOT_DIR = TASK_ROOT / "results" diff --git a/campaigns/qe/kpoints/scripts/submit_initial.py b/campaigns/qe/kpoints/scripts/submit_initial.py deleted file mode 100644 index 6dd36f7..0000000 --- a/campaigns/qe/kpoints/scripts/submit_initial.py +++ /dev/null @@ -1,279 +0,0 @@ -from __future__ import annotations - -import argparse -from dataclasses import dataclass -from pathlib import Path - -import pandas as pd -from aiida import load_profile -from aiida.orm import StructureData -from ase.io import read -from pymatgen.io.ase import AseAtomsAdaptor - -from goldilocks_data.aiida import AiidaScfConfig -from goldilocks_data.aiida.cleanup import cleanup_finished -from goldilocks_data.aiida.registry import failed_source_ids -from goldilocks_data.aiida.submit import existing_kindices_by_source, submit_jobs -from goldilocks_data.codes import DftCode -from goldilocks_data.intents import CalculationIntent -from goldilocks_data.sweeps import AiidaJobSpec -from goldilocks_data.sweeps.kindex import kindex_points -from goldilocks_data.sweeps.kmesh import generate_candidate_k_distances, k_distance_to_mesh - -SUMMARY_NOTE_COLUMN = "Convergence Notes" -FULL_NOTE_COLUMN = "Convergence Notes (modified)" -DEFAULT_SOURCE_DB = "Goldilocks_scf_K_mesh_convergence_nospin" - - -@dataclass(frozen=True, slots=True) -class LocalJobInput: - source_db_id: str - kindex_max: int - - -def main() -> int: - args = parse_args() - load_profile() - - config = AiidaScfConfig( - code_label=args.code_label, - pseudo_family_label=args.pseudo_family_label, - group_label=args.group_label, - degauss_ry=args.degauss_ry, - num_machines=args.num_machines, - num_mpiprocs_per_machine=args.num_mpiprocs_per_machine, - max_wallclock_seconds=args.max_wallclock_seconds, - ) - - active = active_workchain_count(config.group_label) - print(f"active workchains: {active}") - - if args.cleanup_limit != 0: - cleanup_failures = cleanup_finished( - config.group_label, - dry_run=not args.execute, - limit=args.cleanup_limit, - ) - print(f"cleanup failures: {len(cleanup_failures)}") - if not cleanup_failures.empty: - print(cleanup_failures.to_string(index=False)) - - if active >= args.max_active: - print(f"skip submit: active workchains >= max_active ({args.max_active})") - return 0 - - local_inputs = next_unprocessed_inputs( - config=config, - cif_dir=args.cif_dir, - convergence_summary=args.convergence_summary, - full_scf_summary=args.full_scf_summary, - batch_size=args.batch_size, - ultra_extra_kindex=args.ultra_extra_kindex, - ) - print(f"next unprocessed source_db_ids: {len(local_inputs)}") - for item in local_inputs[:10]: - print(f" {item.source_db_id}: submit 0..{item.kindex_max}") - - if not args.execute: - print("dry run: pass --execute to submit jobs and delete cleanup files") - return 0 - - summary = submit_jobs( - (build_job_spec(item, args.cif_dir, args.source_db) for item in local_inputs), - config, - ) - print( - f"submitted={len(summary.submitted)} " - f"skipped_existing={len(summary.skipped_existing)} " - f"failed_source_db_ids={len(summary.failed_sources)}" - ) - if summary.failed_sources: - for failed in summary.failed_sources: - print(f" failed {failed.source_db_id} {failed.stage}: {failed.reason}") - return 0 - - -def parse_args() -> argparse.Namespace: - parser = argparse.ArgumentParser() - parser.add_argument("--cif-dir", type=Path, required=True) - parser.add_argument("--convergence-summary", type=Path, required=True) - parser.add_argument("--full-scf-summary", type=Path, required=True) - parser.add_argument("--source-db", default=DEFAULT_SOURCE_DB) - parser.add_argument("--execute", action="store_true") - parser.add_argument("--batch-size", type=int, default=100) - parser.add_argument("--max-active", type=int, default=1000) - parser.add_argument("--cleanup-limit", type=int, default=500) - parser.add_argument("--ultra-extra-kindex", type=int, default=2) - parser.add_argument("--code-label", default="qe-7.5-pw-admin@scarf") - parser.add_argument("--pseudo-family-label", default="PseudoDojo/0.4/PBEsol/SR/standard/upf") - parser.add_argument("--group-label", default="goldilocks/qe-scf/nospin/pseudodojo") - parser.add_argument("--degauss-ry", type=float, default=0.01) - parser.add_argument("--num-machines", type=int, default=1) - parser.add_argument("--num-mpiprocs-per-machine", type=int, default=32) - parser.add_argument("--max-wallclock-seconds", type=int, default=7200) - return parser.parse_args() - - -def next_unprocessed_inputs( - *, - config: AiidaScfConfig, - cif_dir: Path, - convergence_summary: Path, - full_scf_summary: Path, - batch_size: int, - ultra_extra_kindex: int, -) -> list[LocalJobInput]: - summary = load_summary(convergence_summary, SUMMARY_NOTE_COLUMN) - full_summary = load_summary(full_scf_summary, FULL_NOTE_COLUMN) - failed_ids = failed_source_ids(config) - ultra_rows = ultra_rows_by_source(summary, full_summary) - eligible_ids = eligible_source_ids(ultra_rows, cif_dir) - existing_kindices = existing_kindices_by_source(config.group_label) - print( - f"scan setup: {len(eligible_ids)} eligible source_db_ids, " - f"{len(existing_kindices)} source_db_ids with existing workchains, " - f"{len(failed_ids)} failed source_db_id records" - ) - - items: list[LocalJobInput] = [] - skipped_failed = 0 - skipped_complete = 0 - for scanned, source_db_id in enumerate(eligible_ids, start=1): - if source_db_id in failed_ids: - skipped_failed += 1 - if scanned % 100 == 0: - print_scan_progress(scanned, len(eligible_ids), len(items), skipped_failed, skipped_complete) - continue - row = ultra_rows[source_db_id] - ultra_mesh = (int(row["k1"]), int(row["k2"]), int(row["k3"])) - pmg_structure = load_pmg_structure(source_db_id, cif_dir) - new_ultra = new_kindex_for_mesh(pmg_structure, ultra_mesh) - kindex_max = new_ultra + int(ultra_extra_kindex) - existing = existing_kindices.get(source_db_id, set()) - if source_has_expected_kindices(existing, 0, kindex_max): - skipped_complete += 1 - if scanned % 100 == 0: - print_scan_progress(scanned, len(eligible_ids), len(items), skipped_failed, skipped_complete) - continue - items.append(LocalJobInput(source_db_id, kindex_max)) - if scanned % 100 == 0: - print_scan_progress(scanned, len(eligible_ids), len(items), skipped_failed, skipped_complete) - if len(items) >= batch_size: - break - print_scan_progress(scanned if eligible_ids else 0, len(eligible_ids), len(items), skipped_failed, skipped_complete) - return items - - -def source_has_expected_kindices(existing: set[int], kindex_min: int, kindex_max: int) -> bool: - return set(range(int(kindex_min), int(kindex_max) + 1)).issubset(existing) - - -def print_scan_progress( - scanned: int, - total: int, - selected: int, - skipped_failed: int, - skipped_complete: int, -) -> None: - print( - f"scan progress: {scanned}/{total} scanned, " - f"{selected} selected, " - f"{skipped_complete} complete, " - f"{skipped_failed} failed-record source_db_ids" - ) - - -def load_summary(path: Path, note_column: str) -> pd.DataFrame: - table = pd.read_csv(path) - table["source_db_id"] = table["source_db_id"].astype(str) - table[note_column] = table[note_column].astype(str) - return table - - -def eligible_source_ids(ultra_rows: dict[str, pd.Series], cif_dir: Path) -> list[str]: - cif_ids = {path.stem for path in cif_dir.glob("*.cif")} - return sorted(cif_ids & set(ultra_rows)) - - -def ultra_rows_by_source(summary: pd.DataFrame, full_summary: pd.DataFrame) -> dict[str, pd.Series]: - rows: dict[str, pd.Series] = {} - for table, note_column in [(summary, SUMMARY_NOTE_COLUMN), (full_summary, FULL_NOTE_COLUMN)]: - ultra_table = table.loc[table[note_column] == "ultra"].sort_values(["source_db_id", "k_index"]) - for source_db_id, group in ultra_table.groupby("source_db_id", sort=False): - row = group.iloc[-1] - previous = rows.get(source_db_id) - if previous is None or int(row["k_index"]) > int(previous["k_index"]): - rows[source_db_id] = row - return rows - - -def load_structure(source_db_id: str, cif_dir: Path, source_db: str) -> StructureData: - cif_path = cif_dir / f"{source_db_id}.cif" - atoms = read(cif_path) - structure = StructureData(ase=atoms) - structure.base.extras.set_many( - { - "source_db": source_db, - "source_db_id": source_db_id, - "source_file": str(cif_path), - } - ) - structure.store() - return structure - - -def load_pmg_structure(source_db_id: str, cif_dir: Path): - cif_path = cif_dir / f"{source_db_id}.cif" - atoms = read(cif_path) - return AseAtomsAdaptor.get_structure(atoms) - - -def new_kindex_for_mesh(pmg_structure, mesh: tuple[int, int, int]) -> int: - candidates = generate_candidate_k_distances(pmg_structure) - if not candidates: - raise ValueError(f"Mesh {mesh} not found in gamma-inclusive schedule") - - schedule = [k_distance_to_mesh(pmg_structure, candidates[0] + 1.0)] - schedule.extend( - k_distance_to_mesh(pmg_structure, 0.5 * (upper + lower)) - for upper, lower in zip(candidates[:-1], candidates[1:], strict=True) - ) - - seen: set[tuple[int, int, int]] = set() - for candidate_mesh in schedule: - if candidate_mesh in seen: - continue - if candidate_mesh == mesh: - return len(seen) - seen.add(candidate_mesh) - raise ValueError(f"Mesh {mesh} not found in gamma-inclusive schedule") - - -def build_job_spec(item: LocalJobInput, cif_dir: Path, source_db: str) -> AiidaJobSpec: - structure = load_structure(item.source_db_id, cif_dir, source_db) - pmg_structure = AseAtomsAdaptor.get_structure(structure.get_ase()) - return AiidaJobSpec( - source_db_id=item.source_db_id, - structure=structure, - code=DftCode.QE, - intent=CalculationIntent.SCF, - points=kindex_points(pmg_structure, 0, item.kindex_max), - ) - - -def active_workchain_count(group_label: str) -> int: - from aiida.orm import Group, QueryBuilder, WorkChainNode - - qb = QueryBuilder() - qb.append(Group, filters={"label": group_label}, tag="group") - qb.append( - WorkChainNode, - with_group="group", - filters={"attributes.process_state": {"in": ["created", "waiting", "running"]}}, - project="id", - ) - return qb.count() - - -if __name__ == "__main__": - raise SystemExit(main()) diff --git a/docs/published-records.md b/docs/published-records.md index 4e0bd34..adb1a7f 100644 --- a/docs/published-records.md +++ b/docs/published-records.md @@ -7,16 +7,18 @@ identifier and can be cited. These are snapshots. The AiiDA database remains the authoritative calculation record; a published dataset is a documented view of it at one point in time. -## Quantum ESPRESSO no-spin SCF calculations (SSSP, K-index) +## Quantum ESPRESSO no-spin SCF calculations (SSSP, k-index) -[`d5ds2-64f16`](https://data-collections.psdi.ac.uk/records/d5ds2-64f16) · v1 · +[`52713-55d86`](https://data-collections.psdi.ac.uk/records/52713-55d86) · v1.0 · CC BY 4.0 -The converged k-point mesh for 17,757 MC3D structures. No spin polarisation, -SSSP PBEsol pseudopotentials, every mesh unshifted and therefore -gamma-inclusive. +The current SSSP k-index dataset: the converged k-point mesh for 17,757 MC3D +structures, numbered on the **1-based** ladder (rung 1 the Γ-only `(1, 1, 1)` +mesh) and built with the resolution floor `min_k_distance = 0.03` Å⁻¹ rather +than a per-axis k-point cap. No spin polarisation, SSSP PBEsol pseudopotentials, +every mesh unshifted and therefore gamma-inclusive. -Convergence is the first of three consecutive k-distances whose total energies +Convergence is the first of three consecutive ladder meshes whose total energies agree within **1 meV per atom**. Energy only — no force criterion. | File | Contents | @@ -28,28 +30,9 @@ The archive carries more structures than the table has rows: 463 structures were calculated but never met the criterion within the range of meshes swept, so they have a structure file and no converged answer. -### The k_index in this record - -`k_index` in this record is **0-based, with rung 0 the gamma-only `(1, 1, 1)` -mesh**, and was computed with a per-axis enumeration bound of **50**, the ladder -truncated at the first rung where an axis count would rise by more than one. - -!!! warning "This record is 0-based; the convention since is 1-based" - - Everything produced after this record numbers the same ladder from 1, so - rung *n* here is rung *n + 1* under the current convention. The published - record is not rewritten: it keeps the convention it was published with, and - a consumer reads the base from the record rather than assuming it. - -That bound is part of the definition, not an implementation detail. The change -points of the ladder are `|b_i| / n`, the bound applies per axis, and axes with -different `|b_i|` exhaust their change points at different k-distances — so the -ladder is a complete set of meshes only down to `max(|b_i|) / 50`. Below that a -reachable mesh can be skipped. Any recomputation must state the bound it used or -its `k_index` values are not comparable with these. - -Raising the bound only appends rungs and never renumbers an existing one, so -these values stay valid under a larger enumeration. +`k_index` needs both its base and its floor to mean anything; the record's +`README.md` and `manifest.json` carry both, and a consumer reads them from there +rather than assuming. See [convergence criteria](reference/convergence.md) for how labels are assigned, and the record's own `README.md` for the full definition and reproduction code. diff --git a/docs/publishing.md b/docs/publishing.md index a3c6dde..c276869 100644 --- a/docs/publishing.md +++ b/docs/publishing.md @@ -49,7 +49,7 @@ convention a consumer must share are required fields: "rows": 17757, "columns": [{"name": "k_index", "dtype": "int", "description": "..."}], "conventions": { - "kmesh_ladder": {"base": 0, "rung_0": "gamma_only", "max_kpoints_per_axis": 50} + "kmesh_ladder": {"base": 1, "rung_1": "gamma_only", "min_k_distance": 0.03} }, "provenance": {"code": "quantum_espresso", "calculation": "scf", "spin": "none"} } diff --git a/docs/reference/kmesh.md b/docs/reference/kmesh.md index 126517b..7158c07 100644 --- a/docs/reference/kmesh.md +++ b/docs/reference/kmesh.md @@ -5,7 +5,7 @@ carrying six quantities. They describe the same mesh in different units, and two of them use different reciprocal-lattice conventions, so mixing them silently is easy. This page defines each one. -Built by `goldilocks_data.sweeps.kmesh.build_gamma_kmesh_entries`. +Built by `goldilocks_data.kmesh.build_gamma_kmesh_entries`. ## Two reciprocal conventions, both in use @@ -32,11 +32,9 @@ same scale. The rung's position, **1-based**, with rung 1 the Γ-only `(1, 1, 1)` mesh. Each step up is the next denser mesh the reciprocal lattice admits. -Record `d5ds2-64f16` predates this convention and is 0-based; see -[published records](../published-records.md). - -`kindex` only means something together with the enumeration bound that built the -ladder — see [the ladder](#the-ladder) below. +`kindex` only means something together with the resolution floor +`min_k_distance` that bounded the ladder — see [the ladder](#the-ladder) below. +Read the base and the floor from a record rather than assuming them. ## `mesh` @@ -55,13 +53,13 @@ n_i = max(1, ceil(|b_i| / k_distance)) The half-open range of k-distance that yields this mesh, as `(lower, upper)` in Å⁻¹ on the solid-state lengths. Any k-distance in it gives the same mesh. -A mesh corresponds to an interval, never a single value. Rung 0's upper bound is +A mesh corresponds to an interval, never a single value. Rung 1's upper bound is infinite: every k-distance above `max(|b_i|)` gives `(1, 1, 1)`. ```text -kindex 0 mesh (1, 1, 1) k_distance_interval (1.61061, inf) -kindex 1 mesh (2, 2, 1) k_distance_interval (0.87042, 1.61061) -kindex 20 mesh (14, 14, 8) k_distance_interval (0.11504, 0.12389) +kindex 1 mesh (1, 1, 1) k_distance_interval (1.61061, inf) +kindex 2 mesh (2, 2, 1) k_distance_interval (0.87042, 1.61061) +kindex 21 mesh (14, 14, 8) k_distance_interval (0.11504, 0.12389) ``` If you reduce the interval to one number for training, record which end you @@ -70,17 +68,13 @@ took. The two ends are different numbers for the same mesh. !!! warning "`entry_payload` names the ends the other way round" `entry_payload` emits `k_dist_left` for the interval's **lower** bound and - `k_dist_right` for its **upper** bound, so for `kindex 20` above it gives + `k_dist_right` for its **upper** bound, so for `kindex 21` above it gives `k_dist_left = 0.11504` and `k_dist_right = 0.12389`. - The published record - [`d5ds2-64f16`](https://data-collections.psdi.ac.uk/records/d5ds2-64f16) - uses the opposite orientation — its `k_dist_interval` is written - `[0.123, 0.115)`, larger value first, because a larger k-distance means a - coarser mesh. - - Joining notebook output with that record without checking will swap the - bounds. Compare magnitudes, not column names. Tracked as + Published records write `k_dist_interval` the other way round — larger value + first, `[0.12389, 0.11504)`, because a larger k-distance means a coarser + mesh. Joining notebook output with a record without checking will swap the + bounds: compare magnitudes, not column names. Tracked as [#30](https://github.com/stfc/goldilocks-data/issues/30). ## `k_line_density_interval` @@ -105,7 +99,7 @@ K-points per reciprocal atom: the full mesh size times the number of atoms. k_pra = n_atoms * n1 * n2 * n3 ``` -For `100115` at `kindex 20`: `4 * 14 * 14 * 8 = 6272`. It is a cost-like measure +For `100115` at `kindex 21`: `4 * 14 * 14 * 8 = 6272`. It is a cost-like measure that lets meshes be compared across cells of different sizes, and it uses the **full** mesh, not the symmetry-reduced count. @@ -113,7 +107,7 @@ that lets meshes be compared across cells of different sizes, and it uses the How many k-points survive symmetry reduction of the unshifted mesh, via pymatgen's `SpacegroupAnalyzer.get_ir_reciprocal_mesh`. This is what the -calculation actually costs. For `100115` at `kindex 20`: **120**, against a full +calculation actually costs. For `100115` at `kindex 21`: **120**, against a full mesh of 1568. !!! note "It falls back to the full mesh size" @@ -127,29 +121,30 @@ mesh of 1568. Change points are `|b_i| / n`, because `ceil(|b_i| / k_distance)` steps from `n` to `n + 1` exactly there. Sorted descending, each interval between neighbours is -one rung, and probing its midpoint gives the mesh. - -`max_kpoints_per_axis` (default **50**) bounds `n` — k-points per axis, **not** -the number of rungs, which is roughly the number of distinct axis lengths times -the bound. - -The bound applies per axis, and axes with different `|b_i|` exhaust their change -points at different k-distances, so the change-point list is complete only down -to `max(|b_i|) / max_kpoints_per_axis`. The ladder is therefore built with two -rules: - -- **truncate at the first gap** — stop at the first rung where an axis count - would rise by more than one, which is exactly the signature of a change point - that was never enumerated; -- **skip a repeat** — axes with equal `|b_i|` share change points, so two - consecutive intervals can yield the same mesh; keeping both would give one mesh - two `kindex` values. - -Raising the bound only appends rungs and never renumbers an existing one, so a -recorded `kindex` stays valid under a larger enumeration. A `kindex` computed -under a *different* bound is not comparable unless it lies in the region both -bounds cover — so any published `kindex` column must state its bound. - -This convention is shared with -[goldilocks-core](https://github.com/stfc/goldilocks-core). Changing it -invalidates every `kindex` already recorded or trained on. +one rung, and probing its midpoint gives the mesh. The first probe sits above +`max(|b_i|)` and yields the Γ-only mesh. + +`min_k_distance` (default **0.03 Å⁻¹**, on the solid-state 2π lengths — the +AiiDA-QuantumESPRESSO convention) is the resolution floor: change points denser +than it are not enumerated, so the ladder ends at a mesh of roughly +`ceil(|b_i| / min_k_distance)` per axis. + +The floor is identical for every axis, so the axes run out of change points +together. Every change point on `[min_k_distance, ∞)` is therefore present and +consecutive rungs differ by at most one k-point on each axis — the ladder has no +region where a reachable mesh is silently skipped. + +One rule still shapes it: **skip a repeat**. Axes with equal `|b_i|` share their +change points, and the two rounding precisions (`round(·, 8)` on the candidate +k-distances, `round(·, 5)` inside `k_distance_to_mesh`) can land two adjacent +intervals on one mesh; without the skip, that mesh would take two `kindex` +values. + +Lowering `min_k_distance` only appends rungs and never renumbers an existing +one, so a recorded `kindex` stays valid under a smaller floor. A `kindex` +computed under a *different* floor is comparable only where the two ranges +overlap — so any published `kindex` column must state its floor. + +This convention must stay in step with +[goldilocks-core](https://github.com/stfc/goldilocks-core) and goldilocks-ml. +Changing it invalidates every `kindex` already recorded or trained on. diff --git a/src/goldilocks_data/sweeps/kmesh.py b/src/goldilocks_data/kmesh.py similarity index 58% rename from src/goldilocks_data/sweeps/kmesh.py rename to src/goldilocks_data/kmesh.py index 9f41d4c..7d73d31 100644 --- a/src/goldilocks_data/sweeps/kmesh.py +++ b/src/goldilocks_data/kmesh.py @@ -4,6 +4,14 @@ from dataclasses import dataclass from typing import Any +from goldilocks_data.sweeps.models import SweepAxis, SweepPoint + +# Resolution floor for the k-mesh ladder, in Angstrom^-1 on the solid-state +# (2*pi) reciprocal lattice -- the same convention as AiiDA-QuantumESPRESSO +# k-distance. Change points denser than this are not enumerated, and this value +# is part of every recorded ``kindex``: a rung means nothing without it. +MIN_K_DISTANCE = 0.03 + @dataclass(frozen=True, slots=True) class KMeshEntry: @@ -30,25 +38,28 @@ def k_distance_to_mesh(structure: Any, k_distance: float) -> tuple[int, int, int return tuple(max(1, math.ceil(round(length / k_distance, 5))) for length in lengths) -def generate_candidate_k_distances(structure: Any, max_kpoints_per_axis: int = 50) -> list[float]: - """Return the k-distances at which any axis changes its k-point count. +def generate_candidate_k_distances(structure: Any, min_k_distance: float = MIN_K_DISTANCE) -> list[float]: + """Return every k-distance at which some axis changes its k-point count. - ``mesh_i = ceil(|b_i| / k_distance)`` steps from ``n`` to ``n + 1`` exactly at - ``k_distance = |b_i| / n``, so those quotients are the only distances where - the mesh can change. + ``mesh_i = ceil(|b_i| / k_distance)`` steps from ``n`` to ``n + 1`` exactly + at ``k_distance = |b_i| / n``, so those quotients are the only distances + where the mesh can change. - ``max_kpoints_per_axis`` bounds ``n`` — how many k-points per axis are - enumerated. It is *not* a bound on the number of rungs, which is roughly the - number of distinct axis lengths times the bound. Because the bound applies - per axis and the axes have different ``|b_i|``, they exhaust their quotients - at different distances: the returned list is a complete set of change points - only down to ``max(|b_i|) / max_kpoints_per_axis``. See - ``build_gamma_kmesh_entries`` for what that implies. + ``min_k_distance`` is the resolution floor in Angstrom^-1 on the solid-state + (2*pi) reciprocal lattice; quotients below it are not enumerated. The floor + is identical for every axis, so all axes run out of change points together + and the returned list is a complete set of change points over the whole + ``[min_k_distance, inf)`` range -- there is no region where a reachable mesh + is silently skipped. """ lengths = _reciprocal_lengths(structure) return sorted( - {round(length / index, 8) for length in lengths for index in range(1, max_kpoints_per_axis + 1)}, + { + round(length / index, 8) + for length in lengths + for index in range(1, max(1, math.floor(length / min_k_distance)) + 1) + }, reverse=True, ) @@ -76,7 +87,7 @@ def _n_reduced_kpoints(structure: Any, mesh: tuple[int, int, int]) -> int: return full_mesh_size -def build_gamma_kmesh_entries(structure: Any, max_kpoints_per_axis: int = 50) -> list[KMeshEntry]: +def build_gamma_kmesh_entries(structure: Any, min_k_distance: float = MIN_K_DISTANCE) -> list[KMeshEntry]: """Build the unshifted, Gamma-inclusive k-mesh ladder for a structure. ``kindex`` is 1-based and rung 1 is the Gamma-only ``(1, 1, 1)`` mesh, which @@ -84,23 +95,17 @@ def build_gamma_kmesh_entries(structure: Any, max_kpoints_per_axis: int = 50) -> The rung is an ordinal position on this structure's ladder and counts nothing. It equals the densest axis count only where the axes step together, - as in a cubic cell: for ``|b| = [1.6106, 1.6106, 0.8704]`` rung 3 is - ``(2, 2, 2)``, whose densest axis carries 2. The same rung is a different - mesh for a different cell. - - The ladder is complete and non-repeating: - - * It stops at the first rung where an axis count would rise by more than one. - Once the longest axis has used up its enumerated quotients its count keeps - rising with no candidate marking the change, so adjacent candidates span - several meshes and probing the midpoint keeps only one of them. A jump - greater than one is exactly that condition. - * It skips a mesh already on the ladder. Axes with equal ``|b_i|`` share - their change points, so two consecutive intervals can yield the same mesh; - without this, two ``kindex`` values would name one mesh. + as in a cubic cell. + + The ladder is complete and non-repeating down to ``min_k_distance``: every + change point above the floor is enumerated, so consecutive rungs differ by + at most one k-point on each axis and no reachable mesh is skipped. A mesh + already on the ladder is dropped -- axes with equal ``|b_i|`` share their + change points, so two consecutive intervals can yield the same mesh, and + without the skip two ``kindex`` values would name one mesh. """ - candidates = generate_candidate_k_distances(structure, max_kpoints_per_axis) + candidates = generate_candidate_k_distances(structure, min_k_distance) if not candidates: return [] @@ -110,11 +115,7 @@ def build_gamma_kmesh_entries(structure: Any, max_kpoints_per_axis: int = 50) -> entries: list[KMeshEntry] = [] seen: set[tuple[int, int, int]] = set() - previous: tuple[int, int, int] | None = None for mesh, interval in intervals: - if previous is not None and any(now - before > 1 for before, now in zip(previous, mesh, strict=True)): - break - previous = mesh if mesh in seen: continue seen.add(mesh) @@ -147,3 +148,29 @@ def entry_payload(entry: KMeshEntry) -> dict[str, Any]: "k_dist_left": float(left), "k_dist_right": None if math.isinf(right) else float(right), } + + +def kindex_points(structure: object, kindex_min: int, kindex_max: int) -> tuple[SweepPoint, ...]: + """Build explicit sweep points for a gamma-inclusive kindex range. + + ``kindex_min`` and ``kindex_max`` are rungs, not list positions: rung 1 is + the Gamma-only mesh and lives at index 0. + """ + + entries = build_gamma_kmesh_entries(structure) + selected = [entry for entry in entries if int(kindex_min) <= entry.kindex <= int(kindex_max)] + points: list[SweepPoint] = [] + for entry in selected: + points.append( + SweepPoint( + axis_values={SweepAxis.KINDEX.value: int(entry.kindex)}, + k_mesh=entry.mesh, + extras={ + "kindex": int(entry.kindex), + "k_mesh": list(entry.mesh), + "k_pra": entry.k_pra, + "n_reduced_kpoints": entry.n_reduced_kpoints, + }, + ) + ) + return tuple(points) diff --git a/src/goldilocks_data/sweeps/__init__.py b/src/goldilocks_data/sweeps/__init__.py index c923d74..d1ab497 100644 --- a/src/goldilocks_data/sweeps/__init__.py +++ b/src/goldilocks_data/sweeps/__init__.py @@ -1,20 +1,14 @@ """Sweep definitions independent of data source and execution backend.""" from goldilocks_data.sweeps.extension import ExtensionPlan, KindexExtension, plan_well_not_ultra_extensions -from goldilocks_data.sweeps.kindex import kindex_points -from goldilocks_data.sweeps.kmesh import KMeshEntry, build_gamma_kmesh_entries, entry_payload from goldilocks_data.sweeps.models import AiidaJobSpec, ScfSweepSpec, SweepAxis, SweepPoint __all__ = [ "AiidaJobSpec", "ExtensionPlan", - "KMeshEntry", "KindexExtension", "ScfSweepSpec", "SweepAxis", "SweepPoint", - "build_gamma_kmesh_entries", - "entry_payload", - "kindex_points", "plan_well_not_ultra_extensions", ] diff --git a/src/goldilocks_data/sweeps/kindex.py b/src/goldilocks_data/sweeps/kindex.py deleted file mode 100644 index 89d579c..0000000 --- a/src/goldilocks_data/sweeps/kindex.py +++ /dev/null @@ -1,30 +0,0 @@ -from __future__ import annotations - -from goldilocks_data.sweeps.kmesh import build_gamma_kmesh_entries -from goldilocks_data.sweeps.models import SweepAxis, SweepPoint - - -def kindex_points(structure: object, kindex_min: int, kindex_max: int) -> tuple[SweepPoint, ...]: - """Build explicit sweep points for a gamma-inclusive kindex range. - - ``kindex_min`` and ``kindex_max`` are rungs, not list positions: rung 1 is - the Gamma-only mesh and lives at index 0. - """ - - entries = build_gamma_kmesh_entries(structure) - selected = [entry for entry in entries if int(kindex_min) <= entry.kindex <= int(kindex_max)] - points: list[SweepPoint] = [] - for entry in selected: - points.append( - SweepPoint( - axis_values={SweepAxis.KINDEX.value: int(entry.kindex)}, - k_mesh=entry.mesh, - extras={ - "kindex": int(entry.kindex), - "k_mesh": list(entry.mesh), - "k_pra": entry.k_pra, - "n_reduced_kpoints": entry.n_reduced_kpoints, - }, - ) - ) - return tuple(points) diff --git a/tests/conftest.py b/tests/conftest.py index 5bee9c9..0e209c0 100644 --- a/tests/conftest.py +++ b/tests/conftest.py @@ -18,7 +18,7 @@ {"name": "source_db_id", "dtype": "str"}, {"name": "k_index", "dtype": "int", "definition": "kmesh_ladder_rung"}, ], - "conventions": {"kmesh_ladder": {"base": 0, "max_kpoints_per_axis": 50}}, + "conventions": {"kmesh_ladder": {"base": 1, "min_k_distance": 0.03}}, "provenance": {"code": "quantum_espresso", "calculation": "scf", "spin": "none"}, } diff --git a/tests/test_kmesh.py b/tests/test_kmesh.py index 4e5d740..48be2aa 100644 --- a/tests/test_kmesh.py +++ b/tests/test_kmesh.py @@ -2,7 +2,7 @@ from dataclasses import dataclass -from goldilocks_data.sweeps.kmesh import ( +from goldilocks_data.kmesh import ( build_gamma_kmesh_entries, entry_payload, k_distance_to_mesh, @@ -40,7 +40,7 @@ def test_gamma_kmesh_entries_start_with_gamma_mesh() -> None: natoms=4, ) - entries = build_gamma_kmesh_entries(structure, max_kpoints_per_axis=4) + entries = build_gamma_kmesh_entries(structure, min_k_distance=0.25) assert entries[0].kindex == 1 assert entries[0].mesh == (1, 1, 1) @@ -69,7 +69,7 @@ def test_entry_payload_serializes_infinite_right_bound_as_none() -> None: natoms=2, ) - payload = entry_payload(build_gamma_kmesh_entries(structure, max_kpoints_per_axis=2)[0]) + payload = entry_payload(build_gamma_kmesh_entries(structure, min_k_distance=0.5)[0]) assert payload["kindex"] == 1 assert payload["k_mesh"] == (1, 1, 1) @@ -86,11 +86,12 @@ def _structure(a: float, b: float, c: float, natoms: int = 1) -> Structure: ) -def test_ladder_never_skips_a_reachable_mesh() -> None: - # An anisotropic cell: the long axes exhaust their change points while the - # short one is still stepping, which is where an unbounded ladder starts - # jumping several meshes at once. - entries = build_gamma_kmesh_entries(_structure(2.5547, 2.5547, 0.6485), max_kpoints_per_axis=30) +def test_ladder_is_gap_free_down_to_the_floor() -> None: + # An anisotropic cell: a per-axis k-point cap would let the long axes run + # out of change points while the short one is still stepping. A k-distance + # floor stops every axis at the same place, so the ladder is gap-free by + # construction. + entries = build_gamma_kmesh_entries(_structure(2.5547, 2.5547, 0.6485), min_k_distance=0.09) meshes = [entry.mesh for entry in entries] assert len(meshes) > 1 @@ -110,15 +111,23 @@ def test_ladder_never_repeats_a_mesh_for_degenerate_axes() -> None: assert meshes[:4] == [(1, 1, 1), (1, 2, 1), (2, 2, 2), (2, 3, 2)] -def test_raising_the_axis_bound_only_extends_the_ladder() -> None: - # kindex is recorded in campaign snapshots and published records, so a - # larger enumeration must never renumber a rung that already existed. +def test_lowering_the_floor_only_extends_the_ladder() -> None: + # kindex is recorded in campaign snapshots and published records, so a lower + # floor must never renumber a rung that already existed. structure = _structure(2.5547, 2.5547, 0.6485) - short = [entry.mesh for entry in build_gamma_kmesh_entries(structure, max_kpoints_per_axis=20)] - long = [entry.mesh for entry in build_gamma_kmesh_entries(structure, max_kpoints_per_axis=60)] + coarse = [entry.mesh for entry in build_gamma_kmesh_entries(structure, min_k_distance=0.1)] + fine = [entry.mesh for entry in build_gamma_kmesh_entries(structure, min_k_distance=0.02)] - assert len(long) > len(short) - assert long[: len(short)] == short + assert len(fine) > len(coarse) + assert fine[: len(coarse)] == coarse + + +def test_ladder_stops_at_the_resolution_floor() -> None: + # The floor caps the densest mesh at ceil(|b_i| / min_k_distance) per axis, + # independent of any k-point count. 0.125 = 1/8 exactly, so no float noise. + entries = build_gamma_kmesh_entries(_structure(1.0, 1.0, 1.0), min_k_distance=0.125) + + assert entries[-1].mesh == (8, 8, 8) def test_kindex_is_contiguous_and_one_based() -> None: diff --git a/tests/test_publish_deposit.py b/tests/test_publish_deposit.py index 29d671b..87e89a3 100644 --- a/tests/test_publish_deposit.py +++ b/tests/test_publish_deposit.py @@ -126,11 +126,11 @@ def test_validate_dataset_record_rejects_a_column_without_a_dtype(dataset_record def test_validate_dataset_record_keeps_the_kmesh_convention(dataset_record: dict) -> None: # The convention block is why this file exists: a k_index with no stated - # base and no stated enumeration bound is not reproducible. + # base and no stated resolution floor is not reproducible. ladder = validate_dataset_record(dataset_record)["conventions"]["kmesh_ladder"] - assert ladder["base"] == 0 - assert ladder["max_kpoints_per_axis"] == 50 + assert ladder["base"] == 1 + assert ladder["min_k_distance"] == 0.03 def test_parse_sha256sums_accepts_binary_mode_lines() -> None: diff --git a/tests/test_sweep_models.py b/tests/test_sweep_models.py index ec08d3e..49b9e3f 100644 --- a/tests/test_sweep_models.py +++ b/tests/test_sweep_models.py @@ -4,8 +4,8 @@ from goldilocks_data.codes import DftCode from goldilocks_data.intents import CalculationIntent +from goldilocks_data.kmesh import kindex_points from goldilocks_data.sweeps import AiidaJobSpec, SweepAxis -from goldilocks_data.sweeps.kindex import kindex_points @dataclass(frozen=True, slots=True)