diff --git a/.gitignore b/.gitignore index cc3120282..0e5f65fed 100644 --- a/.gitignore +++ b/.gitignore @@ -46,6 +46,7 @@ __pycache__/ *.py[cod] *$py.class .vscode +.codegraph/ # C extensions *.so dptb/tests/data/esk_orb/sktb.json diff --git a/docs/maintenance/scc_numerics_test_plan.md b/docs/maintenance/scc_numerics_test_plan.md new file mode 100644 index 000000000..f92fef4a1 --- /dev/null +++ b/docs/maintenance/scc_numerics_test_plan.md @@ -0,0 +1,210 @@ +# Carbon BCC Self-Consistent Total-Energy K-Point Convergence Test Plan + +## Status + +In progress. The current scope is the convergence of the self-consistent total +energy of one standalone Carbon BCC calculation with respect to k-point +density. Inner SCC charge convergence is a required diagnostic but is not the +target convergence criterion. This is a validation plan, not an implementation +proposal. + +## Motivation + +DeePTB's SCC-DFTB functionality should give a stable, reproducible total energy +for a fixed Carbon BCC input as the k-point density is increased. An individual +mesh reaching its SCC charge tolerance only proves that mesh was solved +self-consistently; it does not prove Brillouin-zone integration or total-energy +convergence. This work is not a direct test against the existing multi- +structure DFTB+ pytest benchmark. + +The downstream UniSK workflow exposed a practical concern. Its legacy Carbon +Bayesian-optimization script uses per-structure k-point spacings such as +`[0.05, 0.05, 0.05]` for BCC and FCC. For the Carbon BCC cell used by the +DeePTB DFTB+ benchmark, that spacing generates a `75 x 75 x 75` mesh +(`421875` raw points), whereas the existing DeePTB benchmark uses +`20 x 20 x 20` (`8000` raw points). The former was copied from an old +production-oriented workflow with the comment that it was chosen to "ensure +convergence"; no k-point convergence study was preserved with that script. + +This difference must not be interpreted as proof that DeePTB SCC is slow or +incorrect. It is first a numerical-configuration question. A mesh that is +overly fine can make a correct SCC calculation impractical inside a Bayesian +optimization loop. + +## Existing Assets + +The repository contains a broader correctness regression: + +- `dptb/tests/test_dftb_scc.py::test_dftbscc_matches_dftbp_benchmarks` + compares DeePTB SCC electronic energies against DFTB+ reference data for + Carbon BCC, graphene, and dimer EOS scans. +- The scan contains seven lattice/bond-length scales for each structure. +- DFTB+ reference tables and the corresponding structures/SK files live in + `dptb/tests/data/dftb/structs_eos` and `dptb/tests/data/dftb`. +- The current acceptance threshold is a maximum electronic-energy difference + below `1e-2 eV` at every scale. +- Its direct command is: + + ```bash + python -m pytest -q dptb/tests/test_dftb_scc.py -k matches_dftbp_benchmarks + ``` + +That regression is useful background and supplies the canonical Carbon +structure and SK data, but it is not the execution target for the present +study. The current standalone entry point is: + +- `examples/dftb_scc/carbon_bcc_scc_convergence/run.py`; +- `examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR`, containing only the + scale-1.000 primitive Carbon BCC cell; and +- the canonical `C-C.skf` under `dptb/tests/data/dftb`, shared without invoking + any pytest helper or multi-structure EOS scan. + +The standalone example records self-consistent total, electronic, and repulsive +energies; energy changes between successive meshes and against all denser +meshes in the scan; k-point counts; SCC iteration count and charge residual +history; electron count; Mulliken charge; Fermi level; and wall time. Existing +unit tests for SCC state, mixers, API equivalence, and k-point symmetry +reduction remain useful but do not define this study's convergence criterion. + +## Questions To Answer + +The following effects are coupled in an SCC calculation but need separate +tests and diagnostics: + +1. **K-point integration:** What mesh keeps the self-consistent total energy per + atom within the agreed tolerance of every denser mesh tested for the fixed + Carbon BCC cell? +2. **Occupations and smearing:** Are Fermi-Dirac and Gaussian paths correctly + normalized, electron-number conserving, and stable in the relevant + zero-/finite-temperature limits? +3. **Fermi-level solve:** Does the root solve use the same k-point weights and + occupation convention as the band-energy and charge calculations? +4. **Symmetry reduction:** Do full, time-reversal-reduced, and + rotational-symmetry-reduced meshes give equivalent observables with correct + accumulated weights? +5. **Inner SCC loop:** Does every k mesh reach the charge-residual tolerance, + and can SCC numerical noise contaminate the outer total-energy trend? +6. **Performance:** Which component dominates runtime after the numerical + configuration is fixed: Hamiltonian construction, eigensolve, occupation, + Mulliken population, Gamma/SCC shift, or mixing? + +## Important Configuration Distinction + +The current UniSK Carbon workflow and the DeePTB DFTB+ benchmark are not a +direct apples-to-apples timing or numerical comparison. Among other settings, +the legacy UniSK path uses Gaussian smearing (`"G"`) and a `0.05` spacing, +while the DFTB+ benchmark uses Fermi-Dirac smearing (`"FD"`) and fixed meshes +such as Carbon BCC `20 x 20 x 20`. Future comparisons must explicitly record +all physical and numerical settings, including SK files, repulsive parameters, +cutoffs, temperature, mesh, symmetry flags, mixer, tolerance, and maximum +iteration count. + +## Phased Work + +### Phase 0: Lock the standalone Carbon BCC baseline + +- Use only the scale-1.000 Carbon BCC structure in the standalone example. +- Do not call `test_dftbscc_matches_dftbp_benchmarks` or scan graphene, dimer, + or the seven-point Carbon EOS in this phase. +- Record wall time, SCC iteration count, full residual history, final residual, + electron count, Fermi level, total/electronic/repulsive energies, SCC shift + energy, and Mulliken charge for every numerical case. +- Keep the command runnable independently of notebooks and without generating + or duplicating SK tables. + +### Phase 1: K-point convergence matrix + +For the fixed Carbon BCC cell, scan deliberately chosen fixed meshes and, where +useful, equivalent spacing-derived meshes. At minimum record: + +- total-energy change per atom from the preceding mesh; +- maximum absolute total-energy change per atom relative to every denser mesh + in the same run; +- an explicit k-point convergence flag evaluated against a documented energy + threshold; the last and densest point cannot validate itself; +- charge and Fermi-level change relative to the finest accepted mesh; +- SCC iteration count and convergence failures; +- raw and symmetry-reduced k-point counts; +- wall time split by SCC iteration and total calculation. + +Use the results to define separate **smoke**, **development**, and +**production** mesh presets. Do not retain `0.05` as a default solely because +it appeared in a legacy script. + +The initial `20^3, 24^3, 28^3, 32^3, 36^3, 40^3` even-grid scan shows that +`20^3` differs from a denser result by up to `2.252 meV/atom`. It therefore does +not meet a `1 meV/atom` criterion. `24^3` and denser tested even grids remain +within `1 meV/atom`, but neighbouring odd grids and spacing-derived meshes must +still be checked before fixing a production preset. + +### Phase 2: Occupation and Fermi-level regression tests + +- Add small, deterministic tests for occupation normalization and exact + electron-number conservation with non-uniform k-point weights. +- Test Fermi-Dirac and Gaussian implementations independently. +- Verify expected low-temperature and high-temperature limiting behaviour. +- First record Fermi level, band energy, SCC charge, and electron-count error + for the standalone Carbon BCC input. Broader material-class coverage is a + separate follow-up. + +### Phase 3: Symmetry-reduction equivalence + +- For the Carbon BCC cell, compare full meshes with time-reversal and + rotational reductions. +- Assert equivalence of electron count, Fermi level, charge, total/electronic + energy, SCC shift energy, and converged SCC shift within explicit tolerances. +- Assert that the sum of reduced k-point weights equals the full mesh weight. +- Report the reduction factor as a performance metric. + +Symmetry reduction is a performance optimization, not a change to the target +physical result. + +### Phase 4: SCC robustness and performance work + +Only after Phases 0-3 identify the limiting configuration should solver work +begin. Candidate changes include static-data caching across SCC iterations, +batched k-point eigensolves, warm starts, and mixer/preconditioner tuning. +Every change must run the numerical regression matrix and a separate timing +benchmark. + +The primitive Carbon BCC cell has one atom, so charge conservation and symmetry +largely pin its Mulliken charge. It can reveal occupation-noise effects in a +very tight SCC tolerance, but it is not a difficult charge-mixing problem. A +true mixer stress test will require a separate cell with non-equivalent charge +degrees of freedom or an explicitly perturbed initial charge; that must not be +silently substituted into this Carbon BCC baseline. + +## Tetrahedron Integration + +Tetrahedron integration may be useful for zero-temperature Brillouin-zone +integration, DOS, and difficult metallic cases, but it is not a drop-in +replacement for smearing. It requires a separate design covering band +interpolation, degeneracies, weights, self-consistency, and eventually force +consistency. It should be evaluated only after the existing Fermi-Dirac and +Gaussian paths are fully benchmarked and tested. + +## Acceptance Criteria For Future Changes + +A proposed SCC numerical or performance change is acceptable only when: + +1. the standalone Carbon BCC total-energy matrix remains within its documented + meV/atom tolerance against all denser meshes tested, while inner SCC charge + residuals also meet their separate tolerance; +2. symmetry-reduced and full-mesh results agree within documented tolerances; +3. electron-number conservation is verified for every occupation path; +4. convergence failures and iteration counts do not regress without a recorded + numerical reason; and +5. any claimed speedup reports the mesh, symmetry flags, precision, hardware, + warm-up policy, and the same numerical acceptance criteria. + +## Non-Goals + +- Covering graphene, the Carbon dimer, or the seven-point Carbon EOS in the + current standalone SCC convergence study. +- Using `test_dftbscc_matches_dftbp_benchmarks` as the current execution entry + point. +- Changing default UniSK production parameters in this task. +- Declaring a legacy `kmeshspacing` value physically necessary without a + convergence study. +- Replacing the SCC algorithm with a tetrahedron method before the current + occupation and k-point paths have a complete regression suite. diff --git a/dptb/nn/dftb/dftb_scc.py b/dptb/nn/dftb/dftb_scc.py index 5b7c2ce7a..3480ded32 100644 --- a/dptb/nn/dftb/dftb_scc.py +++ b/dptb/nn/dftb/dftb_scc.py @@ -85,7 +85,7 @@ def reset(self) -> None: structures, as it avoids the overhead of reinitializing the expensive model, SK parameters, and transformers for each structure. - The following attributes are reset to None: + The following per-calculation state is cleared or reset: - atomic_numbers: Atomic numbers of the current structure - elec_totE, elec_H0_bandE, elec_bandE: Electronic energies - E_fermi: Fermi energy @@ -96,6 +96,9 @@ def reset(self) -> None: - inv_r: Inverse distance matrix - Gamma: Coulomb interaction matrix - total_energy, total_rep_energy: Total energies + - n_iterations: Number of Mulliken/SCC iterations performed + - charge_residual_history: Maximum charge residual from each SCC update + - final_charge_residual: Last maximum charge residual, if one was evaluated Also resets the Mulliken calculator state via mulliken.reset(). ''' @@ -115,6 +118,9 @@ def reset(self) -> None: self.total_energy = None self.total_rep_energy = None self.is_converged = False # Track SCC convergence status + self.n_iterations = 0 + self.charge_residual_history = [] + self.final_charge_residual = None # Also reset the Mulliken calculator state self.mulliken.reset() @@ -218,6 +224,12 @@ def run_iters(self, - self.elec_totE: Total electronic energy (eV) - self.E_fermi: Fermi energy (eV) - self.mulliken.mul_charge: Final Mulliken charges + - self.is_converged: Whether the charge residual reached ``tol`` + - self.n_iterations: Number of Mulliken/SCC iterations performed + - self.charge_residual_history: Maximum charge residual for each + update after the initial charge evaluation + - self.final_charge_residual: Last value in the residual history, + or ``None`` when no residual was evaluated ''' @@ -328,10 +340,13 @@ def run_iters(self, # Convergence criterion: ||q_diff|| = ||q_out - q_inp|| (DFTBplus convention) diff = np.abs(q_diff) - log.debug(f' ITERATION: {int(iteration)} Max charge difference: {diff.max().item():.15f} e') + max_charge_residual = float(diff.max()) + self.charge_residual_history.append(max_charge_residual) + self.final_charge_residual = max_charge_residual + log.debug(f' ITERATION: {int(iteration)} Max charge difference: {max_charge_residual:.15f} e') iteration += 1 - if diff.max() < tol: # Converged and not exiting due to max_iter + if max_charge_residual < tol: # Converged and not exiting due to max_iter assert self.scc_shift_energy is not None, "scc_shift_energy should not be None when converged." assert self.elec_H0_bandE is not None, "elec_H0_bandE should not be None when converged." log.debug(f'Convergence reached after {int(iteration)} iters.') @@ -374,6 +389,7 @@ def run_iters(self, data[AtomicDataDict.EIGENVECTOR_KEY] = torch.cat(eigvecs, dim=0) self.data = data + self.n_iterations = iteration self.E_fermi = self.mulliken.estimated_E_fermi self.elec_bandE = self.mulliken.elec_bandE self.elec_totE = self.elec_H0_bandE + self.scc_shift_energy diff --git a/dptb/tests/test_dftb_scc.py b/dptb/tests/test_dftb_scc.py index 7203b5809..2875b8987 100644 --- a/dptb/tests/test_dftb_scc.py +++ b/dptb/tests/test_dftb_scc.py @@ -478,6 +478,9 @@ def test_dftbscc_reset_clears_state(rootdir=rootdir): assert dftbscc.expGamma_onsite is not None, "expGamma_onsite should be set after calculation" assert dftbscc.inv_r is not None, "inv_r should be set after calculation" assert dftbscc.Gamma is not None, "Gamma should be set after calculation" + assert dftbscc.n_iterations > 0, "n_iterations should be populated after calculation" + assert len(dftbscc.charge_residual_history) == dftbscc.n_iterations - 1 + assert dftbscc.final_charge_residual == dftbscc.charge_residual_history[-1] # Call reset explicitly dftbscc.reset() @@ -498,6 +501,9 @@ def test_dftbscc_reset_clears_state(rootdir=rootdir): assert dftbscc.Gamma is None, "Gamma should be None after reset" assert dftbscc.total_energy is None, "total_energy should be None after reset" assert dftbscc.total_rep_energy is None, "total_rep_energy should be None after reset" + assert dftbscc.n_iterations == 0, "n_iterations should be zero after reset" + assert dftbscc.charge_residual_history == [], "residual history should be empty after reset" + assert dftbscc.final_charge_residual is None, "final residual should be None after reset" # Verify persistent state (model, skp, etc.) is NOT cleared assert dftbscc.model is not None, "model should persist after reset" diff --git a/examples/dftb_scc/README.md b/examples/dftb_scc/README.md index 16192091a..f82866a69 100644 --- a/examples/dftb_scc/README.md +++ b/examples/dftb_scc/README.md @@ -6,6 +6,10 @@ Imports and obvious absolute/local SK paths were adjusted to use DeePTB modules Recommended notebooks: +- `carbon_bcc_scc_convergence/run.py`: standalone Carbon BCC self-consistent + total-energy convergence against k meshes/spacings. Inner SCC, mixer, and + smearing diagnostics are recorded separately. See the example README for + commands and interpretation. - `benchmark_eos_dftbp.ipynb`: SCC EOS comparison against DFTB+ reference data in `dptb/tests/data/dftb`. - `hBN_scc/test_scc_hBN.ipynb`: hBN SCC workflow and band calculation. - `hBN_scc/test_scc_hBN_nnsk.ipynb`: hBN nonorthogonal NNSK-SCC workflow. diff --git a/examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR b/examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR new file mode 100644 index 000000000..a4dfdcb24 --- /dev/null +++ b/examples/dftb_scc/carbon_bcc_scc_convergence/POSCAR @@ -0,0 +1,9 @@ +C BCC primitive cell, scale 1.000 +1.0000000000000000 + 2.0557030981993387 0.0000005029268204 0.0000008733164838 + -0.6852331938919997 1.9381358756482154 0.0000008681559103 + -0.6852331920572658 -0.9690657046357299 1.6784761976400846 +C +1 +Cartesian + -0.0000000006781663 0.0000000005972072 -0.0000000000918888 diff --git a/examples/dftb_scc/carbon_bcc_scc_convergence/README.md b/examples/dftb_scc/carbon_bcc_scc_convergence/README.md new file mode 100644 index 000000000..781c50fca --- /dev/null +++ b/examples/dftb_scc/carbon_bcc_scc_convergence/README.md @@ -0,0 +1,121 @@ +# Carbon BCC self-consistent total-energy k-point convergence + +This is a standalone numerical diagnostic for one fixed Carbon BCC structure. +It does **not** invoke +`dptb/tests/test_dftb_scc.py::test_dftbscc_matches_dftbp_benchmarks`, scan an +EOS, or include graphene and the Carbon dimer. + +There are two distinct convergence loops: + +1. **Inner SCC convergence**: at one fixed k mesh, the charge residual reaches + `--tol`. This only says that particular discrete Brillouin-zone problem was + solved self-consistently. +2. **Outer k-point convergence**: the converged total energy per atom stops + changing as the k mesh is made denser. This is the primary criterion in this + example. + +An SCC-converged result at `20 x 20 x 20` is not automatically a k-point- +converged result. + +The script reports, for every case: + +- requested and resolved k-point mesh, raw and symmetry-reduced k-point counts, + and the reduced-weight sum; +- self-consistent total energy, electronic energy, and repulsive energy, using + `sigma_rep = {"C": 0.561}`; +- total-energy changes from the preceding mesh and maximum absolute change + against every denser mesh in the same scan; +- `kpoint_energy_converged_within_scan`, evaluated against the configurable + `--energy-tol-mev-per-atom` threshold (default `1 meV/atom`); +- inner SCC status, iteration count, final charge residual, and the full + residual history in JSON output; +- electron count/error, Mulliken and delta charges, Fermi level, electronic + energy, H0 band energy, and SCC shift energy; +- wall time for the SCC calculation. + +The densest point in a scan has no denser validation point and is therefore +reported with `kpoint_energy_converged_within_scan = null`, rather than being +declared converged automatically. For this fixed structure, the repulsive +energy is independent of k sampling, so total- and electronic-energy changes +are identical; the script still calculates and records the actual total energy. + +The committed BCC POSCAR is the scale-1.000 primitive cell used by the existing +Carbon data. The script shares the repository's canonical `C-C.skf` file under +`dptb/tests/data/dftb` by default, avoiding a second copy of the same 104 KiB +table. Override it with `--sk-path` if another Carbon parameter set is being +tested. + +## Important limitation + +This primitive BCC cell contains one Carbon atom. Charge conservation and +symmetry therefore pin its Mulliken population close to 4 electrons, so SCC +charge mixing normally converges in the first residual check. The example is a +valid total-energy k-point convergence test even though its inner SCC loop is +trivial. It is not, by itself, a difficult mixer-robustness benchmark. A later +mixer stress test needs a structure with non-equivalent charge degrees of +freedom or an explicitly perturbed initial charge. + +## Run + +From the repository root: + +```bash +uv run python examples/dftb_scc/carbon_bcc_scc_convergence/run.py \ + --scan kmesh \ + --meshes 20 24 28 32 36 40 \ + --energy-tol-mev-per-atom 1.0 \ + --output /tmp/carbon_bcc_kmesh.json +``` + +If `uv` is unavailable but the repository virtual environment already exists: + +```bash +.venv/bin/python examples/dftb_scc/carbon_bcc_scc_convergence/run.py \ + --scan kmesh --meshes 12 16 20 24 +``` + +Preview spacing-derived meshes without allocating or diagonalizing them: + +```bash +uv run python examples/dftb_scc/carbon_bcc_scc_convergence/run.py \ + --scan spacing --spacings 0.20 0.10 0.05 --dry-run +``` + +The `0.05` case should be previewed first because it resolves to a very large +mesh for this cell. It is intentionally not included in the default spacing +list and must be requested explicitly. + +Compare mixers at one fixed mesh: + +```bash +uv run python examples/dftb_scc/carbon_bcc_scc_convergence/run.py \ + --scan mixer --fixed-mesh 8 \ + --mixers simple pulay broyden diis \ + --mix-rates 0.1 0.3 \ + --output /tmp/carbon_bcc_mixers.json +``` + +Compare Fermi-Dirac and Gaussian occupations: + +```bash +uv run python examples/dftb_scc/carbon_bcc_scc_convergence/run.py \ + --scan smearing --fixed-mesh 8 \ + --smearing-methods FD G --temperatures 0.1 300 +``` + +Use `--help` for convergence tolerance, maximum iterations, symmetry controls, +structure override, and other numerical settings. The default charge-residual +tolerance is `1e-10 e`; tighter tolerances such as `1e-12` or the historical +`1e-14` should be requested explicitly because occupation/charge numerical +noise can otherwise dominate the inner SCC iteration count. This charge +tolerance is independent of `--energy-tol-mev-per-atom`, which controls the +outer k-point convergence assessment. + +## Current observation + +For the tested even-grid sequence `20, 24, 28, 32, 36, 40`, `20^3` differs from +denser meshes by as much as `2.252 meV/atom`, so it is not converged at the +`1 meV/atom` threshold. Starting from `24^3`, all tested denser even grids stay +within `1 meV/atom`. This is evidence only within that scan; neighbouring odd +meshes and the spacing-derived `75^3` case remain separate checks before a +production mesh is declared. diff --git a/examples/dftb_scc/carbon_bcc_scc_convergence/run.py b/examples/dftb_scc/carbon_bcc_scc_convergence/run.py new file mode 100644 index 000000000..c4ed719c5 --- /dev/null +++ b/examples/dftb_scc/carbon_bcc_scc_convergence/run.py @@ -0,0 +1,470 @@ +#!/usr/bin/env python3 +"""Standalone Carbon BCC self-consistent total-energy k-point scan. + +This example deliberately does not call the multi-structure DFTB+ pytest +benchmark. It runs one fixed Carbon BCC structure and tests convergence of the +self-consistent total energy with respect to k-point density. SCC charge +convergence is recorded as an inner-loop diagnostic, not the primary result. +""" + +from __future__ import annotations + +import argparse +import json +import sys +import time +from dataclasses import asdict, dataclass +from pathlib import Path +from typing import Any, Optional + +import numpy as np +import torch +from ase.io import read + + +EXAMPLE_DIR = Path(__file__).resolve().parent +REPO_ROOT = EXAMPLE_DIR.parents[2] +if str(REPO_ROOT) not in sys.path: + sys.path.insert(0, str(REPO_ROOT)) + +from dptb.kpoints.geometry import calculate_reciprocal_vectors +from dptb.kpoints.mesh import build_kmeshgrid +from dptb.nn.dftb.dftb_scc import SKSCC +from dptb.nn.dftb.scc_params import SCCParams +from dptb.nn.dftb.sk_param import SKParam +from dptb.nn.dftbsk import DFTBSK + + +DEFAULT_STRUCTURE = EXAMPLE_DIR / "POSCAR" +DEFAULT_SK_PATH = REPO_ROOT / "dptb" / "tests" / "data" / "dftb" +BASIS = {"C": ["2s", "2p"]} +VALENCE_ELECTRONS = {"C": 4} +R_MAX = {"C": 4.2} +SIGMA_REP = {"C": 0.561} + + +@dataclass(frozen=True) +class ScanCase: + name: str + meshgrid: Optional[list[int]] + meshspacing: Optional[list[float]] + mixer: str + mix_rate: float + smearing_method: str + temperature_K: float + + +def parse_args() -> argparse.Namespace: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "--scan", + choices=("kmesh", "spacing", "mixer", "smearing"), + default="kmesh", + help="Numerical dimension to scan (default: kmesh).", + ) + parser.add_argument( + "--structure", + type=Path, + default=DEFAULT_STRUCTURE, + help="Carbon BCC structure (default: the POSCAR beside this script).", + ) + parser.add_argument( + "--sk-path", + type=Path, + default=DEFAULT_SK_PATH, + help="Directory containing C-C.skf.", + ) + parser.add_argument( + "--meshes", + type=int, + nargs="+", + default=(8, 12, 16, 20, 24, 28, 32), + help="Isotropic N x N x N grids for --scan kmesh.", + ) + parser.add_argument( + "--spacings", + type=float, + nargs="+", + default=(0.20, 0.10), + help=( + "Isotropic reciprocal-space spacings for --scan spacing. " + "Request the expensive legacy 0.05 case explicitly." + ), + ) + parser.add_argument( + "--fixed-mesh", + type=int, + default=8, + help="Isotropic grid used by mixer and smearing scans.", + ) + parser.add_argument("--mixer", default="pulay", help="Mixer used outside a mixer scan.") + parser.add_argument( + "--mixers", + nargs="+", + default=("simple", "pulay", "broyden", "diis"), + help="Mixers used by --scan mixer.", + ) + parser.add_argument("--mix-rate", type=float, default=0.30) + parser.add_argument( + "--mix-rates", + type=float, + nargs="+", + default=(0.10, 0.30), + help="Mixing rates used by --scan mixer.", + ) + parser.add_argument("--smearing-method", default="FD") + parser.add_argument( + "--smearing-methods", + nargs="+", + default=("FD", "G"), + help="Occupation methods used by --scan smearing.", + ) + parser.add_argument("--temperature", type=float, default=0.1, help="Electronic temperature in K.") + parser.add_argument( + "--temperatures", + type=float, + nargs="+", + default=(0.1, 300.0), + help="Temperatures used by --scan smearing.", + ) + parser.add_argument("--tol", type=float, default=1e-10, help="Maximum charge-residual tolerance.") + parser.add_argument("--max-iter", type=int, default=200) + parser.add_argument( + "--energy-tol-mev-per-atom", + type=float, + default=1.0, + help=( + "Total-energy k-point convergence tolerance in meV/atom. A mesh " + "passes only if it agrees with every denser mesh in the same scan." + ), + ) + parser.add_argument( + "--rotational-symmetry", + action=argparse.BooleanOptionalAction, + default=False, + ) + parser.add_argument( + "--time-inversion-symmetry", + action=argparse.BooleanOptionalAction, + default=True, + ) + parser.add_argument( + "--gamma-centered", + action=argparse.BooleanOptionalAction, + default=True, + ) + parser.add_argument( + "--dry-run", + action="store_true", + help="Resolve meshes and raw k-point counts without running SCC.", + ) + parser.add_argument( + "--output", + type=Path, + help="Optional JSON output path; parent directories are created.", + ) + return parser.parse_args() + + +def validate_args(args: argparse.Namespace) -> None: + if any(mesh < 1 for mesh in (*args.meshes, args.fixed_mesh)): + raise ValueError("All mesh dimensions must be positive integers.") + if any(spacing <= 0 for spacing in args.spacings): + raise ValueError("All mesh spacings must be positive.") + if any(rate <= 0 or rate > 1 for rate in (*args.mix_rates, args.mix_rate)): + raise ValueError("All mixing rates must be in (0, 1].") + if args.temperature <= 0 or any(temp <= 0 for temp in args.temperatures): + raise ValueError("All electronic temperatures must be positive.") + if args.tol <= 0: + raise ValueError("--tol must be positive.") + if args.max_iter < 1: + raise ValueError("--max-iter must be at least 1.") + if args.energy_tol_mev_per_atom <= 0: + raise ValueError("--energy-tol-mev-per-atom must be positive.") + if not args.structure.is_file(): + raise FileNotFoundError(f"Structure not found: {args.structure}") + if not (args.sk_path / "C-C.skf").is_file(): + raise FileNotFoundError(f"C-C.skf not found under: {args.sk_path}") + + +def build_cases(args: argparse.Namespace) -> list[ScanCase]: + fixed_grid = [args.fixed_mesh] * 3 + common = { + "mixer": args.mixer, + "mix_rate": args.mix_rate, + "smearing_method": args.smearing_method, + "temperature_K": args.temperature, + } + + if args.scan == "kmesh": + return [ + ScanCase( + name=f"grid_{mesh}x{mesh}x{mesh}", + meshgrid=[mesh] * 3, + meshspacing=None, + **common, + ) + for mesh in args.meshes + ] + if args.scan == "spacing": + return [ + ScanCase( + name=f"spacing_{spacing:g}", + meshgrid=None, + meshspacing=[spacing] * 3, + **common, + ) + for spacing in args.spacings + ] + if args.scan == "mixer": + return [ + ScanCase( + name=f"{mixer}_rate_{rate:g}", + meshgrid=fixed_grid, + meshspacing=None, + mixer=mixer, + mix_rate=rate, + smearing_method=args.smearing_method, + temperature_K=args.temperature, + ) + for mixer in args.mixers + for rate in args.mix_rates + ] + return [ + ScanCase( + name=f"{method}_T_{temperature:g}K", + meshgrid=fixed_grid, + meshspacing=None, + mixer=args.mixer, + mix_rate=args.mix_rate, + smearing_method=method, + temperature_K=temperature, + ) + for method in args.smearing_methods + for temperature in args.temperatures + ] + + +def resolve_mesh(structure_path: Path, case: ScanCase) -> list[int]: + if case.meshgrid is not None: + return case.meshgrid + structure = read(structure_path) + reciprocal = calculate_reciprocal_vectors(*np.asarray(structure.cell)) + return build_kmeshgrid(*reciprocal, case.meshspacing) + + +def build_solver(sk_path: Path) -> SKSCC: + model = DFTBSK( + basis=BASIS, + skdata=str(sk_path), + overlap=True, + smooth_ski=True, + dtype=torch.float64, + r_max=R_MAX, + ) + sk_params = SKParam( + basis=BASIS, + skdata=str(sk_path), + cal_rcuts=True, + dtype=torch.float64, + ) + return SKSCC( + model=model, + params=SCCParams.from_skparam(sk_params), + overlap=True, + scc_dtype=torch.float64, + ) + + +def scalar(value: Any) -> Optional[float]: + if value is None: + return None + if isinstance(value, torch.Tensor): + return float(value.detach().cpu().reshape(-1)[0]) + return float(value) + + +def run_case( + solver: SKSCC, + structure_path: Path, + case: ScanCase, + args: argparse.Namespace, +) -> dict[str, Any]: + resolved_mesh = resolve_mesh(structure_path, case) + started = time.perf_counter() + total_energy = solver.get_total_energy( + data=str(structure_path), + nel_atom=VALENCE_ELECTRONS, + sigma_rep=SIGMA_REP, + kmeshgrid=case.meshgrid, + kmeshspacing=case.meshspacing, + kgamma_center=args.gamma_centered, + krotational_symmetry=args.rotational_symmetry, + ktime_inversion_symmetry=args.time_inversion_symmetry, + Temp=case.temperature_K, + AtomicData_options={"r_max": R_MAX}, + smearing_method=case.smearing_method, + mixer=case.mixer, + tol=args.tol, + mix_rate=case.mix_rate, + max_iter=args.max_iter, + ) + elapsed = time.perf_counter() - started + + mulliken_charge = np.asarray(solver.mulliken.mul_charge, dtype=float) + delta_charge = np.asarray(solver.mulliken.delta_charge, dtype=float) + target_electrons = sum(solver.mulliken.per_atom_charge) + atom_count = len(mulliken_charge) + return { + **asdict(case), + "resolved_mesh": resolved_mesh, + "raw_kpoints": int(np.prod(resolved_mesh)), + "reduced_kpoints": int(len(solver.mulliken.klist)), + "kpoint_weight_sum": float(np.sum(solver.mulliken.wk)), + "scc_converged": bool(solver.is_converged), + "scc_iterations": int(solver.n_iterations), + "scc_final_charge_residual_e": solver.final_charge_residual, + "scc_charge_residual_history_e": list(solver.charge_residual_history), + "electron_count": float(np.sum(mulliken_charge)), + "electron_count_error": float(np.sum(mulliken_charge) - target_electrons), + "mulliken_charge_e": mulliken_charge.tolist(), + "delta_charge_e": delta_charge.tolist(), + "max_abs_delta_charge_e": float(np.max(np.abs(delta_charge))), + "fermi_level_eV": scalar(solver.E_fermi), + "total_energy_eV": float(total_energy), + "total_energy_per_atom_eV": float(total_energy) / atom_count, + "electronic_energy_eV": scalar(solver.elec_totE), + "repulsive_energy_eV": scalar(solver.total_rep_energy), + "h0_band_energy_eV": scalar(solver.elec_H0_bandE), + "scc_shift_energy_eV": scalar(solver.scc_shift_energy), + "wall_time_s": elapsed, + } + + +def dry_run_case(structure_path: Path, case: ScanCase) -> dict[str, Any]: + resolved_mesh = resolve_mesh(structure_path, case) + return { + **asdict(case), + "resolved_mesh": resolved_mesh, + "raw_kpoints": int(np.prod(resolved_mesh)), + } + + +def add_kpoint_energy_convergence( + results: list[dict[str, Any]], + tolerance_meV_per_atom: float, +) -> None: + completed = [row for row in results if row.get("total_energy_per_atom_eV") is not None] + if not completed: + return + completed.sort(key=lambda row: row["raw_kpoints"]) + reference = completed[-1] + reference_charge = np.asarray(reference["mulliken_charge_e"], dtype=float) + for index, row in enumerate(completed): + row["finest_mesh_case"] = reference["name"] + row["total_energy_delta_to_finest_meV_per_atom"] = 1000.0 * ( + row["total_energy_per_atom_eV"] - reference["total_energy_per_atom_eV"] + ) + if index == 0: + row["successive_total_energy_delta_meV_per_atom"] = None + else: + previous = completed[index - 1] + row["successive_total_energy_delta_meV_per_atom"] = 1000.0 * ( + row["total_energy_per_atom_eV"] - previous["total_energy_per_atom_eV"] + ) + denser = completed[index + 1 :] + if denser: + max_delta = 1000.0 * max( + abs(row["total_energy_per_atom_eV"] - denser_row["total_energy_per_atom_eV"]) + for denser_row in denser + ) + row["max_abs_total_energy_delta_to_denser_meV_per_atom"] = max_delta + row["kpoint_energy_converged_within_scan"] = max_delta <= tolerance_meV_per_atom + else: + row["max_abs_total_energy_delta_to_denser_meV_per_atom"] = None + row["kpoint_energy_converged_within_scan"] = None + row["fermi_delta_to_finest_eV"] = row["fermi_level_eV"] - reference["fermi_level_eV"] + charge = np.asarray(row["mulliken_charge_e"], dtype=float) + row["max_charge_delta_to_finest_e"] = float(np.max(np.abs(charge - reference_charge))) + + +def format_value(value: Any, precision: int = 6) -> str: + if value is None: + return "-" + if isinstance(value, float): + return f"{value:.{precision}g}" + return str(value) + + +def print_results(results: list[dict[str, Any]], dry_run: bool) -> None: + if dry_run: + columns = ("name", "resolved_mesh", "raw_kpoints") + else: + columns = [ + "name", + "resolved_mesh", + "raw_kpoints", + "reduced_kpoints", + "total_energy_per_atom_eV", + "scc_converged", + "scc_iterations", + "scc_final_charge_residual_e", + "electron_count_error", + "wall_time_s", + ] + if any("successive_total_energy_delta_meV_per_atom" in row for row in results): + insert_at = columns.index("scc_converged") + columns[insert_at:insert_at] = [ + "successive_total_energy_delta_meV_per_atom", + "max_abs_total_energy_delta_to_denser_meV_per_atom", + "kpoint_energy_converged_within_scan", + ] + widths = { + column: max(len(column), *(len(format_value(row.get(column))) for row in results)) + for column in columns + } + print(" ".join(column.ljust(widths[column]) for column in columns)) + print(" ".join("-" * widths[column] for column in columns)) + for row in results: + print(" ".join(format_value(row.get(column)).ljust(widths[column]) for column in columns)) + + +def write_output(path: Path, args: argparse.Namespace, results: list[dict[str, Any]]) -> None: + payload = { + "structure": str(args.structure.resolve()), + "sk_path": str(args.sk_path.resolve()), + "scan": args.scan, + "tol": args.tol, + "max_iter": args.max_iter, + "energy_tolerance_meV_per_atom": args.energy_tol_mev_per_atom, + "sigma_rep": SIGMA_REP, + "gamma_centered": args.gamma_centered, + "rotational_symmetry": args.rotational_symmetry, + "time_inversion_symmetry": args.time_inversion_symmetry, + "dry_run": args.dry_run, + "results": results, + } + path.parent.mkdir(parents=True, exist_ok=True) + path.write_text(json.dumps(payload, indent=2) + "\n", encoding="utf-8") + print(f"\nWrote {path}") + + +def main() -> int: + args = parse_args() + validate_args(args) + cases = build_cases(args) + if args.dry_run: + results = [dry_run_case(args.structure, case) for case in cases] + else: + solver = build_solver(args.sk_path) + results = [run_case(solver, args.structure, case, args) for case in cases] + if args.scan in {"kmesh", "spacing"}: + add_kpoint_energy_convergence(results, args.energy_tol_mev_per_atom) + print_results(results, args.dry_run) + if args.output is not None: + write_output(args.output, args, results) + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/examples/silicon/tutorial_v2.2/sktb.json b/examples/silicon/tutorial_v2.2/sktb.json new file mode 100644 index 000000000..04441f34c --- /dev/null +++ b/examples/silicon/tutorial_v2.2/sktb.json @@ -0,0 +1,212 @@ +{ + "version": 2, + "unit": "eV", + "model_options": { + "nnsk": { + "onsite": { + "method": "uniform" + }, + "hopping": { + "method": "poly4pow", + "rs": { + "Si-Si": 5.64 + }, + "w": 0.2 + }, + "soc": {}, + "freeze": false, + "push": false, + "std": 0.01, + "atomic_radius": "cov" + } + }, + "common_options": { + "basis": { + "Si": [ + "3s", + "3p", + "d*" + ] + }, + "dtype": "float32", + "device": "cuda", + "overlap": true + }, + "model_params": { + "onsite": { + "Si-3s-0": [ + 0.22577190399169922 + ], + "Si-3p-0": [ + 0.2590517997741699 + ], + "Si-d*-0": [ + 2.2772068977355957 + ] + }, + "hopping": { + "Si-Si-3s-3s-0": [ + -3.251499891281128, + 2.5374953746795654, + 1.614464521408081, + -7.179530143737793, + 2.173645496368408, + -0.49579325318336487 + ], + "Si-Si-3s-3p-0": [ + -3.975768804550171, + 2.448185920715332, + 2.9865596294403076, + -8.961633682250977, + 2.418846845626831, + -0.0003527404915075749 + ], + "Si-Si-3s-d*-0": [ + -2.621565580368042, + -1.9421700239181519, + 3.0903854370117188, + 3.818775177001953, + -2.722308874130249, + -1.7021923065185547 + ], + "Si-Si-3p-3p-0": [ + 3.6266324520111084, + -0.12833693623542786, + -4.806272506713867, + 5.545502185821533, + -0.4553164541721344, + 6.516263238154352e-05 + ], + "Si-Si-3p-3p-1": [ + -1.412111759185791, + 1.0882034301757812, + 0.5235573053359985, + -2.5653207302093506, + 0.7675366997718811, + 0.8608765602111816 + ], + "Si-Si-3p-d*-0": [ + 1.1860305070877075, + 1.5664552450180054, + -1.7998991012573242, + -2.402371406555176, + 1.47646963596344, + -9.084110206458718e-05 + ], + "Si-Si-3p-d*-1": [ + -1.5894290208816528, + 0.8904605507850647, + 0.7706379890441895, + -1.976307988166809, + 0.41529569029808044, + -7.361805182881653e-05 + ], + "Si-Si-d*-d*-0": [ + 0.605339527130127, + 0.45712438225746155, + -3.7858543395996094, + -5.206111907958984, + 4.897231101989746, + 2.0922365188598633 + ], + "Si-Si-d*-d*-1": [ + 0.32882407307624817, + 1.2024805545806885, + -0.1820884495973587, + -4.604873180389404, + 1.9599689245224, + 0.26321470737457275 + ], + "Si-Si-d*-d*-2": [ + -0.40751418471336365, + -0.2584984302520752, + 0.49840906262397766, + 0.7221986651420593, + -0.5389226078987122, + -2.2478110790252686 + ] + }, + "overlap": { + "Si-Si-3s-3s-0": [ + 0.15159258246421814, + -0.14375975728034973, + -0.04423067346215248, + 0.344686359167099, + -0.11389884352684021, + -0.48649898171424866 + ], + "Si-Si-3s-3p-0": [ + 0.23676219582557678, + -0.19545333087444305, + -0.12815618515014648, + 0.5768850445747375, + -0.17450246214866638, + 0.0001684250310063362 + ], + "Si-Si-3s-d*-0": [ + 0.22751925885677338, + -0.0771898627281189, + -0.2530990242958069, + 0.4345456659793854, + -0.0724181979894638, + 0.00017254143313039094 + ], + "Si-Si-3p-3p-0": [ + -0.31967324018478394, + 0.11857856065034866, + 0.3982611298561096, + -0.7833613157272339, + 0.16310521960258484, + -0.0001644319563638419 + ], + "Si-Si-3p-3p-1": [ + 0.09834299236536026, + -0.10546348989009857, + -0.011958927847445011, + 0.24422256648540497, + -0.09109091013669968, + -0.9448649883270264 + ], + "Si-Si-3p-d*-0": [ + -0.2652897238731384, + -0.29591259360313416, + 0.3893950283527374, + 0.7813657522201538, + -0.4985273480415344, + 1.8512097597122192 + ], + "Si-Si-3p-d*-1": [ + 0.1863427758216858, + -0.16239936649799347, + -0.08167143166065216, + 0.44544517993927, + -0.14129622280597687, + -0.48326730728149414 + ], + "Si-Si-d*-d*-0": [ + 0.17181868851184845, + 0.27958017587661743, + -0.2196175456047058, + -0.8320618867874146, + 0.4289088547229767, + -0.9734773635864258 + ], + "Si-Si-d*-d*-1": [ + -0.2813691794872284, + 0.1467258185148239, + 0.2982117235660553, + -0.7385584712028503, + 0.18330776691436768, + -0.00021080633450765163 + ], + "Si-Si-d*-d*-2": [ + 0.0818948969244957, + -0.09360547363758087, + 0.0026824597734957933, + 0.19942183792591095, + -0.07977117598056793, + 0.9735291004180908 + ] + } + } +} \ No newline at end of file