"""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 for
the same already specified χ-CCM finite Hamiltonian. This internal
representation theorem does not identify χ-CCM with the separately
constructed union-and-weight Γ-CCM approach.
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, field
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,
Functional,
PeriodicKSOptions,
PeriodicRHFOptions,
PeriodicSystem,
)
from .kpoints import KPoints
from .pbc_bipole import PBCBipoleRHFResult, run_pbc_bipole_rhf
from .pbc_bipole_common import resolve_fock_mixing
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
from .periodic_screened_exchange import resolve_periodic_exchange
__all__ = [
"AICCM2026DevBChargeBookkeeping",
"AICCM2026DevBDiagnostics",
"AICCM2026DevBBackend",
"AICCM2026DevBExperimentalWarning",
"AICCM2026DevBFiniteTorusConvention",
"AICCM2026DevBLatticeExtension",
"WignerSeitzRepresentative",
"aiccm2026dev_b_charge_bookkeeping",
"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"
_CCM_APPROACH = "chi-ccm"
_CCM_CONSTRUCTION = "finite-translation-group-character"
_EVALUATION_REPRESENTATION = "gamma-centred-character-mesh"
[docs]
@dataclass(frozen=True)
class AICCM2026DevBFiniteTorusConvention:
"""Finite-N Hamiltonian convention attached to every B-stream result.
The descriptor is intentionally explicit about the exchange ``q=0`` seam:
finite-N HF, MP2, CCSD(T), and local-correlation numbers are comparable
only when this convention matches. A strict-zero-mode exchange reference
is a different finite torus, even though it has the same thermodynamic
target.
"""
coulomb_kernel: str
exchange_q0: str
exchange_q0_applicability: str
boundary_model: str
periodic_dimension: int
character_mesh_shape: tuple[int, int, int]
bvk_madelung_supercell_repetitions: tuple[int, int, int]
bvk_madelung_supercell_lattice_bohr: tuple[tuple[float, float, float], ...]
lattice_vector_convention: str = "columns"
ccm_approach: str = _CCM_APPROACH
ccm_construction: str = _CCM_CONSTRUCTION
evaluation_representation: str = _EVALUATION_REPRESENTATION
orbital_energy_convention: str = field(init=False)
def __post_init__(self) -> None:
identity = {
"ccm_approach": (self.ccm_approach, _CCM_APPROACH),
"ccm_construction": (self.ccm_construction, _CCM_CONSTRUCTION),
"evaluation_representation": (
self.evaluation_representation,
_EVALUATION_REPRESENTATION,
),
}
for field_name, (actual, expected) in identity.items():
if actual != expected:
raise ValueError(
f"{field_name} must be {expected!r}; got {actual!r}"
)
if self.lattice_vector_convention != "columns":
raise ValueError(
"lattice_vector_convention must be 'columns'; "
f"got {self.lattice_vector_convention!r}"
)
if self.exchange_q0_applicability == "active":
orbital_energy_convention = (
"HF/Kohn-Sham eigenvalues include the declared exchange_q0 seam"
)
elif self.exchange_q0_applicability == "inactive":
orbital_energy_convention = (
"Kohn-Sham eigenvalues do not include the declared exchange_q0 "
"seam because the full-range exact-exchange coefficient is zero"
)
else:
raise ValueError(
"exchange_q0_applicability must be 'active' or 'inactive'; "
f"got {self.exchange_q0_applicability!r}"
)
object.__setattr__(
self,
"orbital_energy_convention",
orbital_energy_convention,
)
@dataclass(frozen=True)
class AICCM2026DevBChargeBookkeeping:
"""Effective charge/electron convention used by a χ-CCM SCF run.
``physical_electrons`` is the full all-electron primitive-cell count from
:class:`PeriodicSystem`. When a periodic ECP is attached to the SCF
options, ``ecp_total_ncore`` electrons are removed from the variational
space and the Coulomb background must use the corresponding
``effective_nuclear_charges``. The effective net charge is therefore
``sum(effective_nuclear_charges) - effective_electrons``; production
χ-CCM routes require it to vanish.
"""
physical_electrons: int
effective_electrons: int
ecp_total_ncore: int
effective_nuclear_charges: tuple[float, ...]
effective_nuclear_charge: float
effective_net_charge: float
has_ecp: bool
@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
scf_trace_length: int = 0
final_scf_delta_e_ha: float | None = None
final_scf_grad_norm: float | None = None
final_scf_diis_subspace: int | None = None
physical_electron_count: int | None = None
effective_electron_count: int | None = None
ecp_total_ncore: int = 0
effective_nuclear_charge: float | None = None
effective_nuclear_charges: tuple[float, ...] | None = None
effective_net_charge: float | None = None
use_diis: bool | None = None
diis_start_iter: int | None = None
diis_subspace_size: int | None = None
scf_accelerator: str | None = None
damping: float | None = None
dynamic_damping: bool | None = None
# Executed effective previous-Fock weight, after backend defaults.
fock_mixing: float | None = None
level_shift: float | None = None
level_shift_warmup_cycles: int | None = None
smearing_temperature: float | None = None
variational_space: str = "translation-invariant closed-shell determinants"
coulomb_gauge: str = "neutral BvK-periodic Coulomb; G=0 removed"
finite_torus_convention: AICCM2026DevBFiniteTorusConvention | None = None
@property
def coulomb_kernel(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.coulomb_kernel
@property
def exchange_q0(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.exchange_q0
@property
def exchange_q0_applicability(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.exchange_q0_applicability
@property
def boundary_model(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.boundary_model
@property
def ccm_approach(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.ccm_approach
@property
def ccm_construction(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.ccm_construction
@property
def evaluation_representation(self) -> str:
if self.finite_torus_convention is None:
return "not-recorded"
return self.finite_torus_convention.evaluation_representation
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,
)
lattice = np.asarray(system.lattice, dtype=float)
supercell = lattice @ np.diag(np.asarray(repetitions, dtype=float))
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)),
)
def _finite_torus_convention(
system: PeriodicSystem,
mesh: tuple[int, int, int],
extension: AICCM2026DevBLatticeExtension | None = None,
*,
full_range_exchange_coefficient: float,
) -> AICCM2026DevBFiniteTorusConvention:
"""Return the declared finite-N Hamiltonian convention for B production."""
dim = int(system.dim)
if extension is None:
extension = cyclic_lattice_extension(system, mesh)
c_full = float(full_range_exchange_coefficient)
if not np.isfinite(c_full):
raise ValueError(
"full_range_exchange_coefficient must be finite; "
f"got {full_range_exchange_coefficient!r}"
)
exchange_q0_applicability = "active" if c_full != 0.0 else "inactive"
boundary_model = "3d-periodic" if dim == 3 else "3d-periodic-vacuum"
lattice_matrix = tuple(
tuple(float(value) for value in row)
for row in np.asarray(extension.supercell_lattice_bohr, dtype=float)
)
return AICCM2026DevBFiniteTorusConvention(
coulomb_kernel="3d-periodic-g0",
exchange_q0="bvk-ewald",
exchange_q0_applicability=exchange_q0_applicability,
boundary_model=boundary_model,
periodic_dimension=dim,
character_mesh_shape=mesh,
bvk_madelung_supercell_repetitions=extension.repetitions,
bvk_madelung_supercell_lattice_bohr=lattice_matrix,
)
def _full_range_exchange_coefficient(functional: str, *, spin: int) -> float:
"""Return the live full-range exact-exchange coefficient for KS metadata."""
func = Functional(str(functional), int(spin))
assembly = resolve_periodic_exchange(
func,
where="aiccm2026dev-b finite-torus convention",
)
return float(assembly.c_full)
[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 columns. Restrict the
# closest-vector generator to the active primitive directions.
active_lattice = lattice[:, :dim]
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 = lattice @ np.asarray(translation, dtype=float) + 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 _inverse_bloch_transform_native(
matrices: np.ndarray,
kfrac: np.ndarray,
translations: np.ndarray,
weights: np.ndarray,
) -> np.ndarray | None:
"""Return native inverse Bloch blocks, or ``None`` for older extensions."""
try:
from . import _vibeqc_core
except ImportError:
return None
kernel = getattr(_vibeqc_core, "aiccm2026dev_b_inverse_bloch_transform", None)
if kernel is None:
return None
return np.asarray(
kernel(
np.ascontiguousarray(matrices, dtype=np.complex128),
np.ascontiguousarray(kfrac, dtype=float),
np.ascontiguousarray(translations, dtype=np.int64),
np.ascontiguousarray(weights, dtype=float),
),
dtype=np.complex128,
)
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 _require_supported_scf_dimension(system: PeriodicSystem) -> None:
"""Reject dimensional Coulomb gauges that every current SCF path lacks.
The 3D BIPOLE path has a neutral Ewald split. The current 1D/2D
fallback is a direct-truncated interaction, so its absolute energy is
not the declared finite-torus Hamiltonian. The lower-dimensional fitted
mesh is also unavailable because collapsing every transverse reciprocal
component does not define a Coulomb kernel. Returning a converged number
from either family is worse than failing.
"""
if int(system.dim) == 3:
return
raise NotImplementedError(
"aiccm2026dev-b SCF is currently enabled only for 3D periodicity. "
"The 1D/2D four-centre path lacks a neutral wire/slab "
"Coulomb gauge and was found to over-bind chain benchmarks by a "
"Madelung-scale constant shift. RI and RIJCOSX are not alternatives: "
"their shared lower-dimensional mesh collapses transverse reciprocal "
"structure and is not a Coulomb kernel. All backends remain blocked "
"until one neutral wire/slab Hamiltonian is implemented."
)
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 _option_sequence(options: object | None, name: str) -> tuple[object, ...]:
if options is None or not hasattr(options, name):
return ()
value = getattr(options, name)
if value is None:
return ()
try:
return tuple(value)
except TypeError:
return ()
def aiccm2026dev_b_charge_bookkeeping(
system: PeriodicSystem,
options: object | None = None,
) -> AICCM2026DevBChargeBookkeeping:
"""Return the effective charge/electron convention for a χ-CCM run."""
atoms = list(system.unit_cell)
physical_electrons = int(system.n_electrons())
bare_charges = tuple(float(atom.Z) for atom in atoms)
ecp_total_ncore = int(getattr(options, "ecp_total_ncore", 0) or 0)
effective_charges_raw = _option_sequence(options, "ecp_effective_charges")
ecp_blocks = _option_sequence(options, "ecp_primitive_blocks")
ecp_centers = _option_sequence(options, "ecp_home_centers")
has_ecp = bool(effective_charges_raw or ecp_blocks or ecp_centers or ecp_total_ncore)
if ecp_total_ncore < 0:
raise ValueError("aiccm2026dev-b ECP core electron count must be non-negative")
if ecp_total_ncore > physical_electrons:
raise ValueError(
"aiccm2026dev-b ECP core electron count exceeds the physical "
f"electron count ({ecp_total_ncore} > {physical_electrons})"
)
if has_ecp:
if len(effective_charges_raw) != len(atoms):
raise ValueError(
"aiccm2026dev-b ECP runs require one effective nuclear charge "
"per primitive-cell atom"
)
effective_charges = tuple(float(charge) for charge in effective_charges_raw)
expected_charge = float(sum(bare_charges) - ecp_total_ncore)
observed_charge = float(sum(effective_charges))
if abs(observed_charge - expected_charge) > 1.0e-8:
raise ValueError(
"aiccm2026dev-b ECP effective nuclear charges are inconsistent "
"with ecp_total_ncore: "
f"sum(Z_eff)={observed_charge:.12g}, expected {expected_charge:.12g}"
)
else:
effective_charges = bare_charges
observed_charge = float(sum(effective_charges))
effective_electrons = physical_electrons - ecp_total_ncore
net_charge = observed_charge - float(effective_electrons)
return AICCM2026DevBChargeBookkeeping(
physical_electrons=physical_electrons,
effective_electrons=effective_electrons,
ecp_total_ncore=ecp_total_ncore,
effective_nuclear_charges=effective_charges,
effective_nuclear_charge=observed_charge,
effective_net_charge=net_charge,
has_ecp=has_ecp,
)
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}"
)
bookkeeping = aiccm2026dev_b_charge_bookkeeping(system, options)
if abs(bookkeeping.effective_net_charge) > 1.0e-8:
raise ValueError(
"aiccm2026dev-b requires a neutral effective primitive cell; "
f"effective net charge is {bookkeeping.effective_net_charge:.12g}. "
"Charged cells need an explicitly selected background convention"
)
if int(system.multiplicity) != 1:
raise ValueError("aiccm2026dev-b implements closed-shell RHF/RKS only")
if bookkeeping.effective_electrons % 2:
raise ValueError(
"aiccm2026dev-b requires an even effective electron count; "
f"got {bookkeeping.effective_electrons}"
)
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}"
)
bookkeeping = aiccm2026dev_b_charge_bookkeeping(system, options)
if abs(bookkeeping.effective_net_charge) > 1.0e-8:
raise ValueError(
"aiccm2026dev-b requires a neutral effective primitive cell; "
f"effective net charge is {bookkeeping.effective_net_charge:.12g}. "
"Charged cells need an explicitly selected background convention"
)
multiplicity = int(system.multiplicity)
n_electrons = bookkeeping.effective_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 _trace_field(item: object, name: str, index: int) -> object | None:
if hasattr(item, name):
return getattr(item, name)
if isinstance(item, (tuple, list)) and len(item) > index:
return item[index]
return None
def _final_scf_trace_values(
result: object,
) -> tuple[int, float | None, float | None, int | None]:
trace = list(getattr(result, "scf_trace", []) or [])
if not trace:
return 0, None, None, None
last = trace[-1]
delta_e = _trace_field(last, "delta_e", 2)
grad_norm = _trace_field(last, "grad_norm", 3)
diis_subspace = _trace_field(last, "diis_subspace", 4)
return (
len(trace),
None if delta_e is None else float(delta_e),
None if grad_norm is None else float(grad_norm),
None if diis_subspace is None else int(diis_subspace),
)
def _scf_accelerator_name(options: object | None) -> str | None:
if options is None or not hasattr(options, "scf_accelerator"):
return None
accelerator = getattr(options, "scf_accelerator")
name = getattr(accelerator, "name", None)
if isinstance(name, str):
return name
return str(accelerator).split(".")[-1]
def _executed_fock_mixing(
result: object,
options: object | None,
) -> float | None:
"""Return the backend-resolved previous-Fock weight when available."""
executed = getattr(result, "fock_mixing", None)
if executed is not None:
return float(executed)
if options is None:
return None
return float(getattr(options, "fock_mixing", 0.0))
def _attach_diagnostics(
result: object,
system: PeriodicSystem,
basis: BasisSet,
kpoints: KPoints,
mesh: tuple[int, int, int],
backend: AICCM2026DevBBackend,
electronic_method: str,
full_range_exchange_coefficient: float,
elapsed: float,
symmetry_plan: AICCM2026DevBSymmetryPlan | None = None,
options: object | 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")]
charge_bookkeeping = aiccm2026dev_b_charge_bookkeeping(system, options)
effective_electrons = charge_bookkeeping.effective_electrons
is_unrestricted = hasattr(result, "density_alpha")
n_alpha_error = None
n_beta_error = None
if is_unrestricted:
two_s = int(system.multiplicity) - 1
n_alpha = (effective_electrons + two_s) // 2
n_beta = (effective_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, effective_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)
convention = _finite_torus_convention(
system,
mesh,
extension,
full_range_exchange_coefficient=full_range_exchange_coefficient,
)
(
scf_trace_length,
final_scf_delta_e,
final_scf_grad_norm,
final_scf_diis_subspace,
) = _final_scf_trace_values(result)
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),
effective_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
),
scf_trace_length=scf_trace_length,
final_scf_delta_e_ha=final_scf_delta_e,
final_scf_grad_norm=final_scf_grad_norm,
final_scf_diis_subspace=final_scf_diis_subspace,
physical_electron_count=charge_bookkeeping.physical_electrons,
effective_electron_count=effective_electrons,
ecp_total_ncore=charge_bookkeeping.ecp_total_ncore,
effective_nuclear_charge=charge_bookkeeping.effective_nuclear_charge,
effective_nuclear_charges=charge_bookkeeping.effective_nuclear_charges,
effective_net_charge=charge_bookkeeping.effective_net_charge,
use_diis=(
None if options is None else bool(getattr(options, "use_diis", False))
),
diis_start_iter=(
None if options is None else int(getattr(options, "diis_start_iter", 0))
),
diis_subspace_size=(
None if options is None else int(getattr(options, "diis_subspace_size", 0))
),
scf_accelerator=_scf_accelerator_name(options),
damping=(
None if options is None else float(getattr(options, "damping", 0.0))
),
dynamic_damping=(
None
if options is None
else bool(getattr(options, "dynamic_damping", False))
),
fock_mixing=_executed_fock_mixing(result, options),
level_shift=(
None if options is None else float(getattr(options, "level_shift", 0.0))
),
level_shift_warmup_cycles=(
None
if options is None
else int(getattr(options, "level_shift_warmup_cycles", 0))
),
smearing_temperature=(
None
if options is None
else float(getattr(options, "smearing_temperature", 0.0))
),
variational_space=variational_space,
finite_torus_convention=convention,
)
setattr(result, "backend", f"aiccm2026dev-b-{backend.value}")
setattr(result, "aiccm2026dev_b", diagnostics)
setattr(result, "finite_torus_convention", convention)
setattr(result, "ccm_approach", convention.ccm_approach)
setattr(result, "ccm_construction", convention.ccm_construction)
setattr(
result,
"evaluation_representation",
convention.evaluation_representation,
)
setattr(result, "coulomb_kernel", convention.coulomb_kernel)
setattr(result, "exchange_q0", convention.exchange_q0)
setattr(
result,
"exchange_q0_applicability",
convention.exchange_q0_applicability,
)
setattr(result, "boundary_model", convention.boundary_model)
setattr(result, "effective_n_electrons", effective_electrons)
setattr(result, "effective_nuclear_charge", charge_bookkeeping.effective_nuclear_charge)
setattr(result, "effective_nuclear_charges", charge_bookkeeping.effective_nuclear_charges)
setattr(result, "ecp_total_ncore", charge_bookkeeping.ecp_total_ncore)
setattr(result, "kpoints_frac", np.asarray(kpoints.kpoints_frac, dtype=float))
setattr(result, "kpoints_cart", np.asarray(kpoints.kpoints_cart, dtype=float))
setattr(result, "kpoint_weights", np.asarray(kpoints.weights, dtype=float))
_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()
requested_fock_mixing = resolve_fock_mixing(
opts,
fock_mixing,
where="run_aiccm2026dev_b_rhf",
)
_validate_closed_shell_problem(system, opts)
_require_supported_scf_dimension(system)
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"
)
if (
selected is AICCM2026DevBBackend.RI
and int(system.dim) == 3
and int(np.prod(mesh_tuple)) == 1
and requested_fock_mixing != 0.0
):
raise NotImplementedError(
"aiccm2026dev-b Gamma-only RI RHF cannot execute previous-Fock "
"mixing without switching to the legacy molecular-limit GDF "
"operator; use fock_mixing=0 or at least two cyclic cells"
)
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,
fock_mixing=requested_fock_mixing,
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=requested_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",
1.0,
perf_counter() - started,
symmetry_plan,
opts,
)
[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)
requested_fock_mixing = resolve_fock_mixing(
opts,
fock_mixing,
where="run_aiccm2026dev_b_rks",
)
_validate_closed_shell_problem(system, opts)
_require_supported_scf_dimension(system)
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)
c_full = _full_range_exchange_coefficient(str(functional), spin=1)
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),
fock_mixing=requested_fock_mixing,
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=requested_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}",
c_full,
perf_counter() - started,
symmetry_plan,
opts,
)
[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()
requested_fock_mixing = resolve_fock_mixing(
opts,
fock_mixing,
where="run_aiccm2026dev_b_uhf",
)
_validate_open_shell_problem(system, opts)
_require_supported_scf_dimension(system)
if (
selected is not AICCM2026DevBBackend.FOUR_CENTER
and requested_fock_mixing != 0.0
):
raise NotImplementedError(
f"aiccm2026dev-b {selected.value} UHF does not implement "
"previous-Fock mixing; use fock_mixing=0 or select "
"backend='four_center'"
)
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,
fock_mixing=requested_fock_mixing,
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",
1.0,
perf_counter() - started,
symmetry_plan,
opts,
)
[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)
requested_fock_mixing = resolve_fock_mixing(
opts,
fock_mixing,
where="run_aiccm2026dev_b_uks",
)
_validate_open_shell_problem(system, opts)
_require_supported_scf_dimension(system)
if (
selected is not AICCM2026DevBBackend.FOUR_CENTER
and requested_fock_mixing != 0.0
):
raise NotImplementedError(
f"aiccm2026dev-b {selected.value} UKS does not implement "
"previous-Fock mixing; use fock_mixing=0 or select "
"backend='four_center'"
)
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)
c_full = _full_range_exchange_coefficient(str(functional), spin=2)
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),
fock_mixing=requested_fock_mixing,
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}",
c_full,
perf_counter() - started,
symmetry_plan,
opts,
)