"""Build (na × nb × nc) supercells from a :class:`CrystalData`."""
from __future__ import annotations
import logging
from typing import List, Optional, Sequence, Tuple
import numpy as np
from mag4.constants import CrystalData
from mag4.dimers import can_assign_unique_dimers, find_dimers
from mag4.lattice import lattice_matrix
log = logging.getLogger(__name__)
[docs]
def expand_supercell(crystal: CrystalData, na: int, nb: int, nc: int) -> CrystalData:
"""Expand a :class:`CrystalData` by ``(na × nb × nc)``.
Returns the original instance unchanged when na = nb = nc = 1.
"""
if na == nb == nc == 1:
return crystal
log.info("Building %d×%d×%d supercell ...", na, nb, nc)
a0, b0, c0, al, be, ga = crystal.cell
new_a, new_b, new_c = na*a0, nb*b0, nc*c0
new_cell = [new_a, new_b, new_c, al, be, ga]
new_lat_mat = lattice_matrix(new_a, new_b, new_c, al, be, ga)
new_lat_vecs = np.array([new_lat_mat[:, 0],
new_lat_mat[:, 1],
new_lat_mat[:, 2]])
target_group_idx = crystal.all_species.index(crystal.target_species)
new_all_species: List[list] = []
total = 1
for group_idx, group in enumerate(crystal.all_species):
el = crystal.elements[group_idx]
new_group: list = []
for ia in range(na):
for ib_ in range(nb):
for ic in range(nc):
for orig in group:
xf = (orig[1] + ia) / na
yf = (orig[2] + ib_) / nb
zf = (orig[3] + ic) / nc
orig_local = orig[0].replace(el, "")
label = f"{el}{orig_local}_{ia}{ib_}{ic}"
new_group.append([label, xf, yf, zf, 0, total])
total += 1
new_all_species.append(new_group)
new_nb_species = [len(g) for g in new_all_species]
new_target = new_all_species[target_group_idx]
log.info("Supercell: %d magnetic atoms (was %d)",
len(new_target), len(crystal.target_species))
return CrystalData(new_cell, new_all_species, new_target,
new_nb_species, total, crystal.elements,
new_lat_mat, new_lat_vecs)
[docs]
def shortest_periodic_image(lat_mat: np.ndarray, max_n: int = 3) -> float:
"""Shortest non-zero supercell lattice vector length.
Iterates over ``(n_a, n_b, n_c) ∈ [-max_n, max_n]^3 \\ {(0,0,0)}``
and returns the minimum ``|n_a a + n_b b + n_c c|``. For cubic /
orthorhombic / hexagonal cells this reduces to ``min(|a|,|b|,|c|)``;
for general triclinic cells the enumeration catches the short
diagonal combinations.
Parameters
----------
lat_mat
``3×3`` matrix whose **columns** are the Cartesian lattice
vectors ``a, b, c`` (mag4's ``CrystalData.lattice_matrix``
convention).
max_n
Range of integer coefficients to try. ``3`` is enough for
all standard Bravais types — only severely skewed triclinic
lattices could need more, and mag4 doesn't build those.
"""
a, b, c = lat_mat[:, 0], lat_mat[:, 1], lat_mat[:, 2]
shortest = float("inf")
for na in range(-max_n, max_n + 1):
for nb in range(-max_n, max_n + 1):
for nc in range(-max_n, max_n + 1):
if na == 0 and nb == 0 and nc == 0:
continue
v = na * a + nb * b + nc * c
d = float(np.linalg.norm(v))
if d < shortest:
shortest = d
return shortest
[docs]
def find_dimer_direction(
conv_lat_mat: np.ndarray,
target_species: List,
distance: float,
*,
dist_tol: float = 1e-2,
max_image: int = 1,
) -> Optional[np.ndarray]:
"""Return a representative unit vector along a J-shell dimer in the
conventional cell.
For each pair ``(i, j)`` of magnetic atoms in ``target_species`` and
each periodic-image translation ``T`` within ``[-max_image, max_image]^3``
of the conv lattice, returns the first ``(r_j + T − r_i)/|.|`` whose
length matches ``distance`` within ``dist_tol``. ``None`` if no such
pair exists in the conv cell (rare: would mean the J shell only
appears across multiple conv-cell copies in the supercell).
"""
a, b, c = conv_lat_mat[:, 0], conv_lat_mat[:, 1], conv_lat_mat[:, 2]
positions = []
for atom in target_species:
frac = np.array([float(atom[1]), float(atom[2]), float(atom[3])])
cart = frac[0] * a + frac[1] * b + frac[2] * c
positions.append(cart)
n = len(positions)
rng = list(range(-max_image, max_image + 1))
for i in range(n):
for j in range(n):
for ni in rng:
for nj in rng:
for nk in rng:
if i == j and ni == 0 and nj == 0 and nk == 0:
continue
T = ni * a + nj * b + nk * c
v = positions[j] + T - positions[i]
d = float(np.linalg.norm(v))
if abs(d - distance) < dist_tol:
return v / d
return None
[docs]
def find_all_dimer_directions(
conv_lat_mat: np.ndarray,
target_species: List,
distance: float,
*,
dist_tol: float = 1e-2,
max_image: int = 1,
) -> List[np.ndarray]:
"""All distinct unit bond directions of a J shell in the conventional cell.
Like :func:`find_dimer_direction` but returns **every** symmetry-equivalent
direction (deduplicated by sign, since ``±u`` give the same periodic-image
distance). A single coupling can run along several directions at once —
e.g. a honeycomb J1 has three — and in an *anisotropic* supercell they
isolate differently, so the isolation check must test all of them, not just
one representative. Returns ``[]`` if the shell has no in-cell pair.
"""
a, b, c = conv_lat_mat[:, 0], conv_lat_mat[:, 1], conv_lat_mat[:, 2]
positions = [float(at[1]) * a + float(at[2]) * b + float(at[3]) * c
for at in target_species]
rng = range(-max_image, max_image + 1)
out: List[np.ndarray] = []
for i in range(len(positions)):
for j in range(len(positions)):
for ni in rng:
for nj in rng:
for nk in rng:
if i == j and ni == 0 and nj == 0 and nk == 0:
continue
v = positions[j] + ni * a + nj * b + nk * c - positions[i]
d = float(np.linalg.norm(v))
if abs(d - distance) >= dist_tol:
continue
u = v / d
if not any(np.allclose(u, s, atol=1e-3)
or np.allclose(u, -s, atol=1e-3) for s in out):
out.append(u)
return out
[docs]
def check_four_state_isolation(
conv_lat_mat: np.ndarray,
multipliers: Tuple[int, int, int],
couplings: List[tuple],
target_species: Optional[List] = None,
*,
dimer_vectors: Optional[dict] = None,
rel_tol: float = 1e-3,
max_multiplicity: Optional[int] = None,
dimer_multiplicity: Optional[dict] = None,
) -> dict:
"""Per-J **algebraic** isolation check for the four-state method.
A Jk dimer ``(i, j)`` is *isolable* in the chosen supercell when its
four-state combination ``E_uu + E_dd − E_ud − E_du`` returns **only** Jk.
Flipping ``(i, j)`` touches one bond for every periodic image of ``j`` that
``i`` couples to (all other atoms are a fixed bath that cancels), so the
combination equals ``4·(Σ over those bonds of their J)``. Therefore:
* **clean** — none of those images sits at a *different* coupling distance
``Jm (m≠k)``. If one does, the formula mixes Jk with Jm and the dimer is
contaminated (``cross`` lists the offending Jm). Cleanliness is the whole
criterion — this is the "check the cancellation of the other dimers in the
4-state equation" rule.
* **multiplicity** ``M`` — the number of images at ``d_k`` itself (``≥ 1``,
the bond included). When a compact cell wraps the bond onto its own image
through the *same* Jk (e.g. an axial bond along a doubled axis links ``±a``,
``M = 2``), the combination yields ``M·Jk·S²`` and the extractor divides by
the self-coefficient ``4·M`` instead of the textbook ``4``.
This is purely a property of which images fall at which coupling distance,
so it is frame-independent and handles oblique cells. It replaces the older
geometric surrogates — bounding box ``Σ|u_i|L_i``, per-axis ``n_i·L_i ≥ 3·d``,
and the "3 copies of ``d`` along the bond" chain rule — all of which demanded
a *fully isolated* (``M = 1``) dimer and so over-sized the cell. The
algebraic rule accepts the smallest cell in which the cross terms cancel,
multiplicity and all.
Worked examples (``--cutoff`` giving the J set):
tetragonal La₂CuO₄ (J1∥a,b; J2∥[110]) → **2×2×1** (J1 ``M=2`` → /8, J2
``M=4`` → /16); CsVI₃ (J1, J2 both ∥c) → **1×1×2** (J2 ``M=2``), J1 alone
→ **1×1×1** (``M=2``); SrFeO₂ (J1,J2,J3) → **2×2×3** (J2,J3 ``M=2``).
The bond directions come from ``target_species`` (all symmetry-equivalent
directions) or a single ``dimer_vectors`` entry; with neither, a shell has
no testable direction and is passed optimistically (``M=1``).
``max_multiplicity`` (opt-in) caps the accepted dimer multiplicity: a shell
is accepted only when it is clean **and** its multiplicity ``M`` (the number
of Jk images linking the dimer atoms — the ``4·M`` divisor) does not exceed
the cap. ``max_multiplicity=1`` demands a *fully* isolated dimer (textbook
divisor 4) — the geometric-style behaviour exposed on the CLI as
``--full-isolation`` — and grows the cell (La₂CuO₄ 2×2×1 → 3×3×1, CsVI₃
1×1×2 → 1×1×3, SrFeO₂ 2×2×3 → 3×3×3) in exchange for larger, less
noise-sensitive cells; intermediate caps (e.g. ``2``) trade cell size
against the largest divisor. A clean shell whose ``M`` exceeds the cap is
reported ``ok=False`` with an *empty* ``cross`` (the failure is high
multiplicity, not contamination).
``dimer_multiplicity`` (post-selection only) maps each coupling to the
multiplicity of the *actually selected* dimer — the exact ``4·M`` the report
divides by, counted from that dimer's supercell images. When given it
overrides the direction-based worst case for clean shells, so the displayed
multiplicity and the cap decision match the formula exactly. The supercell
*search* omits it (no dimer is chosen yet) and stays pessimistic.
Returns a dict with keys:
* ``axes`` — tuple ``(L_a, L_b, L_c)`` of supercell axis lengths.
* ``d_max`` — longest coupling distance (the in-range threshold).
* ``per_shell`` — list of ``{name, d, dimer_dir, multiplicity, cross, ok}``
dicts, one per coupling. ``dimer_dir`` is the chosen clean direction (the
one of smallest multiplicity), ``cross`` the sorted list of contaminating
Jm (empty when clean), ``ok = (cross is empty)``. A shell with several
symmetry-equivalent directions is OK if *any* one is clean (mag4 selects
that representative dimer).
* ``ok`` — ``True`` iff every shell is clean.
* ``failing`` — ``(name, d)`` tuples for contaminated shells.
* ``recommended`` — smallest ``(na, nb, nc)`` (≥ ``multipliers``, by atom
count) for which every shell is clean.
"""
na, nb, nc = multipliers
a_vec = np.asarray(conv_lat_mat[:, 0], dtype=float)
b_vec = np.asarray(conv_lat_mat[:, 1], dtype=float)
c_vec = np.asarray(conv_lat_mat[:, 2], dtype=float)
conv_lens = (float(np.linalg.norm(a_vec)),
float(np.linalg.norm(b_vec)),
float(np.linalg.norm(c_vec)))
axes = (na * conv_lens[0], nb * conv_lens[1], nc * conv_lens[2])
L_max = max(axes)
shell_dist = [(c[0], c[1]) for c in couplings]
d_max = max((d for _, d in shell_dist), default=0.0)
# Non-zero supercell translations: ±2 cells along each axis is enough to
# reach the nearest image of any single bond, even in an oblique cell.
_shifts = np.array([(i, j, k)
for i in (-2, -1, 0, 1, 2)
for j in (-2, -1, 0, 1, 2)
for k in (-2, -1, 0, 1, 2)
if (i, j, k) != (0, 0, 0)], dtype=float)
def _dirs_for(name: str, d: float) -> List[np.ndarray]:
"""Every distinct bond direction of a shell. ALL symmetry-equivalent
directions must be isolated (a honeycomb J1 runs three ways at once, and
in an anisotropic supercell they isolate differently)."""
if target_species is not None:
dirs = find_all_dimer_directions(conv_lat_mat, target_species, d)
if dirs:
return dirs
if dimer_vectors is not None and name in dimer_vectors:
v = dimer_vectors[name]
if v is not None:
v = np.asarray(v, dtype=float)
nv = float(np.linalg.norm(v))
if nv > 0:
return [v / nv]
return []
all_d = [dd for _, dd in shell_dist]
DTOL = 1e-2 # Å: match an image distance to a coupling shell
# All non-zero translations PLUS T=0 (the bond itself).
_trans0 = np.vstack([np.zeros((1, 3)), _shifts])
def _dir_clean(d_k: float, u: np.ndarray, mult: Tuple[int, int, int]):
"""``(clean, multiplicity, cross)`` for one Jk dimer direction
(bond ``b = d_k·u``).
The four-state combination of the dimer ``(i, j)`` gains a term for every
bond linking ``i`` to an image of ``j``. ``clean`` is ``True`` when none
of those images sits at a *different* coupling distance (Jm, m≠k), so the
formula returns only Jk; ``cross`` lists the offending Jm otherwise.
``multiplicity`` counts the images at ``d_k`` itself (≥1, incl. the
bond): the combination then yields ``multiplicity·Jk·S²`` and the
extractor divides by it.
"""
b = d_k * u
sc = np.array([mult[0] * a_vec, mult[1] * b_vec, mult[2] * c_vec])
m_k, cross = 0, set()
for r in np.linalg.norm(b + (_trans0 @ sc), axis=1):
if r > d_max + DTOL:
continue
name_near, d_near = min(shell_dist, key=lambda nd: abs(nd[1] - r))
if abs(r - d_near) < DTOL:
if abs(d_near - d_k) < DTOL:
m_k += 1
else:
cross.add(name_near)
return (not cross), m_k, cross
def _shell(name: str, d: float, dirs: List[np.ndarray],
mult: Tuple[int, int, int]):
"""Cleanliness/multiplicity of a shell. Returns
``(ok, multiplicity, dir, cross)``.
Isolation (cleanliness) is optimistic — a shell is clean if SOME
symmetry-equivalent dimer is clean, since mag4 only needs one clean
representative. Multiplicity, however, is reported as the WORST
(largest) case over the clean directions: in an anisotropic supercell
the equivalent dimers carry different multiplicities (NiO J2 runs M=1
along the long axes but M=2 along a short one), and the greedy dimer
assignment may pick any of them — so ``--max-multiplicity`` /
``--full-isolation`` must hold for the worst one, or an
over-multiplicity dimer slips through. A clean shell whose multiplicity
exceeds the cap is returned ``ok=False`` with an empty ``cross`` (the
failure is high multiplicity, not contamination)."""
if not dirs: # no direction → optimistic pass
return True, 1, None, set()
per = [(_dir_clean(d, u, mult), u) for u in dirs]
clean = [(cm[1], u) for cm, u in per if cm[0] and cm[1] >= 1]
if clean:
# Worst-case multiplicity over the clean equivalent dimers.
m_k, wdir = max(clean, key=lambda x: x[0])
if max_multiplicity is not None and m_k > max_multiplicity:
return False, m_k, wdir, set() # clean but exceeds the cap
return True, m_k, wdir, set()
(cm, u) = min(per, key=lambda v: len(v[0][2])) # fewest contaminants
return False, None, u, cm[2]
dirs_of = {name: _dirs_for(name, d) for name, d in shell_dist}
def _accept(name: str, d: float, dirs, mult: Tuple[int, int, int]) -> bool:
return _shell(name, d, dirs, mult)[0]
per_shell, failing = [], []
for name, d in shell_dist:
ok, m_k, wdir, cross = _shell(name, d, dirs_of[name], (na, nb, nc))
# Post-selection callers know the multiplicity of the *actually chosen*
# dimer (counted from its supercell images) and pass it in. It is the
# exact ``4·M`` the report divides by, so prefer it over the direction-
# based worst case for both the displayed value and the cap decision —
# the search stays pessimistic (it has no chosen dimer yet).
if (dimer_multiplicity is not None and name in dimer_multiplicity
and not cross):
m_k = int(dimer_multiplicity[name])
ok = (max_multiplicity is None) or (m_k <= max_multiplicity)
per_shell.append({"name": name, "d": d, "dimer_dir": wdir,
"multiplicity": m_k, "cross": sorted(cross), "ok": ok})
if not ok:
failing.append((name, d))
# Smallest enlargement (by atom count, then longest axis) for which every
# shell is accepted; searched over a bounded window above the request.
def _search() -> Tuple[int, int, int]:
base = tuple(max(1, m) for m in multipliers)
best, best_key, window = None, None, 6
for da in range(window + 1):
for db in range(window + 1):
for dc in range(window + 1):
cand = (base[0] + da, base[1] + db, base[2] + dc)
if all(_accept(name, d, dirs_of[name], cand)
for name, d in shell_dist):
key = (cand[0] * cand[1] * cand[2], max(cand))
if best_key is None or key < best_key:
best, best_key = cand, key
return best or base
rec = _search() if failing else tuple(max(1, m) for m in multipliers)
return {
"axes": axes,
"L_max": L_max,
"d_max": d_max,
"per_shell": per_shell,
"ok": not failing,
"failing": failing,
"recommended": tuple(rec),
}
[docs]
def find_optimal_supercell(
crystal: CrystalData,
couplings: list,
tol: float,
max_mult: int = 4,
*,
ligand_elements: Optional[Sequence[str]] = None,
bond_cutoff: float = 2.5,
angle_tol: float = 1.0,
ll_cutoff: float = 0.0,
signed_dihedral: bool = False,
isolation: bool = True,
max_multiplicity: Optional[int] = None,
) -> Tuple[Tuple[int, int, int], CrystalData]:
"""Return the smallest ``(na, nb, nc)`` supercell that **resolves and
isolates** every coupling for the four-state method.
Candidates are enumerated over ``[1..max_mult]^3`` and ranked by
``(volume, max(na,nb,nc), (na,nb,nc))`` so that (1,1,1) is tried first
and the most compact / isotropic cells win ties.
Two acceptance criteria are applied to each candidate, in order of size:
1. **Resolve** — every coupling has a distinct intra-cell dimer pair
(:func:`can_assign_unique_dimers`).
2. **Isolate** (when ``isolation=True``, the default) — every J shell is
also four-state-*clean* (:func:`check_four_state_isolation`): no periodic
image of the partner couples to it through a *different* Jm, so the
four-state combination returns only Jk (mixed in with its own multiplicity
``M``, which the extractor divides out as ``4·M``). This is the algebraic
cancellation rule — much weaker than demanding a fully isolated dimer — so
the chosen cell is the smallest in which the cross terms vanish.
The smallest candidate passing both is returned. If none isolates within
``max_mult`` the function **falls back** to the smallest *resolving* cell
and logs a warning; the CLI's post-selection isolation check then **errors**
(the four-state method is only valid for an isolated dimer), telling the
user to raise ``--max-supercell`` or lower ``--cutoff``. Pass
``isolation=False`` to restore the legacy resolve-only selection, or
``max_multiplicity=N`` to cap the accepted dimer multiplicity — ``1`` demands
a *fully* isolated dimer (divisor 4, the old geometric "3 copies of d" rule),
larger values trade cell size against the largest ``4·M`` divisor.
When ``ligand_elements`` is supplied the geometric fingerprint of each
candidate pair must also match the target — sub-shells like ``J1a`` /
``J1b`` are then resolved as distinct couplings, which may require a
larger supercell than distance-only mode.
Raises
------
RuntimeError
When no candidate up to ``max_mult`` in each direction yields a
complete, conflict-free dimer assignment. The message lists the
couplings that remained unresolvable.
"""
geom_mode = bool(ligand_elements)
candidates: List[Tuple[int, int, int]] = [
(na, nb, nc)
for na in range(1, max_mult + 1)
for nb in range(1, max_mult + 1)
for nc in range(1, max_mult + 1)
]
candidates.sort(key=lambda t: (t[0] * t[1] * t[2], max(t), t))
best_partial: Tuple[Tuple[int, int, int], dict] | None = None
resolved_fallback: Tuple[Tuple[int, int, int], CrystalData] | None = None
# Dimer directions are a property of the conv-cell geometry (independent of
# the multipliers), so derive them once and reuse for every candidate's
# isolation check.
conv_lat = np.asarray(crystal.lattice_matrix, dtype=float)
conv_targets = list(crystal.target_species)
dimer_dirs = None
if isolation:
dimer_dirs = {c[0]: find_dimer_direction(conv_lat, conv_targets, c[1])
for c in couplings}
for na, nb, nc in candidates:
sc = expand_supercell(crystal, na, nb, nc) if (na, nb, nc) != (1, 1, 1) else crystal
sc_struct = None
if geom_mode:
from mag4.geometry import crystal_to_pymatgen
sc_struct = crystal_to_pymatgen(sc)
dimers = find_dimers(
sc.target_species, sc.lattice_matrix, couplings, tol,
verbose=False,
structure=sc_struct,
ligand_elements=ligand_elements,
bond_cutoff=bond_cutoff,
angle_tol=angle_tol,
ll_cutoff=ll_cutoff,
signed_dihedral=signed_dihedral,
)
if can_assign_unique_dimers(dimers, couplings):
# Resolves every coupling. Remember the smallest such cell as the
# fallback in case nothing isolates within max_mult.
if resolved_fallback is None:
resolved_fallback = ((na, nb, nc), sc)
if isolation:
# No-self-image-folding guard: every coupling must be shorter
# than the cell's shortest lattice translation. Otherwise a
# dimer bond spans a full cell — its atoms link to their own
# images — and couplings BEYOND the cutoff fold onto the dimer
# and contaminate the four-state formula (e.g. La₂CuO₄'s 7.56 Å
# collinear J3 leaking into J2 in a 1×2×1 cell). This is the same
# criterion energy-mapping uses, so the two methods agree on the
# cell.
d_max_c = max((c[1] for c in couplings), default=0.0)
if shortest_periodic_image(sc.lattice_matrix) <= d_max_c + 1e-3:
continue # a coupling folds onto a self-image → grow
iso = check_four_state_isolation(
conv_lat, (na, nb, nc), couplings,
target_species=conv_targets, dimer_vectors=dimer_dirs,
max_multiplicity=max_multiplicity,
)
if not iso["ok"]:
continue # resolved but not isolated → keep growing
return (na, nb, nc), sc
# Track the candidate with the most couplings found, for diagnostics.
found_count = sum(1 for c in couplings if dimers.get(c[0]))
if best_partial is None or found_count > best_partial[1].get("_count", -1):
best_partial = ((na, nb, nc), {"_count": found_count, "dimers": dimers})
# No candidate isolated every shell. Fall back to the smallest cell that at
# least resolves them, and warn (the CLI prints the per-shell detail).
if resolved_fallback is not None:
if isolation:
log.warning(
"No supercell in [1..%d]^3 ISOLATES every J shell (every "
"four-state combination free of other-Jm contamination). Falling "
"back to the smallest cell that resolves them: %dx%dx%d — some "
"shells will be contaminated by periodic-image coupling. Raise "
"--max-supercell to search larger cells, or pass --supercell.",
max_mult, *resolved_fallback[0])
return resolved_fallback
# Nothing even resolved → assemble a helpful error.
missing: list[str] = []
if best_partial is not None:
d = best_partial[1]["dimers"]
missing = [c[0] for c in couplings if not d.get(c[0])]
raise RuntimeError(
f"No supercell in [1..{max_mult}]^3 yields a distinct intra-cell pair "
f"for every coupling. Couplings still missing at the best candidate "
f"{best_partial[0] if best_partial else '(none)'}: "
f"{', '.join(missing) if missing else 'unknown'}. "
f"Re-run with --supercell N1 N2 N3 to override, or relax --cutoff."
)
[docs]
def find_supercell_for_enumeration(
crystal: CrystalData,
couplings: list,
max_mult: int = 4,
*,
tol: float = 1e-3,
) -> Tuple[Tuple[int, int, int], CrystalData]:
"""Smallest supercell with no coupling folding onto a periodic self-image.
This is the **energy-mapping** sizing criterion (distinct from the
four-state :func:`find_optimal_supercell`). The all-orders enumeration
needs every coupling distance ``d`` to stay *strictly below* the shortest
lattice translation; otherwise ``s_i·s_j`` collapses to ``s_i² = +1`` in
every order, the coefficient becomes a constant absorbed into ``E0``, and
the coupling is indeterminate in the fit.
Returns the smallest ``(na, nb, nc)`` (size-ordered) for which
``shortest_periodic_image(super) > d_max + tol``. If none within
``max_mult`` qualifies, returns the largest candidate (which has the
longest lattice translations, i.e. the fewest folds) and warns.
"""
d_max = max(d for _, d in couplings)
candidates = sorted(
((na, nb, nc)
for na in range(1, max_mult + 1)
for nb in range(1, max_mult + 1)
for nc in range(1, max_mult + 1)),
key=lambda t: (t[0] * t[1] * t[2], max(t), t),
)
largest: Tuple[Tuple[int, int, int], CrystalData] | None = None
for na, nb, nc in candidates:
sc = expand_supercell(crystal, na, nb, nc) if (na, nb, nc) != (1, 1, 1) else crystal
largest = ((na, nb, nc), sc)
if shortest_periodic_image(sc.lattice_matrix) > d_max + tol:
return (na, nb, nc), sc
log.warning(
"No supercell in [1..%d]^3 keeps every coupling below the shortest "
"lattice translation (d_max = %.3f Å): some shells fold onto "
"self-images and will be indeterminate in the fit. Using the largest "
"candidate %dx%dx%d — pass an explicit --supercell to choose another.",
max_mult, d_max, *largest[0])
return largest