"""Independent finite-torus HF/KS formulation for ``aiccm2026dev-b``.
The module deliberately does not import :mod:`vibeqc.periodic.ccm`. Its
finite cyclic cluster is the Born--von Karman translation group
``Z/N1 x Z/N2 x Z/N3``. The SCF is evaluated in the group's irreducible
representations (a Gamma-centred k mesh), which is an exact unitary change
of basis from a translation-invariant Gamma-point supercell calculation.
Wigner--Seitz weights below resolve tied minimum-image representatives.
They form a partition of unity for a translation class; they are not
centre-dependent multipliers on ordinary free-space electron-repulsion
integrals. That distinction preserves the permutation symmetries required
by a single RHF energy functional.
"""
from __future__ import annotations
from dataclasses import dataclass
from enum import Enum
from itertools import product
from time import perf_counter
from typing import Sequence
import warnings
import numpy as np
from ._vibeqc_core import (
BasisSet,
PeriodicKSOptions,
PeriodicRHFOptions,
PeriodicSystem,
)
from .kpoints import KPoints
from .pbc_bipole import PBCBipoleRHFResult, run_pbc_bipole_rhf
from .pbc_bipole_rks import PBCBipoleRKSResult, run_pbc_bipole_rks
from .pbc_bipole_uhf import PBCBipoleUHFResult, run_pbc_bipole_uhf
from .pbc_bipole_uks import PBCBipoleUKSResult, run_pbc_bipole_uks
from .periodic_k_gdf import (
PeriodicKRHFGDFResult,
PeriodicKRKSGDFResult,
PeriodicKUHFGDFResult,
run_krhf_periodic_gdf,
run_krks_periodic_gdf,
run_kuhf_periodic_gdf,
run_kuks_periodic_gdf,
)
from .periodic_aiccm2026dev_b_symmetry import (
AICCM2026DevBSymmetryDiagnostics,
AICCM2026DevBSymmetryMode,
AICCM2026DevBSymmetryPlan,
build_aiccm2026dev_b_symmetry_plan,
gamma_matrix_symmetry_residual,
shell_pair_orbits,
shell_quartet_orbits,
)
from .progress import ProgressLogger
__all__ = [
"AICCM2026DevBDiagnostics",
"AICCM2026DevBBackend",
"AICCM2026DevBExperimentalWarning",
"AICCM2026DevBLatticeExtension",
"WignerSeitzRepresentative",
"cyclic_lattice_extension",
"cyclic_gamma_mesh",
"inverse_bloch_transform",
"pair_wigner_seitz_representatives",
"rhf_idempotency_error",
"run_aiccm2026dev_b_rhf",
"run_aiccm2026dev_b_rks",
"run_aiccm2026dev_b_uhf",
"run_aiccm2026dev_b_uks",
"uhf_idempotency_error",
"wigner_seitz_representatives",
]
class AICCM2026DevBExperimentalWarning(UserWarning):
"""The independently derived B stream is not production-certified."""
def _warn_experimental() -> None:
warnings.warn(
"aiccm2026dev-b is experimental: match the dimensional Coulomb gauge, "
"converge the cyclic lattice extension, and inspect the attached "
"invariants before "
"using an energy quantitatively",
category=AICCM2026DevBExperimentalWarning,
stacklevel=3,
)
[docs]
class AICCM2026DevBBackend(str, Enum):
"""Electron-repulsion backend for the same finite-torus functional."""
FOUR_CENTER = "four_center"
RI = "ri"
RIJCOSX = "rijcosx"
@dataclass(frozen=True)
class WignerSeitzRepresentative:
"""One tied minimum-image representative of a cyclic translation.
``residue`` labels the element of the finite translation group.
``translation`` is an integer primitive-cell translation representing
that residue, and all representatives for one residue have equal
``weight = 1 / multiplicity``.
"""
residue: tuple[int, int, int]
translation: tuple[int, int, int]
displacement_bohr: tuple[float, float, float]
weight: float
@dataclass(frozen=True)
class AICCM2026DevBLatticeExtension:
"""Real-space definition of the finite Born--von Karman torus.
``repetitions`` gives the number of primitive translations in each
lattice direction. The cyclic supercell vectors are
``A_i = repetitions[i] * a_i`` and its Wigner--Seitz half-width is
``repetitions[i] / 2`` in primitive-vector coordinates. The reciprocal
character group is therefore exactly the Gamma-centred uniform net with
the same integer tuple; it is a representation of this real-space
object, not an independent approximation parameter.
"""
repetitions: tuple[int, int, int]
wigner_seitz_half_extent: tuple[float, float, float]
supercell_lattice_bohr: np.ndarray
n_cells: int
[docs]
@dataclass(frozen=True)
class AICCM2026DevBDiagnostics:
"""Internal-consistency data attached to an AICCM2026DEV_B result."""
mesh: tuple[int, int, int]
n_cyclic_cells: int
n_kpoints: int
wigner_seitz_partition_error: float
density_idempotency_error: float
electron_count_error: float
inverse_bloch_imaginary_residual: float
backend: str
electronic_method: str
wall_time_seconds: float
lattice_extension: tuple[int, int, int] | None = None
wigner_seitz_half_extent: tuple[float, float, float] | None = None
n_alpha_error: float | None = None
n_beta_error: float | None = None
s_squared: float | None = None
s_squared_ideal: float | None = None
variational_space: str = "translation-invariant closed-shell determinants"
coulomb_gauge: str = "neutral BvK-periodic Coulomb; G=0 removed"
def _normalise_mesh(
system: PeriodicSystem,
mesh: int | Sequence[int],
) -> tuple[int, int, int]:
dim = int(system.dim)
if dim not in (1, 2, 3):
raise ValueError(f"aiccm2026dev-b requires dim=1, 2, or 3; got {dim}")
if isinstance(mesh, (int, np.integer)):
values = [int(mesh)] * dim
else:
values = [int(value) for value in mesh]
if len(values) == dim:
values += [1] * (3 - dim)
elif len(values) != 3:
raise ValueError(
f"aiccm2026dev-b mesh must have length {dim} or 3 for dim={dim}; "
f"got {values!r}"
)
if any(value < 1 for value in values):
raise ValueError(f"aiccm2026dev-b mesh entries must be >= 1; got {values!r}")
if any(values[axis] != 1 for axis in range(dim, 3)):
raise ValueError(
"aiccm2026dev-b inactive lattice directions must have mesh size 1; "
f"got {values!r} for dim={dim}"
)
return tuple(values) # type: ignore[return-value]
def _resolve_lattice_extension(
system: PeriodicSystem,
lattice_extension: int | Sequence[int] | None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
) -> tuple[int, int, int]:
"""Resolve mutually exclusive real- and reciprocal-space controls.
``wigner_seitz_shells=s`` requests ``s`` complete primitive-cell layers
on either side of the central cell, hence the odd cyclic order
``N=2s+1``. ``lattice_extension=N`` also admits even orders; then the
boundary at ``+/-N/2`` is shared and receives the exact tied-image
weights returned by :func:`wigner_seitz_representatives`.
"""
supplied = sum(
value is not None for value in (lattice_extension, mesh, wigner_seitz_shells)
)
if supplied > 1:
raise ValueError(
"choose exactly one of lattice_extension, mesh, or "
"wigner_seitz_shells"
)
if wigner_seitz_shells is not None:
dim = int(system.dim)
if isinstance(wigner_seitz_shells, (int, np.integer)):
shells = [int(wigner_seitz_shells)] * dim
else:
shells = [int(value) for value in wigner_seitz_shells]
if len(shells) == dim:
shells += [0] * (3 - dim)
if len(shells) != 3 or any(value < 0 for value in shells):
raise ValueError("wigner_seitz_shells must contain nonnegative integers")
if any(shells[axis] != 0 for axis in range(dim, 3)):
raise ValueError("inactive directions must have zero Wigner--Seitz shells")
return _normalise_mesh(
system,
[2 * shells[axis] + 1 if axis < dim else 1 for axis in range(3)],
)
selected = mesh if mesh is not None else lattice_extension
if selected is None:
selected = [1] * int(system.dim)
return _normalise_mesh(system, selected)
def cyclic_lattice_extension(
system: PeriodicSystem,
lattice_extension: int | Sequence[int] | None = None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
) -> AICCM2026DevBLatticeExtension:
"""Construct the user-facing real-space cyclic extension descriptor."""
repetitions = _resolve_lattice_extension(
system,
lattice_extension,
mesh=mesh,
wigner_seitz_shells=wigner_seitz_shells,
)
supercell = np.asarray(system.lattice, dtype=float).copy()
for axis, count in enumerate(repetitions):
supercell[axis] *= count
return AICCM2026DevBLatticeExtension(
repetitions=repetitions,
wigner_seitz_half_extent=tuple(0.5 * value for value in repetitions),
supercell_lattice_bohr=supercell,
n_cells=int(np.prod(repetitions)),
)
[docs]
def cyclic_gamma_mesh(
system: PeriodicSystem,
mesh: int | Sequence[int],
) -> KPoints:
"""Return the reciprocal irreps of a finite cyclic cluster.
A cluster with ``N1*N2*N3`` primitive cells has exactly the uniform,
unreduced Gamma-centred ``(N1,N2,N3)`` reciprocal mesh. No classical
even-mesh half shift and no symmetry reduction are allowed here because
either would change the finite translation group.
"""
mesh_tuple = _normalise_mesh(system, mesh)
return KPoints.gamma_centred(system, mesh_tuple, symmetry=False)
def _nearest_representatives(
active_lattice: np.ndarray,
mesh: tuple[int, int, int],
residue: tuple[int, int, int],
dim: int,
tolerance: float,
offset_bohr: np.ndarray,
) -> list[tuple[int, int, int]]:
"""Solve the small closest-vector problem with a proved stopping bound."""
sigma_min = float(np.linalg.svd(active_lattice, compute_uv=False)[-1])
if sigma_min <= 0.0:
raise ValueError("aiccm2026dev-b requires linearly independent active vectors")
best_sq = np.inf
winners: list[tuple[int, int, int]] = []
radius = 0
while True:
shifts = product(range(-radius, radius + 1), repeat=dim)
for shift_active in shifts:
translation = [0, 0, 0]
for axis in range(dim):
translation[axis] = residue[axis] + mesh[axis] * shift_active[axis]
translation_tuple = tuple(translation)
displacement = (
active_lattice @ np.asarray(translation[:dim], dtype=float)
+ offset_bohr
)
norm_sq = float(displacement @ displacement)
scale = max(1.0, best_sq if np.isfinite(best_sq) else 1.0)
if norm_sq < best_sq - tolerance * scale:
best_sq = norm_sq
winners = [translation_tuple]
elif abs(norm_sq - best_sq) <= tolerance * scale:
if translation_tuple not in winners:
winners.append(translation_tuple)
# If an unseen shift has a component outside [-radius, radius], its
# integer coefficient has this lower bound. Multiplication by the
# smallest singular value gives a Cartesian-distance lower bound.
unseen_component = min(
mesh[axis] * (radius + 1) - abs(residue[axis])
for axis in range(dim)
)
unseen_lower = max(
0.0,
sigma_min * max(0, unseen_component) - float(np.linalg.norm(offset_bohr)),
)
unseen_lower_sq = unseen_lower**2
if unseen_lower_sq > best_sq + tolerance * max(1.0, best_sq):
return sorted(winners)
radius += 1
if radius > 64:
raise RuntimeError(
"aiccm2026dev-b Wigner-Seitz closest-vector search did not "
"reach its stopping bound"
)
[docs]
def wigner_seitz_representatives(
system: PeriodicSystem,
mesh: int | Sequence[int],
*,
tolerance: float = 1e-12,
offset_bohr: Sequence[float] = (0.0, 0.0, 0.0),
) -> tuple[WignerSeitzRepresentative, ...]:
"""Construct the exact minimum-image partition of the cyclic group.
``offset_bohr`` is the intra-cell displacement from the centre defining
the Wigner--Seitz region to the translated centre. Boundary points have
more than one equally short representative. Each receives the inverse
multiplicity, so the weights for every group element sum to one. This is
the only Wigner--Seitz weighting used by ``aiccm2026dev-b``.
"""
if tolerance <= 0.0:
raise ValueError("aiccm2026dev-b Wigner-Seitz tolerance must be positive")
mesh_tuple = _normalise_mesh(system, mesh)
dim = int(system.dim)
lattice = np.asarray(system.lattice, dtype=float)
if lattice.shape != (3, 3):
raise ValueError(f"aiccm2026dev-b lattice must have shape (3,3); got {lattice.shape}")
# PeriodicSystem stores lattice vectors as rows. The closest-vector
# problem uses a 3 x dim column generator, hence the transpose here.
active_lattice = lattice[:dim, :].T
offset = np.asarray(offset_bohr, dtype=float)
if offset.shape != (3,) or not np.all(np.isfinite(offset)):
raise ValueError("aiccm2026dev-b Wigner-Seitz offset must be a finite 3-vector")
output: list[WignerSeitzRepresentative] = []
residue_ranges = [range(mesh_tuple[axis]) for axis in range(dim)]
for active_residue in product(*residue_ranges):
residue = tuple(active_residue) + (0,) * (3 - dim)
representatives = _nearest_representatives(
active_lattice,
mesh_tuple,
residue,
dim,
tolerance,
offset,
)
weight = 1.0 / len(representatives)
for translation in representatives:
displacement = np.asarray(translation, dtype=float) @ lattice + offset
output.append(
WignerSeitzRepresentative(
residue=residue,
translation=translation,
displacement_bohr=tuple(float(value) for value in displacement),
weight=weight,
)
)
return tuple(output)
[docs]
def pair_wigner_seitz_representatives(
system: PeriodicSystem,
mesh: int | Sequence[int],
centre_a_bohr: Sequence[float],
centre_b_bohr: Sequence[float],
*,
tolerance: float = 1e-12,
) -> tuple[WignerSeitzRepresentative, ...]:
"""Wigner--Seitz representatives for translations of centre B about A."""
centre_a = np.asarray(centre_a_bohr, dtype=float)
centre_b = np.asarray(centre_b_bohr, dtype=float)
if centre_a.shape != (3,) or centre_b.shape != (3,):
raise ValueError("aiccm2026dev-b pair centres must be Cartesian 3-vectors")
return wigner_seitz_representatives(
system,
mesh,
tolerance=tolerance,
offset_bohr=centre_b - centre_a,
)
def rhf_idempotency_error(
density_k: Sequence[np.ndarray],
overlap_k: Sequence[np.ndarray],
) -> float:
"""Maximum Frobenius residual of ``D(k) S(k) D(k) = 2 D(k)``."""
if len(density_k) != len(overlap_k) or not density_k:
raise ValueError("rhf_idempotency_error requires matching non-empty k blocks")
return max(
float(np.linalg.norm(density @ overlap @ density - 2.0 * density))
for density, overlap in zip(density_k, overlap_k)
)
def uhf_idempotency_error(
density_alpha_k: Sequence[np.ndarray],
density_beta_k: Sequence[np.ndarray],
overlap_k: Sequence[np.ndarray],
) -> float:
"""Maximum spin-projector residual ``P_sigma S P_sigma - P_sigma``."""
if not density_alpha_k or not (
len(density_alpha_k) == len(density_beta_k) == len(overlap_k)
):
raise ValueError("uhf_idempotency_error requires matching non-empty blocks")
return max(
float(np.linalg.norm(density @ overlap @ density - density))
for blocks in (density_alpha_k, density_beta_k)
for density, overlap in zip(blocks, overlap_k)
)
def _electron_count_error(
density_k: Sequence[np.ndarray],
overlap_k: Sequence[np.ndarray],
weights: np.ndarray,
expected: int,
) -> float:
count = sum(
float(weight * np.trace(density @ overlap).real)
for weight, density, overlap in zip(weights, density_k, overlap_k)
)
return abs(count - expected)
def _resolve_backend(
backend: str | AICCM2026DevBBackend,
) -> AICCM2026DevBBackend:
if isinstance(backend, AICCM2026DevBBackend):
return backend
try:
return AICCM2026DevBBackend(str(backend).strip().lower().replace("-", "_"))
except ValueError as exc:
choices = ", ".join(member.value for member in AICCM2026DevBBackend)
raise ValueError(
f"unknown aiccm2026dev-b backend {backend!r}; choose {choices}"
) from exc
def _resolve_symmetry_mode(
mode: str | AICCM2026DevBSymmetryMode,
) -> AICCM2026DevBSymmetryMode:
if isinstance(mode, AICCM2026DevBSymmetryMode):
return mode
try:
return AICCM2026DevBSymmetryMode(str(mode).strip().lower())
except ValueError as exc:
choices = ", ".join(member.value for member in AICCM2026DevBSymmetryMode)
raise ValueError(
f"unknown aiccm2026dev-b symmetry mode {mode!r}; choose {choices}"
) from exc
def _build_symmetry_plan(
system: PeriodicSystem,
mesh: tuple[int, int, int],
mode: str | AICCM2026DevBSymmetryMode,
*,
symmetry_precision: float,
symmetry_require_full_group: bool,
) -> AICCM2026DevBSymmetryPlan | None:
selected = _resolve_symmetry_mode(mode)
if selected is AICCM2026DevBSymmetryMode.OFF:
return None
if selected is AICCM2026DevBSymmetryMode.INTEGRALS:
raise NotImplementedError(
"aiccm2026dev-b symmetry integral acceleration is disabled: "
"general-k AO sewing phases and libint shell-quartet petite-list "
"scattering have not passed the symmetry-off parity gate"
)
return build_aiccm2026dev_b_symmetry_plan(
system,
mesh,
symprec=symmetry_precision,
require_full_space_group=symmetry_require_full_group,
)
def _attach_symmetry_diagnostics(
result: object,
system: PeriodicSystem,
basis: BasisSet,
symmetry_plan: AICCM2026DevBSymmetryPlan | None,
density_k: Sequence[np.ndarray],
) -> None:
if symmetry_plan is None:
return
fock_residual: float | None = None
density_residual: float | None = None
if symmetry_plan.n_kpoints_full == 1:
fock = getattr(result, "fock", None)
if isinstance(fock, (list, tuple)):
fock = fock[0]
if fock is not None:
fock_residual = gamma_matrix_symmetry_residual(
np.asarray(fock), system, basis, symmetry_plan
)
density_residual = gamma_matrix_symmetry_residual(
np.asarray(density_k[0]), system, basis, symmetry_plan
)
pair_orbits = shell_pair_orbits(basis, symmetry_plan)
n_shells = len(list(basis.shells()))
quartet_orbits = (
shell_quartet_orbits(basis, symmetry_plan)
if n_shells <= 24
else None
)
setattr(
result,
"aiccm2026dev_b_symmetry",
AICCM2026DevBSymmetryDiagnostics(
plan=symmetry_plan,
gamma_fock_residual=fock_residual,
gamma_density_residual=density_residual,
n_shell_pairs=sum(len(orbit) for orbit in pair_orbits),
n_unique_shell_pairs=len(pair_orbits),
n_shell_quartets=(
None
if quartet_orbits is None
else sum(len(orbit) for orbit in quartet_orbits)
),
n_unique_shell_quartets=(
None if quartet_orbits is None else len(quartet_orbits)
),
),
)
def _validate_closed_shell_problem(system: PeriodicSystem, options: object) -> None:
if int(system.dim) not in (1, 2, 3):
raise ValueError(
"aiccm2026dev-b requires a 1D, 2D, or 3D periodic system; "
f"got dim={system.dim}"
)
if int(system.charge) != 0:
raise ValueError(
"aiccm2026dev-b requires a neutral primitive cell; charged cells "
"need an explicitly selected background convention"
)
if int(system.multiplicity) != 1:
raise ValueError("aiccm2026dev-b implements closed-shell RHF/RKS only")
if int(system.n_electrons()) % 2:
raise ValueError("aiccm2026dev-b requires an even electron count")
if float(getattr(options, "smearing_temperature", 0.0) or 0.0) != 0.0:
raise ValueError(
"aiccm2026dev-b is a zero-temperature idempotent variational "
"problem; smearing is not implemented"
)
def _validate_open_shell_problem(system: PeriodicSystem, options: object) -> None:
if int(system.dim) not in (1, 2, 3):
raise ValueError(
"aiccm2026dev-b requires a 1D, 2D, or 3D periodic system; "
f"got dim={system.dim}"
)
if int(system.charge) != 0:
raise ValueError(
"aiccm2026dev-b requires a neutral primitive cell; charged cells "
"need an explicitly selected background convention"
)
multiplicity = int(system.multiplicity)
n_electrons = int(system.n_electrons())
if (
multiplicity < 1
or multiplicity > n_electrons + 1
or (n_electrons + multiplicity - 1) % 2
):
raise ValueError(
"aiccm2026dev-b electron count and multiplicity do not define "
"integer alpha/beta occupations"
)
if float(getattr(options, "smearing_temperature", 0.0) or 0.0) != 0.0:
raise ValueError(
"aiccm2026dev-b UHF/UKS currently minimizes idempotent zero-"
"temperature spin projectors; ensemble smearing is not enabled"
)
def _density_blocks_per_k(
result: object,
n_electrons: int,
) -> list[np.ndarray]:
density = getattr(result, "density")
if isinstance(density, (list, tuple)):
return [np.asarray(block) for block in density]
coefficients = list(getattr(result, "mo_coeffs"))
occupations = getattr(result, "occupations", None)
if isinstance(occupations, (list, tuple)) and len(occupations) == len(coefficients):
return [
(np.asarray(coeff) * np.asarray(occ)[None, :]) @ np.asarray(coeff).conj().T
for coeff, occ in zip(coefficients, occupations)
]
n_occ = n_electrons // 2
return [
2.0 * np.asarray(coeff)[:, :n_occ] @ np.asarray(coeff)[:, :n_occ].conj().T
for coeff in coefficients
]
def _spin_density_blocks_per_k(
result: object,
spin: str,
n_occupied: int,
) -> list[np.ndarray]:
density = getattr(result, f"density_{spin}")
if isinstance(density, (list, tuple)):
return [np.asarray(block) for block in density]
coefficients = list(getattr(result, f"mo_coeffs_{spin}"))
occupations = getattr(result, f"occupations_{spin}", None)
if isinstance(occupations, (list, tuple)) and len(occupations) == len(
coefficients
):
return [
(np.asarray(coeff) * np.asarray(occ)[None, :])
@ np.asarray(coeff).conj().T
for coeff, occ in zip(coefficients, occupations)
]
return [
np.asarray(coeff)[:, :n_occupied]
@ np.asarray(coeff)[:, :n_occupied].conj().T
for coeff in coefficients
]
def _attach_diagnostics(
result: object,
system: PeriodicSystem,
basis: BasisSet,
kpoints: KPoints,
mesh: tuple[int, int, int],
backend: AICCM2026DevBBackend,
electronic_method: str,
elapsed: float,
symmetry_plan: AICCM2026DevBSymmetryPlan | None = None,
) -> object:
representatives = wigner_seitz_representatives(system, mesh)
residue_sums: dict[tuple[int, int, int], float] = {}
for representative in representatives:
residue_sums[representative.residue] = (
residue_sums.get(representative.residue, 0.0) + representative.weight
)
partition_error = max(abs(total - 1.0) for total in residue_sums.values())
overlap_k = [np.asarray(block) for block in getattr(result, "overlap")]
is_unrestricted = hasattr(result, "density_alpha")
n_alpha_error = None
n_beta_error = None
if is_unrestricted:
n_electrons = int(system.n_electrons())
two_s = int(system.multiplicity) - 1
n_alpha = (n_electrons + two_s) // 2
n_beta = (n_electrons - two_s) // 2
density_alpha_k = _spin_density_blocks_per_k(result, "alpha", n_alpha)
density_beta_k = _spin_density_blocks_per_k(result, "beta", n_beta)
density_k = [
alpha + beta
for alpha, beta in zip(density_alpha_k, density_beta_k)
]
idempotency_error = uhf_idempotency_error(
density_alpha_k, density_beta_k, overlap_k
)
n_alpha_error = _electron_count_error(
density_alpha_k,
overlap_k,
np.asarray(kpoints.weights, dtype=float),
n_alpha,
)
n_beta_error = _electron_count_error(
density_beta_k,
overlap_k,
np.asarray(kpoints.weights, dtype=float),
n_beta,
)
transformed = [
inverse_bloch_transform(
blocks,
kpoints.kpoints_frac,
[representative.residue for representative in representatives],
kpoints.weights,
)
for blocks in (density_alpha_k, density_beta_k)
]
imaginary_residual = max(
float(np.max(np.abs(blocks.imag))) for blocks in transformed
)
variational_space = "translation-invariant unrestricted determinants"
else:
density_k = _density_blocks_per_k(result, int(system.n_electrons()))
idempotency_error = rhf_idempotency_error(density_k, overlap_k)
density_blocks = inverse_bloch_transform(
density_k,
kpoints.kpoints_frac,
[representative.residue for representative in representatives],
kpoints.weights,
)
imaginary_residual = float(np.max(np.abs(density_blocks.imag)))
variational_space = "translation-invariant closed-shell determinants"
extension = cyclic_lattice_extension(system, mesh)
diagnostics = AICCM2026DevBDiagnostics(
mesh=mesh,
n_cyclic_cells=int(np.prod(mesh)),
n_kpoints=len(kpoints.weights),
wigner_seitz_partition_error=partition_error,
density_idempotency_error=idempotency_error,
electron_count_error=_electron_count_error(
density_k,
overlap_k,
np.asarray(kpoints.weights, dtype=float),
int(system.n_electrons()),
),
inverse_bloch_imaginary_residual=imaginary_residual,
backend=backend.value,
electronic_method=electronic_method,
wall_time_seconds=elapsed,
lattice_extension=extension.repetitions,
wigner_seitz_half_extent=extension.wigner_seitz_half_extent,
n_alpha_error=n_alpha_error,
n_beta_error=n_beta_error,
s_squared=(
float(getattr(result, "s_squared"))
if hasattr(result, "s_squared")
else None
),
s_squared_ideal=(
float(getattr(result, "s_squared_ideal"))
if hasattr(result, "s_squared_ideal")
else None
),
variational_space=variational_space,
)
setattr(result, "backend", f"aiccm2026dev-b-{backend.value}")
setattr(result, "aiccm2026dev_b", diagnostics)
_attach_symmetry_diagnostics(
result,
system,
basis,
symmetry_plan,
density_k,
)
return result
[docs]
def run_aiccm2026dev_b_rhf(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
backend: str | AICCM2026DevBBackend = AICCM2026DevBBackend.FOUR_CENTER,
aux_basis: str | None = None,
gdf_method: str = "rsgdf",
rsgdf_ke_cutoff: float = 200.0,
mdf_ke_cutoff: float = 40.0,
fock_mixing: float | None = None,
symmetry_mode: str | AICCM2026DevBSymmetryMode = (
AICCM2026DevBSymmetryMode.OFF
),
symmetry_precision: float = 1.0e-5,
symmetry_require_full_group: bool = False,
progress: bool | ProgressLogger | None = None,
verbose: int | None = None,
) -> PeriodicKRHFGDFResult | PBCBipoleRHFResult:
"""Minimise the finite-torus closed-shell RHF energy per cell."""
_warn_experimental()
mesh_tuple = _resolve_lattice_extension(
system,
lattice_extension,
mesh=mesh,
wigner_seitz_shells=wigner_seitz_shells,
)
selected = _resolve_backend(backend)
opts = options if options is not None else PeriodicRHFOptions()
_validate_closed_shell_problem(system, opts)
if (
selected is not AICCM2026DevBBackend.FOUR_CENTER
and gdf_method not in ("rsgdf", "mdf")
):
raise ValueError(
"aiccm2026dev-b requires gdf_method='rsgdf' or 'mdf'; the q-only "
"compcell fit is not a consistent finite-torus Hamiltonian"
)
kpoints = cyclic_gamma_mesh(system, mesh_tuple)
symmetry_plan = _build_symmetry_plan(
system,
mesh_tuple,
symmetry_mode,
symmetry_precision=symmetry_precision,
symmetry_require_full_group=symmetry_require_full_group,
)
if selected is AICCM2026DevBBackend.RIJCOSX and int(np.prod(mesh_tuple)) == 1:
raise NotImplementedError(
"aiccm2026dev-b RIJCOSX requires at least two cyclic cells; the "
"native multi-k COSX bridge has no Gamma-only implementation"
)
started = perf_counter()
if selected is AICCM2026DevBBackend.FOUR_CENTER:
use_3d_ewald = int(system.dim) == 3
result = run_pbc_bipole_rhf(
system,
basis,
kpoints.to_bloch_kmesh(),
opts,
use_ewald_j_split=use_3d_ewald,
use_exchange_ewald_split=use_3d_ewald,
exchange_exxdiv=("ewald" if use_3d_ewald else "none"),
progress=progress,
verbose=verbose,
)
else:
result = run_krhf_periodic_gdf(
system,
basis,
kpoints,
opts,
functional=None,
aux_basis=aux_basis,
fock_mixing=fock_mixing,
use_compcell=selected is AICCM2026DevBBackend.RIJCOSX,
k_exchange=(
"cosx" if selected is AICCM2026DevBBackend.RIJCOSX else "gdf"
),
gdf_method=gdf_method,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
check_energy_sanity=True,
progress=progress,
verbose=verbose,
)
return _attach_diagnostics(
result,
system,
basis,
kpoints,
mesh_tuple,
selected,
"RHF",
perf_counter() - started,
symmetry_plan,
)
[docs]
def run_aiccm2026dev_b_rks(
system: PeriodicSystem,
basis: BasisSet,
functional: str,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicKSOptions | None = None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
backend: str | AICCM2026DevBBackend = AICCM2026DevBBackend.FOUR_CENTER,
aux_basis: str | None = None,
gdf_method: str = "rsgdf",
rsgdf_ke_cutoff: float = 200.0,
mdf_ke_cutoff: float = 40.0,
fock_mixing: float | None = None,
symmetry_mode: str | AICCM2026DevBSymmetryMode = (
AICCM2026DevBSymmetryMode.OFF
),
symmetry_precision: float = 1.0e-5,
symmetry_require_full_group: bool = False,
progress: bool | ProgressLogger | None = None,
verbose: int | None = None,
) -> PeriodicKRKSGDFResult | PBCBipoleRKSResult:
"""Minimise the finite-torus closed-shell Kohn--Sham energy per cell."""
_warn_experimental()
if not str(functional).strip():
raise ValueError("aiccm2026dev-b RKS requires a functional name")
mesh_tuple = _resolve_lattice_extension(
system,
lattice_extension,
mesh=mesh,
wigner_seitz_shells=wigner_seitz_shells,
)
selected = _resolve_backend(backend)
opts = options if options is not None else PeriodicKSOptions()
opts.functional = str(functional)
_validate_closed_shell_problem(system, opts)
if (
selected is not AICCM2026DevBBackend.FOUR_CENTER
and gdf_method not in ("rsgdf", "mdf")
):
raise ValueError(
"aiccm2026dev-b requires gdf_method='rsgdf' or 'mdf'; the q-only "
"compcell fit is not a consistent finite-torus Hamiltonian"
)
kpoints = cyclic_gamma_mesh(system, mesh_tuple)
symmetry_plan = _build_symmetry_plan(
system,
mesh_tuple,
symmetry_mode,
symmetry_precision=symmetry_precision,
symmetry_require_full_group=symmetry_require_full_group,
)
if (
selected is not AICCM2026DevBBackend.FOUR_CENTER
and int(np.prod(mesh_tuple)) == 1
):
raise NotImplementedError(
"aiccm2026dev-b RI and RIJCOSX RKS require at least two cyclic "
"cells because the native Gamma RKS GDF path is not pair-resolved"
)
started = perf_counter()
if selected is AICCM2026DevBBackend.FOUR_CENTER:
use_3d_ewald = int(system.dim) == 3
result = run_pbc_bipole_rks(
system,
basis,
kpoints.to_bloch_kmesh(),
opts,
functional=str(functional),
use_ewald_j_split=use_3d_ewald,
use_exchange_ewald_split=use_3d_ewald,
exchange_exxdiv=("ewald" if use_3d_ewald else "none"),
progress=progress,
verbose=verbose,
)
else:
result = run_krks_periodic_gdf(
system,
basis,
kpoints,
opts,
functional=str(functional),
aux_basis=aux_basis,
fock_mixing=fock_mixing,
use_compcell=True,
k_exchange=(
"cosx" if selected is AICCM2026DevBBackend.RIJCOSX else "gdf"
),
gdf_method=gdf_method,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
check_energy_sanity=True,
progress=progress,
verbose=verbose,
)
return _attach_diagnostics(
result,
system,
basis,
kpoints,
mesh_tuple,
selected,
f"RKS/{functional}",
perf_counter() - started,
symmetry_plan,
)
[docs]
def run_aiccm2026dev_b_uhf(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
backend: str | AICCM2026DevBBackend = AICCM2026DevBBackend.FOUR_CENTER,
aux_basis: str | None = None,
gdf_method: str = "rsgdf",
rsgdf_ke_cutoff: float = 200.0,
mdf_ke_cutoff: float = 40.0,
fock_mixing: float | None = None,
symmetry_mode: str | AICCM2026DevBSymmetryMode = (
AICCM2026DevBSymmetryMode.OFF
),
symmetry_precision: float = 1.0e-5,
symmetry_require_full_group: bool = False,
progress: bool | ProgressLogger | None = None,
verbose: int | None = None,
) -> PeriodicKUHFGDFResult | PBCBipoleUHFResult:
"""Minimize the finite-torus unrestricted Hartree--Fock functional."""
_warn_experimental()
extension = _resolve_lattice_extension(
system,
lattice_extension,
mesh=mesh,
wigner_seitz_shells=wigner_seitz_shells,
)
selected = _resolve_backend(backend)
opts = options if options is not None else PeriodicRHFOptions()
_validate_open_shell_problem(system, opts)
if fock_mixing not in (None, 0.0):
raise NotImplementedError(
"aiccm2026dev-b UHF does not silently substitute Fock mixing; "
"use the shared DIIS options and diagnose oscillatory states"
)
if selected is not AICCM2026DevBBackend.FOUR_CENTER and gdf_method not in (
"rsgdf",
"mdf",
):
raise ValueError("aiccm2026dev-b UHF requires gdf_method='rsgdf' or 'mdf'")
if selected is AICCM2026DevBBackend.RIJCOSX and int(np.prod(extension)) == 1:
raise NotImplementedError(
"aiccm2026dev-b RIJCOSX requires at least two cyclic cells"
)
kpoints = cyclic_gamma_mesh(system, extension)
symmetry_plan = _build_symmetry_plan(
system,
extension,
symmetry_mode,
symmetry_precision=symmetry_precision,
symmetry_require_full_group=symmetry_require_full_group,
)
started = perf_counter()
if selected is AICCM2026DevBBackend.FOUR_CENTER:
use_3d_ewald = int(system.dim) == 3
result = run_pbc_bipole_uhf(
system,
basis,
kpoints.to_bloch_kmesh(),
opts,
use_ewald_j_split=use_3d_ewald,
use_exchange_ewald_split=use_3d_ewald,
exchange_exxdiv=("ewald" if use_3d_ewald else "none"),
progress=progress,
verbose=verbose,
)
else:
result = run_kuhf_periodic_gdf(
system,
basis,
kpoints,
opts,
functional=None,
aux_basis=aux_basis,
gdf_method=gdf_method,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
k_exchange=(
"cosx" if selected is AICCM2026DevBBackend.RIJCOSX else "gdf"
),
check_energy_sanity=True,
progress=progress,
verbose=verbose,
)
return _attach_diagnostics(
result,
system,
basis,
kpoints,
extension,
selected,
"UHF",
perf_counter() - started,
symmetry_plan,
)
[docs]
def run_aiccm2026dev_b_uks(
system: PeriodicSystem,
basis: BasisSet,
functional: str,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicKSOptions | None = None,
*,
mesh: int | Sequence[int] | None = None,
wigner_seitz_shells: int | Sequence[int] | None = None,
backend: str | AICCM2026DevBBackend = AICCM2026DevBBackend.FOUR_CENTER,
aux_basis: str | None = None,
gdf_method: str = "rsgdf",
rsgdf_ke_cutoff: float = 200.0,
mdf_ke_cutoff: float = 40.0,
fock_mixing: float | None = None,
symmetry_mode: str | AICCM2026DevBSymmetryMode = (
AICCM2026DevBSymmetryMode.OFF
),
symmetry_precision: float = 1.0e-5,
symmetry_require_full_group: bool = False,
progress: bool | ProgressLogger | None = None,
verbose: int | None = None,
) -> PeriodicKUHFGDFResult | PBCBipoleUKSResult:
"""Minimize the finite-torus spin-polarized Kohn--Sham functional."""
_warn_experimental()
if not str(functional).strip():
raise ValueError("aiccm2026dev-b UKS requires a functional name")
extension = _resolve_lattice_extension(
system,
lattice_extension,
mesh=mesh,
wigner_seitz_shells=wigner_seitz_shells,
)
selected = _resolve_backend(backend)
opts = options if options is not None else PeriodicKSOptions()
opts.functional = str(functional)
_validate_open_shell_problem(system, opts)
if fock_mixing not in (None, 0.0):
raise NotImplementedError(
"aiccm2026dev-b UKS does not silently substitute Fock mixing; "
"use the shared DIIS options and diagnose oscillatory states"
)
if selected is not AICCM2026DevBBackend.FOUR_CENTER and gdf_method not in (
"rsgdf",
"mdf",
):
raise ValueError("aiccm2026dev-b UKS requires gdf_method='rsgdf' or 'mdf'")
if selected is AICCM2026DevBBackend.RIJCOSX and int(np.prod(extension)) == 1:
raise NotImplementedError(
"aiccm2026dev-b RIJCOSX requires at least two cyclic cells"
)
kpoints = cyclic_gamma_mesh(system, extension)
symmetry_plan = _build_symmetry_plan(
system,
extension,
symmetry_mode,
symmetry_precision=symmetry_precision,
symmetry_require_full_group=symmetry_require_full_group,
)
started = perf_counter()
if selected is AICCM2026DevBBackend.FOUR_CENTER:
use_3d_ewald = int(system.dim) == 3
result = run_pbc_bipole_uks(
system,
basis,
kpoints.to_bloch_kmesh(),
opts,
functional=str(functional),
use_ewald_j_split=use_3d_ewald,
use_exchange_ewald_split=use_3d_ewald,
exchange_exxdiv=("ewald" if use_3d_ewald else "none"),
progress=progress,
verbose=verbose,
)
else:
result = run_kuks_periodic_gdf(
system,
basis,
kpoints,
opts,
functional=str(functional),
aux_basis=aux_basis,
gdf_method=gdf_method,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
mdf_ke_cutoff=mdf_ke_cutoff,
k_exchange=(
"cosx" if selected is AICCM2026DevBBackend.RIJCOSX else "gdf"
),
check_energy_sanity=True,
progress=progress,
verbose=verbose,
)
return _attach_diagnostics(
result,
system,
basis,
kpoints,
extension,
selected,
f"UKS/{functional}",
perf_counter() - started,
symmetry_plan,
)