"""Dimer search for the four-state method (direct intra-cell Cartesian distances)."""
from __future__ import annotations
import logging
import math
from typing import Dict, List, Optional, Sequence
import numpy as np
from scipy.spatial.distance import cdist
from mag4.constants import DIST_EQUIV_TOL, Dimer
from mag4.lattice import get_cartesian_coords
log = logging.getLogger(__name__)
[docs]
def check_coupling_spacing(couplings: list, user_tol: float) -> float:
"""Inspect spacing between consecutive J distances and return a safe dimer-search tolerance.
Warns when:
* ``--tol`` is larger than half the minimum spacing (two J shells could share a dimer).
* Two consecutive J values are closer than ``2*DIST_EQUIV_TOL`` (probably the same shell).
When two adjacent couplings are intentionally at the **same** distance
(e.g. sub-shells ``J1a`` / ``J1b`` produced by geometry-aware grouping —
distinguished by their bridging-angle fingerprint, not by distance), the
zero gap between them is **not** a noise issue: it is dropped from the
warning and from the ``safe_tol`` computation so the dimer search can
still use a non-zero distance tolerance. The angle filter in
:func:`find_dimers` keeps the sub-shells apart.
"""
if len(couplings) < 2:
return user_tol
W = 72
def _diff_by_geometry(i: int) -> bool:
"""True if entries i and i+1 share a distance but differ by fingerprint."""
a, b = couplings[i], couplings[i + 1]
if abs(a[1] - b[1]) > 1e-6:
return False
fa = a[2] if len(a) >= 3 else None
fb = b[2] if len(b) >= 3 else None
return (fa is not None) and (fb is not None) and (fa != fb)
# Compute spacings; skip same-distance sub-shells (different fingerprints).
spacings: list[float] = []
for i in range(len(couplings) - 1):
if _diff_by_geometry(i):
continue
spacings.append(couplings[i + 1][1] - couplings[i][1])
if not spacings:
# All adjacent pairs are geom-split sub-shells: no useful spacing constraint.
return user_tol
min_sp = min(spacings)
warned = False
for i in range(len(couplings) - 1):
if _diff_by_geometry(i):
continue
sp = couplings[i + 1][1] - couplings[i][1]
if sp < 2 * DIST_EQUIV_TOL:
if not warned:
print()
print("!" * W)
print(" WARNING: some J values are very close together.")
print(" They may be the same physical shell split by")
print(" numerical noise. Consider merging them with a")
print(" larger --equiv-tol value.")
print("!" * W)
warned = True
n1, d1 = couplings[i][0], couplings[i][1]
n2, d2 = couplings[i + 1][0], couplings[i + 1][1]
print(f" {n1}={d1:.4f} Å and {n2}={d2:.4f} Å"
f" (gap = {sp*1000:.1f} mÅ)")
safe_tol = min_sp / 2.0
if safe_tol < user_tol:
print()
print("!" * W)
print(f" WARNING: --tol {user_tol:.5f} Å is larger than half the")
print(f" minimum J spacing ({min_sp:.4f} Å).")
print(f" Automatically using safe tolerance = {safe_tol:.5f} Å")
print(f" to prevent adjacent J values from sharing the same dimer.")
print("!" * W)
print()
return safe_tol
return user_tol
[docs]
def find_dimers(atoms: list, lat_mat: np.ndarray,
couplings: list, tol: float,
verbose: bool = True,
*,
structure=None,
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) -> Dict[str, List[Dimer]]:
"""Find all intra-cell atom pairs for each coupling distance ± ``tol``.
Set ``verbose=False`` to silence info logs and the missing-coupling
warning block; useful when scanning many trial supercells.
Design principle
----------------
The four-state method requires that the two atoms of a dimer be DISTINCT
sites in the working cell so that they can be assigned independent spins.
Using periodic images to "find" a coupling would conflate atom *j* with
its image *j + (h, k, l)*, which is already a different atom in the supercell.
The correct workflow is therefore:
1. Build a supercell large enough that every coupling distance d_Jn appears
as a genuine intra-cell pair (direct Cartesian distance).
2. Use plain :func:`scipy.spatial.distance.cdist` here — no wrapping, no
image search.
If a coupling is not found as a direct pair the cell is too small along
one axis. The function diagnoses which axis and tells the user exactly
which ``--supercell`` flag to try.
Angle-aware mode
----------------
When ``structure`` (a :class:`pymatgen.core.Structure` of the *same*
supercell, with the magnetic atoms placed at indices ``0..N_mag-1`` —
see :func:`mag4.geometry.crystal_to_pymatgen`) and ``ligand_elements`` are
supplied, each ``couplings`` entry must be a triple
``[label, distance, fingerprint]``. A candidate pair is then accepted
only when both its distance **and** its M–L–M bridging-angle fingerprint
match the target. This lets the four-state pipeline distinguish J
shells that are degenerate in distance but route super-exchange through
different bridging geometries.
"""
coords = get_cartesian_coords(atoms, lat_mat)
n = len(atoms)
dist_matrix = cdist(coords, coords)
cell = [np.linalg.norm(lat_mat[:, k]) for k in range(3)]
dimers: Dict[str, List[Dimer]] = {}
W = 72
geom_mode = (structure is not None) and bool(ligand_elements)
fp_cache: Dict[tuple, tuple] = {}
def _fingerprint(ia: int, ib: int):
key = (ia, ib) if ia < ib else (ib, ia)
if key not in fp_cache:
from mag4.geometry import compute_pair_fingerprint
fp_cache[key] = compute_pair_fingerprint(
structure, key[0], key[1], (0, 0, 0),
ligand_elements, bond_cutoff, angle_tol,
ll_cutoff=ll_cutoff,
signed_dihedral=signed_dihedral,
)
return fp_cache[key]
for entry in couplings:
if len(entry) == 2:
coup_name, coup_dist = entry
target_fp = None
else:
coup_name, coup_dist, target_fp = entry[0], entry[1], entry[2]
pairs = []
for i in range(n):
for j in range(i + 1, n):
if abs(dist_matrix[i, j] - coup_dist) >= tol:
continue
if geom_mode and target_fp is not None:
if _fingerprint(i, j) != target_fp:
continue
pairs.append(Dimer(
idx_a=i, idx_b=j,
label_a=atoms[i][0], label_b=atoms[j][0],
distance=dist_matrix[i, j],
coupling_name=coup_name,
))
dimers[coup_name] = pairs
if not verbose:
continue
if pairs:
log.info(" %s (d=%.3f Å, tol=±%.5f Å): %d pair(s) found",
coup_name, coup_dist, tol, len(pairs))
# Per-pair listing is debug-only — it can run to hundreds of lines
# in a large supercell; the count above is the useful summary.
for d in pairs:
log.debug(" %s -- %s (d=%.4f Å)", d.label_a, d.label_b, d.distance)
else:
axis_names = ["a", "b", "c"]
log.warning(" %s (d=%.3f Å): NO direct pair found in this cell.",
coup_name, coup_dist)
need = 2.0 * coup_dist
short = [(axis_names[k], cell[k], math.ceil(need / cell[k]))
for k in range(3) if cell[k] < need]
if short:
print()
print("!" * W)
print(f" WARNING: {coup_name} (d={coup_dist:.4f} Å) not found as a")
print(f" direct intra-cell pair. Need each cell length > {need:.4f} Å.")
print(f" Short axes in the CURRENT cell:")
suggest = [1, 1, 1]
for ax, length, mult in short:
k = axis_names.index(ax)
suggest[k] = mult
print(f" {ax} = {length:.4f} Å → multiply by {mult}"
f" → {length*mult:.4f} Å")
print(f" Try: --supercell {suggest[0]} {suggest[1]} {suggest[2]}")
print(f" (use --auto-supercell to let mag4 pick the smallest one)")
print("!" * W)
print()
return dimers
[docs]
def can_assign_unique_dimers(dimers: Dict[str, List[Dimer]],
couplings: list) -> bool:
"""Return True iff the greedy assignment used by the four-state pipeline
can give every coupling a distinct (idx_a, idx_b) pair.
Mirrors the assignment logic in :mod:`mag4.cli.fourstate` so the optimal
supercell search rejects cells that find the right distances but cannot
assign them all to non-overlapping site pairs.
"""
used: set = set()
for entry in couplings:
name = entry[0]
candidates = dimers.get(name, [])
chosen = None
for d in candidates:
key = (d.idx_a, d.idx_b)
if key not in used:
chosen = d
break
if chosen is None:
return False
used.add((chosen.idx_a, chosen.idx_b))
return True