"""Momentum-conserving post-HF methods for AICCM2026DEV-B.
The routines in this module consume the same finite-group RI Hamiltonian as
``periodic_aiccm2026dev_b``. They do not use the historical CCM tensors or the
legacy Gamma-supercell HF reference.
"""
from __future__ import annotations
from contextlib import contextmanager
from dataclasses import dataclass, replace
from types import SimpleNamespace
from typing import Sequence
import numpy as np
from ._vibeqc_core import BasisSet, LatticeSumOptions, PeriodicRHFOptions, PeriodicSystem
from .aux_basis import (
build_lpq_bloch_native_fft,
default_aux_for,
make_aux_basis_set,
make_modrho_aux_basis,
)
from .periodic_aiccm2026dev_b import (
AICCM2026DevBDiagnostics,
_normalise_mesh,
_resolve_lattice_extension,
cyclic_gamma_mesh,
run_aiccm2026dev_b_rhf,
run_aiccm2026dev_b_uhf,
)
__all__ = [
"AICCM2026DevBLocalCorrelationSpace",
"AICCM2026DevBDLPNOCCSDResult",
"AICCM2026DevBDLPNOMP2Result",
"AICCM2026DevBMP2Result",
"AICCM2026DevBUCCSDResult",
"AICCM2026DevBUMP2Result",
"run_aiccm2026dev_b_dlpno_ccsd",
"run_aiccm2026dev_b_dlpno_ccsd_t",
"run_aiccm2026dev_b_dlpno_mp2",
"run_aiccm2026dev_b_dlpno_uccsd",
"run_aiccm2026dev_b_dlpno_uccsd_t",
"run_aiccm2026dev_b_dlpno_ump2",
"run_aiccm2026dev_b_ccsd",
"run_aiccm2026dev_b_ccsd_t",
"run_aiccm2026dev_b_mp2",
"run_aiccm2026dev_b_uccsd",
"run_aiccm2026dev_b_uccsd_t",
"run_aiccm2026dev_b_ump2",
]
[docs]
@dataclass(frozen=True)
class AICCM2026DevBLocalCorrelationSpace:
"""Validated finite-torus PAO and occupied-pair bookkeeping.
``representative_pair_acceleration`` stays false until the pair-amplitude
scatter has an energy-parity gate. The orbit count is therefore a
measured reduction opportunity, not a timing claim.
"""
pao_rank: int
canonical_virtual_rank: int
pao_discarded_rank: int
pao_orthonormality_error: float
pao_occupied_leakage_error: float
n_occupied: int
n_pairs: int
n_translation_unique_pairs: int
translation_reduction_factor: float
translation_covariance_error: float
translation_orbits_validated: bool
representative_pair_acceleration: bool = False
[docs]
@dataclass(frozen=True)
class AICCM2026DevBMP2Result:
"""Canonical finite-torus RI-MP2 energy per primitive cell."""
mesh: tuple[int, int, int]
n_cyclic_cells: int
n_kpoints: int
n_occ: int
n_vir: int
n_aux: int
e_hf_per_cell: float
e_corr_per_cell: float
e_corr_ss_per_cell: float
e_corr_os_per_cell: float
e_total_per_cell: float
max_energy_imaginary_residual: float
momentum_conservation_error: float
aux_basis_name: str
hf_diagnostics: AICCM2026DevBDiagnostics
hf_result: object
backend: str = "aiccm2026dev-b-ri-mp2"
@property
def energy(self) -> float:
"""Total RI-MP2 energy per primitive cell in hartree."""
return self.e_total_per_cell
@dataclass(frozen=True)
class AICCM2026DevBUMP2Result:
"""Unrestricted RI-MP2 or local-PNO UMP2 result per primitive cell."""
mesh: tuple[int, int, int]
n_cyclic_cells: int
e_hf_per_cell: float
e_corr_per_cell: float
e_aa_per_cell: float
e_bb_per_cell: float
e_ab_per_cell: float
e_total_per_cell: float
n_pairs_aa: int
n_pairs_bb: int
n_pairs_ab: int
converged: bool
localization: str
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
solver_result: object
hf_result: object
localization_result: object | None = None
backend: str = "aiccm2026dev-b-ri-ump2"
@property
def energy(self) -> float:
return self.e_total_per_cell
@dataclass(frozen=True)
class AICCM2026DevBUCCSDResult:
"""Full-domain or PNO-truncated UCCSD(T) pilot result per cell.
The current solver is the explicitly gated O(N^6) projection oracle. It
establishes the unrestricted finite-torus equations and exact PNO limit;
it is not presented as the later reduced-scaling production kernel.
"""
mesh: tuple[int, int, int]
n_cyclic_cells: int
e_hf_per_cell: float
e_corr_per_cell: float
e_t_per_cell: float
e_total_per_cell: float
n_pairs: int
n_iter: int
converged: bool
t1_norm: float
avg_pno: float
localization: str
with_triples: bool
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
solver_result: object
hf_result: object
localization_result: object | None = None
backend: str = "aiccm2026dev-b-uccsd"
@property
def energy(self) -> float:
return self.e_total_per_cell
[docs]
@dataclass(frozen=True)
class AICCM2026DevBDLPNOMP2Result:
"""Local-PNO MP2 result on the B finite torus, per primitive cell."""
mesh: tuple[int, int, int]
n_cyclic_cells: int
e_hf_per_cell: float
e_corr_per_cell: float
raw_local_e_corr_per_cell: float
complete_space_correction_per_cell: float
e_total_per_cell: float
n_pairs: int
n_pairs_screened: int
n_iter: int
converged: bool
localization: str
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
localization_result: object | None
local_correlation_space: AICCM2026DevBLocalCorrelationSpace | None
solver_result: object
hf_result: object
backend: str = "aiccm2026dev-b-dlpno-mp2"
@property
def energy(self) -> float:
"""Total local-PNO MP2 energy per primitive cell in hartree."""
return self.e_total_per_cell
[docs]
@dataclass(frozen=True)
class AICCM2026DevBDLPNOCCSDResult:
"""Local-PNO CCSD or CCSD(T) result on the B finite torus."""
mesh: tuple[int, int, int]
n_cyclic_cells: int
e_hf_per_cell: float
e_corr_per_cell: float
e_t_per_cell: float
e_total_per_cell: float
n_pairs: int
n_pairs_screened: int
n_iter: int
converged: bool
t1_norm: float
localization: str
with_triples: bool
triples_mode: str
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
localization_result: object | None
local_correlation_space: AICCM2026DevBLocalCorrelationSpace | None
solver_result: object
hf_result: object
backend: str
@property
def energy(self) -> float:
"""Total local-PNO CC energy per primitive cell in hartree."""
return self.e_total_per_cell
def _fractional_key(kfrac: np.ndarray, *, scale: int = 10**10) -> tuple[int, int, int]:
wrapped = np.mod(np.asarray(kfrac, dtype=float), 1.0)
wrapped[np.isclose(wrapped, 1.0, atol=1.0 / scale, rtol=0.0)] = 0.0
return tuple(int(value) % scale for value in np.rint(wrapped * scale))
def _momentum_conservation_table(kpoints_frac: np.ndarray) -> tuple[np.ndarray, float]:
"""Return ``kb[ki,ka,kj]`` for ``ki-ka+kj-kb = G``."""
kfrac = np.asarray(kpoints_frac, dtype=float).reshape(-1, 3)
lookup = {_fractional_key(value): index for index, value in enumerate(kfrac)}
if len(lookup) != len(kfrac):
raise ValueError("aiccm2026dev-b MP2 requires distinct finite-group characters")
n_k = len(kfrac)
table = np.empty((n_k, n_k, n_k), dtype=int)
maximum_error = 0.0
for ki in range(n_k):
for ka in range(n_k):
for kj in range(n_k):
target = kfrac[ki] - kfrac[ka] + kfrac[kj]
key = _fractional_key(target)
if key not in lookup:
raise RuntimeError(
"aiccm2026dev-b MP2 character mesh is not closed under "
"momentum conservation"
)
kb = lookup[key]
table[ki, ka, kj] = kb
residual = target - kfrac[kb]
residual -= np.rint(residual)
maximum_error = max(maximum_error, float(np.linalg.norm(residual)))
return table, maximum_error
def _build_canonical_lpq_cache(
system: PeriodicSystem,
basis: BasisSet,
kpoints_cart: np.ndarray,
aux_name: str,
*,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
) -> tuple[dict[tuple[int, int], np.ndarray], object]:
molecule = system.unit_cell_molecule()
auxiliary = make_aux_basis_set(molecule, aux_name=aux_name)
auxiliary_modrho = make_modrho_aux_basis(auxiliary, molecule)
lattice_options = LatticeSumOptions()
lattice_options.cutoff_bohr = float(lattice_cutoff_bohr)
kcart = np.asarray(kpoints_cart, dtype=float).reshape(-1, 3)
cache: dict[tuple[int, int], np.ndarray] = {}
for ki, bra_k in enumerate(kcart):
for ka, ket_k in enumerate(kcart):
cache[(ki, ka)] = build_lpq_bloch_native_fft(
system,
basis,
auxiliary_modrho,
bra_k,
ket_k,
ke_cutoff=float(rsgdf_ke_cutoff),
lat_opts=lattice_options,
linear_dep_thr=float(gdf_linear_dep_threshold),
canonical_auxiliary_basis=True,
)
return cache, auxiliary
def _cyclic_translations(mesh: tuple[int, int, int]) -> np.ndarray:
return np.asarray(
[
(i, j, k)
for i in range(mesh[0])
for j in range(mesh[1])
for k in range(mesh[2])
],
dtype=int,
)
def _periodic_matrix_to_real_supercell(
matrices_k: Sequence[np.ndarray],
phases: np.ndarray,
) -> tuple[np.ndarray, float]:
"""Transform character blocks to the AO matrix of the finite torus."""
matrices = np.asarray(matrices_k, dtype=complex)
n_cells = phases.shape[1]
blocks = np.einsum(
"kr,kmn,ks->rmsn",
phases,
matrices,
phases.conj(),
optimize=True,
) / n_cells
matrix = blocks.reshape(
n_cells * matrices.shape[1],
n_cells * matrices.shape[2],
)
imaginary = float(np.max(np.abs(matrix.imag)))
matrix = 0.5 * (matrix + matrix.conj().T)
return np.asarray(np.real_if_close(matrix, tol=1000).real), imaginary
def _lpq_cache_to_real_supercell(
lpq_cache: dict[tuple[int, int], np.ndarray],
kpoints_frac: np.ndarray,
mesh: tuple[int, int, int],
) -> tuple[np.ndarray, float, float]:
"""Inverse-transform pair-resolved RI factors to the real finite torus."""
translations = _cyclic_translations(mesh)
kfrac = np.asarray(kpoints_frac, dtype=float).reshape(-1, 3)
n_cells = len(translations)
if len(kfrac) != n_cells:
raise ValueError("AICCM2026DEV-B cderi transform needs one character per cell")
phases = np.exp(2j * np.pi * (kfrac @ translations.T))
sample = lpq_cache[(0, 0)]
n_aux, nbf, nbf_right = sample.shape
if nbf != nbf_right:
raise ValueError("AICCM2026DEV-B pair cderi matrices must be square")
factors = np.zeros(
(n_cells, n_aux, n_cells, nbf, n_cells, nbf),
dtype=complex,
)
scale = 1.0 / (n_cells * n_cells)
for t_index, translation in enumerate(translations):
for ki in range(n_cells):
for ka in range(n_cells):
auxiliary_phase = np.exp(
2j * np.pi * ((kfrac[ka] - kfrac[ki]) @ translation)
)
cell_phase = (
phases[ki, :, None]
* phases[ka, None, :].conj()
* auxiliary_phase
)
factors[t_index] += scale * np.einsum(
"rs,Pmn->Prmsn",
cell_phase,
lpq_cache[(ki, ka)],
optimize=True,
)
factors = factors.reshape(n_cells * n_aux, n_cells * nbf, n_cells * nbf)
imaginary = float(np.max(np.abs(factors.imag)))
symmetry = float(np.max(np.abs(factors - factors.swapaxes(1, 2))))
factors = 0.5 * (factors + factors.swapaxes(1, 2))
return np.asarray(np.real_if_close(factors, tol=1000).real), imaginary, symmetry
def _build_supercell_system(
system: PeriodicSystem,
mesh: tuple[int, int, int],
*,
multiplicity: int = 1,
) -> PeriodicSystem:
from ._vibeqc_core import Atom
lattice = np.asarray(system.lattice, dtype=float)
super_lattice = lattice.copy()
for axis, count in enumerate(mesh):
super_lattice[axis] *= count
atoms = []
for translation in _cyclic_translations(mesh):
shift = np.asarray(translation, dtype=float) @ lattice
for atom in system.unit_cell:
position = np.asarray(atom.xyz, dtype=float) + shift
atoms.append(Atom(int(atom.Z), position.tolist()))
return PeriodicSystem(
int(system.dim),
super_lattice,
atoms,
charge=int(system.charge) * int(np.prod(mesh)),
multiplicity=int(multiplicity),
)
class _BRealDFAdapter:
def __init__(self, factors: np.ndarray, auxiliary_basis: object):
self.B = factors
self.three_center = factors
self.metric = np.eye(factors.shape[0])
self.aux_basis = auxiliary_basis
self.n_aux = factors.shape[0]
self.n_orb = factors.shape[1]
def mo_transform(
self,
coefficients_left: np.ndarray,
coefficients_right: np.ndarray | None = None,
) -> np.ndarray:
right = coefficients_left if coefficients_right is None else coefficients_right
return np.einsum(
"mi,Pmn,nj->Pij",
np.asarray(coefficients_left),
self.B,
np.asarray(right),
optimize=True,
)
@contextmanager
def _patched_overlap(overlap: np.ndarray):
import vibeqc._vibeqc_core as core
original = core.compute_overlap
core.compute_overlap = lambda _basis: overlap
try:
yield
finally:
core.compute_overlap = original
@dataclass(frozen=True)
class _BRealReference:
super_system: PeriodicSystem
molecule: object
basis: BasisSet
df: _BRealDFAdapter
hf: object
n_cells: int
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
def _build_b_real_reference(
system: PeriodicSystem,
basis: BasisSet,
mesh: tuple[int, int, int],
hf_result: object,
lpq_cache: dict[tuple[int, int], np.ndarray],
auxiliary_name: str,
) -> _BRealReference:
kpoints = cyclic_gamma_mesh(system, mesh)
translations = _cyclic_translations(mesh)
phases = np.exp(
2j * np.pi * (np.asarray(kpoints.kpoints_frac) @ translations.T)
)
overlap, overlap_imaginary = _periodic_matrix_to_real_supercell(
hf_result.overlap,
phases,
)
fock, fock_imaginary = _periodic_matrix_to_real_supercell(
hf_result.fock,
phases,
)
factors, cderi_imaginary, cderi_symmetry = _lpq_cache_to_real_supercell(
lpq_cache,
np.asarray(kpoints.kpoints_frac),
mesh,
)
if max(cderi_imaginary, cderi_symmetry) > 1e-7:
raise RuntimeError(
"AICCM2026DEV-B real RI transform violated Hermitian symmetry: "
f"imaginary={cderi_imaginary:.3e}, "
f"symmetry={cderi_symmetry:.3e}"
)
overlap_eigenvalues, overlap_eigenvectors = np.linalg.eigh(overlap)
keep = overlap_eigenvalues > 1e-10 * float(overlap_eigenvalues[-1])
if int(np.count_nonzero(keep)) != overlap.shape[0]:
raise RuntimeError(
"AICCM2026DEV-B real-torus post-HF reference is linearly dependent"
)
orthogonalizer = overlap_eigenvectors[:, keep] / np.sqrt(
overlap_eigenvalues[keep]
)[None, :]
reduced_fock = orthogonalizer.T @ fock @ orthogonalizer
orbital_energies, reduced_coefficients = np.linalg.eigh(
0.5 * (reduced_fock + reduced_fock.T)
)
coefficients = orthogonalizer @ reduced_coefficients
super_system = _build_supercell_system(system, mesh)
molecule = super_system.unit_cell_molecule()
super_basis = BasisSet(molecule, basis.name)
if int(super_basis.nbasis) != factors.shape[1]:
raise RuntimeError(
"AICCM2026DEV-B supercell basis ordering does not match the "
"finite-group AO transform"
)
super_auxiliary = make_aux_basis_set(molecule, aux_name=auxiliary_name)
if int(super_auxiliary.nbasis) != factors.shape[0]:
raise RuntimeError(
"AICCM2026DEV-B supercell auxiliary ordering does not match the "
"finite-group cderi transform"
)
df = _BRealDFAdapter(factors, super_auxiliary)
n_cells = int(np.prod(mesh))
hf = SimpleNamespace(
energy=float(hf_result.energy) * n_cells,
converged=True,
mo_coeffs=coefficients,
mo_energies=orbital_energies,
fock=fock,
overlap=overlap,
)
return _BRealReference(
super_system=super_system,
molecule=molecule,
basis=super_basis,
df=df,
hf=hf,
n_cells=n_cells,
cderi_imaginary_residual=cderi_imaginary,
cderi_symmetry_residual=cderi_symmetry,
matrix_imaginary_residual=max(overlap_imaginary, fock_imaginary),
)
@dataclass(frozen=True)
class _BRealUReference:
super_system: PeriodicSystem
molecule: object
basis: BasisSet
df: _BRealDFAdapter
hf: object
n_cells: int
cderi_imaginary_residual: float
cderi_symmetry_residual: float
matrix_imaginary_residual: float
def _metric_orbitals(overlap: np.ndarray, fock: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
overlap_eigenvalues, overlap_eigenvectors = np.linalg.eigh(overlap)
keep = overlap_eigenvalues > 1e-10 * float(overlap_eigenvalues[-1])
if int(np.count_nonzero(keep)) != overlap.shape[0]:
raise RuntimeError(
"AICCM2026DEV-B real-torus post-HF reference is linearly dependent"
)
orthogonalizer = overlap_eigenvectors[:, keep] / np.sqrt(
overlap_eigenvalues[keep]
)[None, :]
reduced_fock = orthogonalizer.T @ fock @ orthogonalizer
energies, reduced_coefficients = np.linalg.eigh(
0.5 * (reduced_fock + reduced_fock.T)
)
return energies, orthogonalizer @ reduced_coefficients
def _build_b_real_ureference(
system: PeriodicSystem,
basis: BasisSet,
mesh: tuple[int, int, int],
hf_result: object,
lpq_cache: dict[tuple[int, int], np.ndarray],
auxiliary_name: str,
) -> _BRealUReference:
"""Transform the finite-character UHF Hamiltonian to one real torus."""
kpoints = cyclic_gamma_mesh(system, mesh)
translations = _cyclic_translations(mesh)
phases = np.exp(
2j * np.pi * (np.asarray(kpoints.kpoints_frac) @ translations.T)
)
overlap, overlap_imaginary = _periodic_matrix_to_real_supercell(
hf_result.overlap,
phases,
)
fock_alpha, fock_alpha_imaginary = _periodic_matrix_to_real_supercell(
hf_result.fock_alpha,
phases,
)
fock_beta, fock_beta_imaginary = _periodic_matrix_to_real_supercell(
hf_result.fock_beta,
phases,
)
factors, cderi_imaginary, cderi_symmetry = _lpq_cache_to_real_supercell(
lpq_cache,
np.asarray(kpoints.kpoints_frac),
mesh,
)
if max(cderi_imaginary, cderi_symmetry) > 1e-7:
raise RuntimeError(
"AICCM2026DEV-B real RI transform violated Hermitian symmetry: "
f"imaginary={cderi_imaginary:.3e}, "
f"symmetry={cderi_symmetry:.3e}"
)
energies_alpha, coefficients_alpha = _metric_orbitals(overlap, fock_alpha)
energies_beta, coefficients_beta = _metric_orbitals(overlap, fock_beta)
n_cells = int(np.prod(mesh))
super_multiplicity = n_cells * (int(system.multiplicity) - 1) + 1
super_system = _build_supercell_system(
system,
mesh,
multiplicity=super_multiplicity,
)
molecule = super_system.unit_cell_molecule()
super_basis = BasisSet(molecule, basis.name)
if int(super_basis.nbasis) != factors.shape[1]:
raise RuntimeError(
"AICCM2026DEV-B supercell basis ordering does not match the "
"finite-group AO transform"
)
super_auxiliary = make_aux_basis_set(molecule, aux_name=auxiliary_name)
if int(super_auxiliary.nbasis) != factors.shape[0]:
raise RuntimeError(
"AICCM2026DEV-B supercell auxiliary ordering does not match the "
"finite-group cderi transform"
)
df = _BRealDFAdapter(factors, super_auxiliary)
hf = SimpleNamespace(
energy=float(hf_result.energy) * n_cells,
converged=True,
mo_coeffs_alpha=coefficients_alpha,
mo_coeffs_beta=coefficients_beta,
mo_energies_alpha=energies_alpha,
mo_energies_beta=energies_beta,
fock_alpha=fock_alpha,
fock_beta=fock_beta,
overlap=overlap,
)
return _BRealUReference(
super_system=super_system,
molecule=molecule,
basis=super_basis,
df=df,
hf=hf,
n_cells=n_cells,
cderi_imaginary_residual=cderi_imaginary,
cderi_symmetry_residual=cderi_symmetry,
matrix_imaginary_residual=max(
overlap_imaginary,
fock_alpha_imaginary,
fock_beta_imaginary,
),
)
def _prepare_b_real_reference(
system: PeriodicSystem,
basis: BasisSet,
mesh: tuple[int, int, int],
options: PeriodicRHFOptions | None,
*,
auxiliary_name: str,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> tuple[_BRealReference, object]:
if int(system.dim) != 3:
raise NotImplementedError(
"aiccm2026dev-b local correlation is currently restricted to 3D "
"because the 1D/2D fitted and direct long-range gauges are not "
"yet matched"
)
hf_result = run_aiccm2026dev_b_rhf(
system,
basis,
mesh,
options,
backend="ri",
aux_basis=auxiliary_name,
gdf_method="rsgdf",
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
progress=progress,
verbose=verbose,
)
if not bool(hf_result.converged):
raise RuntimeError(
"aiccm2026dev-b local correlation requires a converged RI-RHF reference"
)
kpoints = cyclic_gamma_mesh(system, mesh)
lpq_cache, _ = _build_canonical_lpq_cache(
system,
basis,
np.asarray(kpoints.kpoints_cart, dtype=float),
auxiliary_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
)
return (
_build_b_real_reference(
system,
basis,
mesh,
hf_result,
lpq_cache,
auxiliary_name,
),
hf_result,
)
def _prepare_b_real_ureference(
system: PeriodicSystem,
basis: BasisSet,
mesh: tuple[int, int, int],
options: PeriodicRHFOptions | None,
*,
auxiliary_name: str,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> tuple[_BRealUReference, object]:
if int(system.dim) != 3:
raise NotImplementedError(
"aiccm2026dev-b unrestricted local correlation is currently "
"restricted to 3D because the 1D/2D fitted and direct long-range "
"gauges are not yet matched"
)
hf_result = run_aiccm2026dev_b_uhf(
system,
basis,
mesh,
options,
backend="ri",
aux_basis=auxiliary_name,
gdf_method="rsgdf",
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
progress=progress,
verbose=verbose,
)
if not bool(hf_result.converged):
raise RuntimeError(
"aiccm2026dev-b unrestricted correlation requires a converged RI-UHF reference"
)
kpoints = cyclic_gamma_mesh(system, mesh)
lpq_cache, _ = _build_canonical_lpq_cache(
system,
basis,
np.asarray(kpoints.kpoints_cart, dtype=float),
auxiliary_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
)
return (
_build_b_real_ureference(
system,
basis,
mesh,
hf_result,
lpq_cache,
auxiliary_name,
),
hf_result,
)
def _localize_b_real_reference(
reference: _BRealReference,
method: str,
*,
hf_result: object,
system: PeriodicSystem,
primitive_basis: BasisSet,
mesh: tuple[int, int, int],
) -> tuple[object, str, object | None]:
if method == "none":
return reference.hf, method, None
if method in {"wannier", "iao"}:
from .periodic_aiccm2026dev_b_localization import (
localize_aiccm2026dev_b_occupied,
)
localization = localize_aiccm2026dev_b_occupied(
hf_result,
system,
primitive_basis,
mesh,
method=method,
)
n_occ = int(reference.molecule.n_electrons()) // 2
occupied = np.asarray(localization.coefficients)
if occupied.shape != (int(reference.basis.nbasis), n_occ):
raise RuntimeError(
"AICCM2026DEV-B localized occupied coefficients do not match "
"the real-torus reference"
)
if np.iscomplexobj(occupied):
real_trial = occupied.real
gram = real_trial.T @ reference.hf.overlap @ real_trial
eigenvalues, eigenvectors = np.linalg.eigh(
0.5 * (gram + gram.T)
)
if eigenvalues[0] <= 1e-10 * eigenvalues[-1]:
raise RuntimeError(
"AICCM2026DEV-B localization has no full-rank real "
"time-reversal gauge for the real local solver"
)
occupied_real = real_trial @ (
eigenvectors / np.sqrt(eigenvalues)[None, :]
) @ eigenvectors.T
projector_error = float(
np.linalg.norm(
occupied @ occupied.conj().T
- occupied_real @ occupied_real.T
)
)
if projector_error > 1e-7:
raise RuntimeError(
"AICCM2026DEV-B real-gauge projection changed the "
f"occupied projector; error={projector_error:.3e}"
)
occupied = occupied_real
coefficients = np.asarray(reference.hf.mo_coeffs)
localized_coefficients = np.concatenate(
[occupied, coefficients[:, n_occ:]],
axis=1,
)
metric_error = float(
np.max(
np.abs(
localized_coefficients.conj().T
@ reference.hf.overlap
@ localized_coefficients
- np.eye(localized_coefficients.shape[1])
)
)
)
if metric_error > 1e-7:
raise RuntimeError(
"AICCM2026DEV-B localized occupied/virtual union is not "
f"metric orthonormal; error={metric_error:.3e}"
)
proxy = SimpleNamespace(
energy=reference.hf.energy,
converged=True,
mo_coeffs=np.asarray(localized_coefficients, dtype=float),
mo_energies=reference.hf.mo_energies,
fock=reference.hf.fock,
overlap=reference.hf.overlap,
)
return proxy, method, localization
if method != "pipek-mezey":
raise ValueError(
"aiccm2026dev-b local correlation supports localise='wannier', "
"'iao', 'pipek-mezey', or 'none'; molecular Boys localization "
"is not periodic-boundary safe"
)
from .periodic_localise import localise_periodic_gamma
n_occ = int(reference.molecule.n_electrons()) // 2
localized = localise_periodic_gamma(
reference.hf,
reference.basis,
reference.super_system,
method="pipek-mezey",
n_occ=n_occ,
)
coefficients = np.asarray(reference.hf.mo_coeffs)
localized_coefficients = np.concatenate(
[localized.C_loc, coefficients[:, n_occ:]],
axis=1,
)
proxy = SimpleNamespace(
energy=reference.hf.energy,
converged=True,
mo_coeffs=localized_coefficients,
mo_energies=reference.hf.mo_energies,
fock=reference.hf.fock,
overlap=reference.hf.overlap,
)
return proxy, method, None
def _real_occupied_gauge(
occupied: np.ndarray,
overlap: np.ndarray,
*,
label: str,
) -> np.ndarray:
if not np.iscomplexobj(occupied):
return np.asarray(occupied, dtype=float)
real_trial = occupied.real
gram = real_trial.T @ overlap @ real_trial
eigenvalues, eigenvectors = np.linalg.eigh(0.5 * (gram + gram.T))
if eigenvalues[0] <= 1e-10 * eigenvalues[-1]:
raise RuntimeError(
f"AICCM2026DEV-B {label} localization has no full-rank real "
"time-reversal gauge"
)
real_occupied = real_trial @ (
eigenvectors / np.sqrt(eigenvalues)[None, :]
) @ eigenvectors.T
projector_error = float(
np.linalg.norm(
occupied @ occupied.conj().T
- real_occupied @ real_occupied.T
)
)
if projector_error > 1e-7:
raise RuntimeError(
f"AICCM2026DEV-B {label} real-gauge projection changed the "
f"occupied projector; error={projector_error:.3e}"
)
return real_occupied
def _localize_b_real_ureference(
reference: _BRealUReference,
method: str,
*,
hf_result: object,
system: PeriodicSystem,
primitive_basis: BasisSet,
mesh: tuple[int, int, int],
) -> tuple[object, str, object | None]:
if method == "none":
return reference.hf, method, None
if method not in {"wannier", "iao"}:
raise ValueError(
"aiccm2026dev-b unrestricted local correlation supports the "
"PBC-safe localise='wannier' or 'iao' gauges, or exact-limit "
"localise='none'"
)
from .periodic_aiccm2026dev_b_localization import (
localize_aiccm2026dev_b_unrestricted_occupied,
)
localization = localize_aiccm2026dev_b_unrestricted_occupied(
hf_result,
system,
primitive_basis,
mesh,
method=method,
)
ne = int(reference.molecule.n_electrons())
two_s = int(reference.molecule.multiplicity) - 1
n_alpha = (ne + two_s) // 2
n_beta = (ne - two_s) // 2
occupied_alpha = _real_occupied_gauge(
np.asarray(localization.alpha.coefficients),
np.asarray(reference.hf.overlap),
label="alpha",
)
if occupied_alpha.shape[1] != n_alpha:
raise RuntimeError("AICCM2026DEV-B alpha localized rank is inconsistent")
if n_beta:
if localization.beta is None:
raise RuntimeError("AICCM2026DEV-B beta localization is missing")
occupied_beta = _real_occupied_gauge(
np.asarray(localization.beta.coefficients),
np.asarray(reference.hf.overlap),
label="beta",
)
else:
occupied_beta = np.empty((int(reference.basis.nbasis), 0))
if occupied_beta.shape[1] != n_beta:
raise RuntimeError("AICCM2026DEV-B beta localized rank is inconsistent")
coefficients_alpha = np.concatenate(
[occupied_alpha, np.asarray(reference.hf.mo_coeffs_alpha)[:, n_alpha:]],
axis=1,
)
coefficients_beta = np.concatenate(
[occupied_beta, np.asarray(reference.hf.mo_coeffs_beta)[:, n_beta:]],
axis=1,
)
for label, coefficients in (
("alpha", coefficients_alpha),
("beta", coefficients_beta),
):
error = float(
np.max(
np.abs(
coefficients.T @ reference.hf.overlap @ coefficients
- np.eye(coefficients.shape[1])
)
)
)
if error > 1e-7:
raise RuntimeError(
f"AICCM2026DEV-B {label} localized occupied/virtual union "
f"is not metric orthonormal; error={error:.3e}"
)
proxy = SimpleNamespace(
energy=reference.hf.energy,
converged=True,
mo_coeffs_alpha=coefficients_alpha,
mo_coeffs_beta=coefficients_beta,
mo_energies_alpha=reference.hf.mo_energies_alpha,
mo_energies_beta=reference.hf.mo_energies_beta,
fock_alpha=reference.hf.fock_alpha,
fock_beta=reference.hf.fock_beta,
overlap=reference.hf.overlap,
)
return proxy, method, localization
def _local_correlation_space(
reference: _BRealReference,
localized_hf: object,
localization_result: object | None,
) -> AICCM2026DevBLocalCorrelationSpace | None:
"""Build exact-limit PAOs and finite-translation pair orbits.
A translation permutation is only available from the B-specific
Wannier/IAO path. Pipek-Mezey and canonical real-torus orbitals may mix
degenerate translation sectors, so manufacturing a cell label for them
would be unsafe.
"""
if localization_result is None:
return None
from .periodic_aiccm2026dev_b_pno import (
enumerate_pair_orbits,
projected_atomic_orbitals,
)
n_occ = int(reference.molecule.n_electrons()) // 2
occupied = np.asarray(localized_hf.mo_coeffs)[:, :n_occ]
paos = projected_atomic_orbitals(
occupied,
np.asarray(reference.hf.overlap),
np.asarray(reference.hf.fock),
)
permutations = np.asarray(localization_result.translation_permutations)
orbits = enumerate_pair_orbits(permutations, n_cells=reference.n_cells)
n_pairs = n_occ * (n_occ + 1) // 2
if sum(orbit.multiplicity for orbit in orbits) != n_pairs:
raise RuntimeError("AICCM2026DEV-B translation pair orbits are incomplete")
return AICCM2026DevBLocalCorrelationSpace(
pao_rank=paos.rank,
canonical_virtual_rank=int(reference.basis.nbasis) - n_occ,
pao_discarded_rank=paos.discarded_rank,
pao_orthonormality_error=paos.orthonormality_error,
pao_occupied_leakage_error=paos.occupied_leakage_error,
n_occupied=n_occ,
n_pairs=n_pairs,
n_translation_unique_pairs=len(orbits),
translation_reduction_factor=n_pairs / len(orbits),
translation_covariance_error=float(
localization_result.translation_covariance_error
),
translation_orbits_validated=bool(
localization_result.translation_covariance_error < 1e-8
),
)
def _complete_space_mp2_audit(
reference: _BRealReference,
) -> float:
"""Evaluate the real canonical gauge of the complete-PAO MP2 limit."""
from .periodic_aiccm2026dev_b_pno import (
full_space_noncanonical_mp2_energy,
projected_atomic_orbitals,
)
n_occ = int(reference.molecule.n_electrons()) // 2
occupied = np.asarray(reference.hf.mo_coeffs)[:, :n_occ]
paos = projected_atomic_orbitals(
occupied,
np.asarray(reference.hf.overlap),
np.asarray(reference.hf.fock),
)
return full_space_noncanonical_mp2_energy(
occupied,
paos.coefficients,
np.asarray(reference.hf.fock),
np.asarray(reference.df.B),
)
[docs]
def run_aiccm2026dev_b_mp2(
system: PeriodicSystem,
basis: BasisSet,
mesh: int | Sequence[int] = (1, 1, 1),
options: PeriodicRHFOptions | None = None,
*,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBMP2Result:
"""Run canonical closed-shell RI-MP2 on the B finite translation group.
The RHF reference is the B-stream pair-resolved RI calculation on the full
unreduced character mesh. Three crystal momenta are summed independently;
the fourth is fixed by ``ki - ka + kj - kb = G``. The returned energy is
per primitive cell.
"""
mesh_tuple = _normalise_mesh(system, mesh)
if int(system.dim) != 3:
raise NotImplementedError(
"aiccm2026dev-b MP2 is currently restricted to 3D because the "
"1D/2D fitted and direct long-range gauges are not yet matched"
)
if lattice_cutoff_bohr <= 0.0:
raise ValueError("aiccm2026dev-b MP2 lattice cutoff must be positive")
if rsgdf_ke_cutoff <= 0.0:
raise ValueError("aiccm2026dev-b MP2 RSGDF cutoff must be positive")
if gdf_linear_dep_threshold <= 0.0:
raise ValueError("aiccm2026dev-b MP2 linear-dependence threshold must be positive")
aux_name = aux_basis or default_aux_for(basis.name)
hf_result = run_aiccm2026dev_b_rhf(
system,
basis,
mesh_tuple,
options,
backend="ri",
aux_basis=aux_name,
gdf_method="rsgdf",
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
progress=progress,
verbose=verbose,
)
if not bool(hf_result.converged):
raise RuntimeError("aiccm2026dev-b MP2 requires a converged RI-RHF reference")
kpoints = cyclic_gamma_mesh(system, mesh_tuple)
kfrac = np.asarray(kpoints.kpoints_frac, dtype=float)
kcart = np.asarray(kpoints.kpoints_cart, dtype=float)
n_k = len(kcart)
n_occ = int(system.n_electrons()) // 2
n_orb = int(basis.nbasis)
n_vir = n_orb - n_occ
if n_vir < 1:
raise RuntimeError("aiccm2026dev-b MP2 requires at least one virtual orbital")
coefficients = [np.asarray(value) for value in hf_result.mo_coeffs]
energies = [np.asarray(value, dtype=float) for value in hf_result.mo_energies]
lpq_cache, auxiliary = _build_canonical_lpq_cache(
system,
basis,
kcart,
aux_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
)
lov: dict[tuple[int, int], np.ndarray] = {}
n_aux = int(auxiliary.nbasis)
for ki in range(n_k):
occupied = coefficients[ki][:, :n_occ]
for ka in range(n_k):
virtual = coefficients[ka][:, n_occ:]
lpq = lpq_cache[(ki, ka)]
if lpq.shape[0] != n_aux:
raise RuntimeError("canonical AICCM2026DEV-B MP2 cderi rank mismatch")
lov[(ki, ka)] = np.einsum(
"mi,Pmn,na->Pia",
occupied.conj(),
lpq,
virtual,
optimize=True,
)
kconserv, momentum_error = _momentum_conservation_table(kfrac)
e_ss = 0.0
e_os = 0.0
maximum_imaginary = 0.0
oovv = [
np.empty((n_occ, n_occ, n_vir, n_vir), dtype=complex)
for _ in range(n_k)
]
for ki in range(n_k):
eps_i = energies[ki][:n_occ]
for kj in range(n_k):
eps_j = energies[kj][:n_occ]
for ka in range(n_k):
kb = int(kconserv[ki, ka, kj])
oovv[ka] = (1.0 / n_k) * np.einsum(
"Pia,Pjb->ijab",
lov[(ki, ka)],
lov[(kj, kb)],
optimize=True,
)
for ka in range(n_k):
kb = int(kconserv[ki, ka, kj])
denominator = (
eps_i[:, None, None, None]
+ eps_j[None, :, None, None]
- energies[ka][None, None, n_occ:, None]
- energies[kb][None, None, None, n_occ:]
)
if np.any(denominator >= -1e-12):
raise RuntimeError(
"aiccm2026dev-b MP2 encountered a non-negative orbital "
"energy denominator"
)
amplitudes = np.conj(oovv[ka] / denominator)
direct_complex = 2.0 * np.einsum(
"ijab,ijab->", amplitudes, oovv[ka], optimize=True
)
exchange_complex = -np.einsum(
"ijab,ijba->", amplitudes, oovv[kb], optimize=True
)
maximum_imaginary = max(
maximum_imaginary,
abs(float(np.imag(direct_complex))),
abs(float(np.imag(exchange_complex))),
)
direct = float(np.real(direct_complex))
exchange = float(np.real(exchange_complex))
e_ss += 0.5 * direct + exchange
e_os += 0.5 * direct
e_ss /= n_k
e_os /= n_k
e_corr = e_ss + e_os
e_hf = float(hf_result.energy)
return AICCM2026DevBMP2Result(
mesh=mesh_tuple,
n_cyclic_cells=int(np.prod(mesh_tuple)),
n_kpoints=n_k,
n_occ=n_occ,
n_vir=n_vir,
n_aux=n_aux,
e_hf_per_cell=e_hf,
e_corr_per_cell=e_corr,
e_corr_ss_per_cell=e_ss,
e_corr_os_per_cell=e_os,
e_total_per_cell=e_hf + e_corr,
max_energy_imaginary_residual=maximum_imaginary,
momentum_conservation_error=momentum_error,
aux_basis_name=aux_name,
hf_diagnostics=hf_result.aiccm2026dev_b,
hf_result=hf_result,
)
[docs]
def run_aiccm2026dev_b_dlpno_mp2(
system: PeriodicSystem,
basis: BasisSet,
mesh: int | Sequence[int] = (1, 1, 1),
options: PeriodicRHFOptions | None = None,
*,
dlpno_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBDLPNOMP2Result:
"""Run local-PNO MP2 on the exact real representation of the B torus.
The default uses PBC-safe Pipek-Mezey occupied localization. Pair-distance
screening and local fitting remain disabled because their current
molecular distance definitions are not invariant across a torus boundary.
"""
from .dlpno.mp2 import DLPNOMP2Options, run_dlpno_mp2
mesh_tuple = _normalise_mesh(system, mesh)
aux_name = aux_basis or default_aux_for(basis.name)
local_options = (
DLPNOMP2Options(
localise="pipek-mezey",
tcut_pairs=0.0,
tcut_pairs_weak=0.0,
)
if dlpno_options is None
else dlpno_options
)
if not isinstance(local_options, DLPNOMP2Options):
raise TypeError("dlpno_options must be a DLPNOMP2Options instance")
if local_options.localise not in {"wannier", "iao", "pipek-mezey", "none"}:
raise ValueError(
"aiccm2026dev-b DLPNO-MP2 requires PBC-safe "
"localise='wannier', 'iao', 'pipek-mezey', or exact-limit "
"localise='none'"
)
if local_options.tcut_pairs != 0.0 or local_options.tcut_pairs_weak != 0.0:
raise NotImplementedError(
"aiccm2026dev-b DLPNO-MP2 pair-distance screening awaits a "
"minimum-image periodic distance and transition-moment estimator"
)
if local_options.local_df:
raise NotImplementedError(
"aiccm2026dev-b DLPNO-MP2 local fitting awaits minimum-image "
"auxiliary domains"
)
if min(lattice_cutoff_bohr, rsgdf_ke_cutoff, gdf_linear_dep_threshold) <= 0.0:
raise ValueError("aiccm2026dev-b local-correlation cutoffs must be positive")
reference, hf_result = _prepare_b_real_reference(
system,
basis,
mesh_tuple,
options,
auxiliary_name=aux_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
localized_hf, localization, localization_result = _localize_b_real_reference(
reference,
local_options.localise,
hf_result=hf_result,
system=system,
primitive_basis=basis,
mesh=mesh_tuple,
)
local_space = _local_correlation_space(
reference,
localized_hf,
localization_result,
)
solver_options = replace(local_options, localise="none")
with _patched_overlap(reference.hf.overlap):
result = run_dlpno_mp2(
reference.molecule,
reference.basis,
localized_hf,
reference.df,
solver_options,
)
raw_local_correlation = float(result.e_corr)
complete_space_correction = 0.0
complete_space_limit = (
local_options.tcut_pno == 0.0
and local_options.tcut_pno_weak == 0.0
and local_options.tcut_mkn == 0.0
)
if complete_space_limit:
exact_correlation = _complete_space_mp2_audit(reference)
complete_space_correction = exact_correlation - raw_local_correlation
result.e_pno_correction += complete_space_correction
result.e_corr = exact_correlation
result.e_total = float(result.e_hf) + exact_correlation
n_cells = reference.n_cells
return AICCM2026DevBDLPNOMP2Result(
mesh=mesh_tuple,
n_cyclic_cells=n_cells,
e_hf_per_cell=float(result.e_hf) / n_cells,
e_corr_per_cell=float(result.e_corr) / n_cells,
raw_local_e_corr_per_cell=raw_local_correlation / n_cells,
complete_space_correction_per_cell=complete_space_correction / n_cells,
e_total_per_cell=float(result.e_total) / n_cells,
n_pairs=int(result.n_pairs),
n_pairs_screened=int(result.n_pairs_screened),
n_iter=int(result.n_iter),
converged=bool(result.converged),
localization=localization,
cderi_imaginary_residual=reference.cderi_imaginary_residual,
cderi_symmetry_residual=reference.cderi_symmetry_residual,
matrix_imaginary_residual=reference.matrix_imaginary_residual,
localization_result=localization_result,
local_correlation_space=local_space,
solver_result=result,
hf_result=hf_result,
)
def _run_aiccm2026dev_b_dlpno_cc(
system: PeriodicSystem,
basis: BasisSet,
mesh: int | Sequence[int],
options: PeriodicRHFOptions | None,
*,
cc_options: object | None,
with_triples: bool,
aux_basis: str | None,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> AICCM2026DevBDLPNOCCSDResult:
from .dlpno.ccsd_local_solver import LocalCCSDOptions, run_local_dlpno_ccsd
mesh_tuple = _normalise_mesh(system, mesh)
aux_name = aux_basis or default_aux_for(basis.name)
local_options = (
LocalCCSDOptions(
localise="pipek-mezey",
tcut_pairs=0.0,
coupling_radius=0.0,
compute_triples=with_triples,
)
if cc_options is None
else cc_options
)
if not isinstance(local_options, LocalCCSDOptions):
raise TypeError("cc_options must be a LocalCCSDOptions instance")
if local_options.localise not in {"wannier", "iao", "pipek-mezey", "none"}:
raise ValueError(
"aiccm2026dev-b DLPNO-CC requires PBC-safe "
"localise='wannier', 'iao', 'pipek-mezey', or exact-limit "
"localise='none'"
)
if local_options.coupling_radius != 0.0:
raise NotImplementedError(
"aiccm2026dev-b DLPNO-CC occupied-distance screening awaits "
"minimum-image periodic Wannier distances"
)
if min(lattice_cutoff_bohr, rsgdf_ke_cutoff, gdf_linear_dep_threshold) <= 0.0:
raise ValueError("aiccm2026dev-b local-correlation cutoffs must be positive")
reference, hf_result = _prepare_b_real_reference(
system,
basis,
mesh_tuple,
options,
auxiliary_name=aux_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
localized_hf, localization, localization_result = _localize_b_real_reference(
reference,
local_options.localise,
hf_result=hf_result,
system=system,
primitive_basis=basis,
mesh=mesh_tuple,
)
local_space = _local_correlation_space(
reference,
localized_hf,
localization_result,
)
solver_options = replace(
local_options,
localise="none",
compute_triples=with_triples,
)
with _patched_overlap(reference.hf.overlap):
result = run_local_dlpno_ccsd(
reference.molecule,
reference.basis,
localized_hf,
reference.df,
solver_options,
)
n_cells = reference.n_cells
triples = float(result.e_t) if with_triples else 0.0
total = float(result.e_hf) + float(result.e_corr) + triples
method = "dlpno-ccsd(t)" if with_triples else "dlpno-ccsd"
return AICCM2026DevBDLPNOCCSDResult(
mesh=mesh_tuple,
n_cyclic_cells=n_cells,
e_hf_per_cell=float(result.e_hf) / n_cells,
e_corr_per_cell=float(result.e_corr) / n_cells,
e_t_per_cell=triples / n_cells,
e_total_per_cell=total / n_cells,
n_pairs=int(result.n_pairs),
n_pairs_screened=int(result.n_screened),
n_iter=int(result.n_iter),
converged=bool(result.converged),
t1_norm=float(result.t1_norm),
localization=localization,
with_triples=with_triples,
triples_mode=str(solver_options.triples_mode),
cderi_imaginary_residual=reference.cderi_imaginary_residual,
cderi_symmetry_residual=reference.cderi_symmetry_residual,
matrix_imaginary_residual=reference.matrix_imaginary_residual,
localization_result=localization_result,
local_correlation_space=local_space,
solver_result=result,
hf_result=hf_result,
backend=f"aiccm2026dev-b-{method}",
)
[docs]
def run_aiccm2026dev_b_dlpno_ccsd(
system: PeriodicSystem,
basis: BasisSet,
mesh: int | Sequence[int] = (1, 1, 1),
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBDLPNOCCSDResult:
"""Run local-PNO CCSD on the real finite-torus B Hamiltonian."""
return _run_aiccm2026dev_b_dlpno_cc(
system,
basis,
mesh,
options,
cc_options=cc_options,
with_triples=False,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_dlpno_ccsd_t(
system: PeriodicSystem,
basis: BasisSet,
mesh: int | Sequence[int] = (1, 1, 1),
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBDLPNOCCSDResult:
"""Run local-PNO CCSD(T) on the real finite-torus B Hamiltonian."""
return _run_aiccm2026dev_b_dlpno_cc(
system,
basis,
mesh,
options,
cc_options=cc_options,
with_triples=True,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
def _run_aiccm2026dev_b_canonical_cc(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None,
options: PeriodicRHFOptions | None,
*,
with_triples: bool,
cc_options: object | None,
aux_basis: str | None,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> AICCM2026DevBDLPNOCCSDResult:
"""Select the complete-domain limit of the restricted local solver."""
from .dlpno.ccsd_local_solver import LocalCCSDOptions
selected = (
LocalCCSDOptions(
localise="none",
tcut_pno=0.0,
tcut_mkn=0.0,
tcut_pairs=0.0,
coupling_radius=0.0,
compute_triples=with_triples,
triples_mode="exact",
)
if cc_options is None
else cc_options
)
if not isinstance(selected, LocalCCSDOptions):
raise TypeError("cc_options must be a LocalCCSDOptions instance")
if (
selected.localise != "none"
or selected.tcut_pno != 0.0
or selected.tcut_mkn != 0.0
or selected.tcut_pairs != 0.0
or selected.coupling_radius != 0.0
):
raise ValueError(
"canonical aiccm2026dev-b CC requires canonical occupieds, "
"complete PNO domains, all pairs, and complete occupied coupling"
)
result = _run_aiccm2026dev_b_dlpno_cc(
system,
basis,
_resolve_lattice_extension(system, lattice_extension),
options,
cc_options=replace(selected, compute_triples=with_triples),
with_triples=with_triples,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
method = "ccsd(t)" if with_triples else "ccsd"
return replace(result, backend=f"aiccm2026dev-b-{method}")
[docs]
def run_aiccm2026dev_b_ccsd(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBDLPNOCCSDResult:
"""Run complete-domain finite-torus restricted CCSD."""
return _run_aiccm2026dev_b_canonical_cc(
system,
basis,
lattice_extension,
options,
with_triples=False,
cc_options=cc_options,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_ccsd_t(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBDLPNOCCSDResult:
"""Run complete-domain finite-torus restricted CCSD(T)."""
return _run_aiccm2026dev_b_canonical_cc(
system,
basis,
lattice_extension,
options,
with_triples=True,
cc_options=cc_options,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
def _run_aiccm2026dev_b_ump2(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None,
options: PeriodicRHFOptions | None,
*,
ump2_options: object,
aux_basis: str | None,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> AICCM2026DevBUMP2Result:
from .dlpno.ump2 import DLPNOUMP2Options, run_dlpno_ump2
if not isinstance(ump2_options, DLPNOUMP2Options):
raise TypeError("ump2_options must be a DLPNOUMP2Options instance")
if ump2_options.localise not in {"none", "wannier", "iao"}:
raise ValueError(
"aiccm2026dev-b UMP2 requires PBC-safe localise='wannier' or "
"'iao', or exact-limit localise='none'"
)
extension = _resolve_lattice_extension(system, lattice_extension)
aux_name = aux_basis or default_aux_for(basis.name)
reference, hf_result = _prepare_b_real_ureference(
system,
basis,
extension,
options,
auxiliary_name=aux_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
localized_hf, localization, localization_result = _localize_b_real_ureference(
reference,
ump2_options.localise,
hf_result=hf_result,
system=system,
primitive_basis=basis,
mesh=extension,
)
solver_options = replace(
ump2_options,
localise=("none" if localization == "none" else "external"),
)
with _patched_overlap(reference.hf.overlap):
result = run_dlpno_ump2(
reference.molecule,
reference.basis,
localized_hf,
reference.df,
solver_options,
)
n_cells = reference.n_cells
local = localization != "none"
return AICCM2026DevBUMP2Result(
mesh=extension,
n_cyclic_cells=n_cells,
e_hf_per_cell=float(result.e_hf) / n_cells,
e_corr_per_cell=float(result.e_corr) / n_cells,
e_aa_per_cell=float(result.e_aa) / n_cells,
e_bb_per_cell=float(result.e_bb) / n_cells,
e_ab_per_cell=float(result.e_ab) / n_cells,
e_total_per_cell=float(result.e_total) / n_cells,
n_pairs_aa=int(result.n_pairs_aa),
n_pairs_bb=int(result.n_pairs_bb),
n_pairs_ab=int(result.n_pairs_ab),
converged=bool(result.converged),
localization=localization,
cderi_imaginary_residual=reference.cderi_imaginary_residual,
cderi_symmetry_residual=reference.cderi_symmetry_residual,
matrix_imaginary_residual=reference.matrix_imaginary_residual,
solver_result=result,
hf_result=hf_result,
localization_result=localization_result,
backend=(
"aiccm2026dev-b-dlpno-ump2" if local else "aiccm2026dev-b-ri-ump2"
),
)
[docs]
def run_aiccm2026dev_b_ump2(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
ump2_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUMP2Result:
"""Run canonical finite-torus RI-UMP2 in the exact PNO limit."""
from .dlpno.ump2 import DLPNOUMP2Options
selected = (
DLPNOUMP2Options(localise="none", tcut_pno=0.0, tcut_pairs=0.0)
if ump2_options is None
else ump2_options
)
if selected.localise != "none" or selected.tcut_pno != 0.0:
raise ValueError(
"canonical aiccm2026dev-b UMP2 requires localise='none' and tcut_pno=0"
)
return _run_aiccm2026dev_b_ump2(
system,
basis,
lattice_extension,
options,
ump2_options=selected,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_dlpno_ump2(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
ump2_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUMP2Result:
"""Run PBC-localized, pair-natural-orbital UMP2 on the B torus."""
from .dlpno.ump2 import DLPNOUMP2Options
selected = (
DLPNOUMP2Options(localise="wannier", tcut_pno=1e-8, tcut_pairs=0.0)
if ump2_options is None
else ump2_options
)
return _run_aiccm2026dev_b_ump2(
system,
basis,
lattice_extension,
options,
ump2_options=selected,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
def _run_aiccm2026dev_b_uccsd(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None,
options: PeriodicRHFOptions | None,
*,
cc_options: object,
with_triples: bool,
aux_basis: str | None,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> AICCM2026DevBUCCSDResult:
from .dlpno.uccsd import DLPNOUCCSDPilotOptions, run_dlpno_uccsd_pilot
if not isinstance(cc_options, DLPNOUCCSDPilotOptions):
raise TypeError("cc_options must be a DLPNOUCCSDPilotOptions instance")
if cc_options.localise not in {"none", "wannier", "iao"}:
raise ValueError(
"aiccm2026dev-b UCCSD requires PBC-safe localise='wannier' or "
"'iao', or exact-limit localise='none'"
)
extension = _resolve_lattice_extension(system, lattice_extension)
aux_name = aux_basis or default_aux_for(basis.name)
reference, hf_result = _prepare_b_real_ureference(
system,
basis,
extension,
options,
auxiliary_name=aux_name,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
localized_hf, localization, localization_result = _localize_b_real_ureference(
reference,
cc_options.localise,
hf_result=hf_result,
system=system,
primitive_basis=basis,
mesh=extension,
)
solver_options = replace(
cc_options,
localise=("none" if localization == "none" else "external"),
compute_triples=with_triples,
)
with _patched_overlap(reference.hf.overlap):
result = run_dlpno_uccsd_pilot(
reference.molecule,
reference.basis,
localized_hf,
reference.df,
solver_options,
)
n_cells = reference.n_cells
triples = float(result.e_t) if with_triples else 0.0
local = localization != "none" or cc_options.tcut_pno > 0.0
method = "uccsd(t)" if with_triples else "uccsd"
prefix = "dlpno-" if local else ""
return AICCM2026DevBUCCSDResult(
mesh=extension,
n_cyclic_cells=n_cells,
e_hf_per_cell=float(result.e_hf) / n_cells,
e_corr_per_cell=float(result.e_corr) / n_cells,
e_t_per_cell=triples / n_cells,
e_total_per_cell=(float(result.e_hf) + float(result.e_corr) + triples)
/ n_cells,
n_pairs=int(result.n_pairs),
n_iter=int(result.n_iter),
converged=bool(result.converged),
t1_norm=float(result.t1_norm),
avg_pno=float(result.avg_pno),
localization=localization,
with_triples=with_triples,
cderi_imaginary_residual=reference.cderi_imaginary_residual,
cderi_symmetry_residual=reference.cderi_symmetry_residual,
matrix_imaginary_residual=reference.matrix_imaginary_residual,
solver_result=result,
hf_result=hf_result,
localization_result=localization_result,
backend=f"aiccm2026dev-b-{prefix}{method}",
)
def _canonical_ucc_options(options: object | None, *, triples: bool) -> object:
from .dlpno.uccsd import DLPNOUCCSDPilotOptions
selected = (
DLPNOUCCSDPilotOptions(
localise="none",
tcut_pno=0.0,
compute_triples=triples,
)
if options is None
else options
)
if selected.localise != "none" or selected.tcut_pno != 0.0:
raise ValueError(
"canonical aiccm2026dev-b UCCSD requires localise='none' and tcut_pno=0"
)
return selected
[docs]
def run_aiccm2026dev_b_uccsd(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUCCSDResult:
"""Run full-domain finite-torus UCCSD with the O(N^6) oracle."""
return _run_aiccm2026dev_b_uccsd(
system,
basis,
lattice_extension,
options,
cc_options=_canonical_ucc_options(cc_options, triples=False),
with_triples=False,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_uccsd_t(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUCCSDResult:
"""Run full-domain finite-torus UCCSD(T) with the O(N^6) oracle."""
return _run_aiccm2026dev_b_uccsd(
system,
basis,
lattice_extension,
options,
cc_options=_canonical_ucc_options(cc_options, triples=True),
with_triples=True,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
def _run_dlpno_ucc(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None,
options: PeriodicRHFOptions | None,
*,
cc_options: object | None,
triples: bool,
aux_basis: str | None,
lattice_cutoff_bohr: float,
rsgdf_ke_cutoff: float,
gdf_linear_dep_threshold: float,
progress: bool | object | None,
verbose: int | None,
) -> AICCM2026DevBUCCSDResult:
from .dlpno.uccsd import DLPNOUCCSDPilotOptions
selected = (
DLPNOUCCSDPilotOptions(
localise="wannier",
tcut_pno=1e-7,
compute_triples=triples,
)
if cc_options is None
else cc_options
)
return _run_aiccm2026dev_b_uccsd(
system,
basis,
lattice_extension,
options,
cc_options=selected,
with_triples=triples,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_dlpno_uccsd(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUCCSDResult:
"""Run the unrestricted local-PNO CCSD correctness pilot."""
return _run_dlpno_ucc(
system,
basis,
lattice_extension,
options,
cc_options=cc_options,
triples=False,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)
[docs]
def run_aiccm2026dev_b_dlpno_uccsd_t(
system: PeriodicSystem,
basis: BasisSet,
lattice_extension: int | Sequence[int] | None = None,
options: PeriodicRHFOptions | None = None,
*,
cc_options: object | None = None,
aux_basis: str | None = None,
lattice_cutoff_bohr: float = 15.0,
rsgdf_ke_cutoff: float = 200.0,
gdf_linear_dep_threshold: float = 1e-9,
progress: bool | object | None = None,
verbose: int | None = None,
) -> AICCM2026DevBUCCSDResult:
"""Run the unrestricted local-PNO CCSD(T) correctness pilot."""
return _run_dlpno_ucc(
system,
basis,
lattice_extension,
options,
cc_options=cc_options,
triples=True,
aux_basis=aux_basis,
lattice_cutoff_bohr=lattice_cutoff_bohr,
rsgdf_ke_cutoff=rsgdf_ke_cutoff,
gdf_linear_dep_threshold=gdf_linear_dep_threshold,
progress=progress,
verbose=verbose,
)