"""CIF reading: ``CrystalData`` construction, M–M pair extraction, distance shells."""
from __future__ import annotations
import logging
from typing import Dict, List
import numpy as np
from pymatgen.core import Structure
from mag4.constants import DIST_EQUIV_TOL, CrystalData
from mag4.lattice import lattice_matrix
log = logging.getLogger(__name__)
[docs]
def crystal_from_cif(cif_path: str, element: str) -> CrystalData:
"""Build a :class:`CrystalData` directly from a CIF file (no POSCAR needed).
Parameters
----------
cif_path : str
Path to the CIF file.
element : str
Magnetic element symbol (e.g. ``"Gd"``).
Raises
------
ValueError
If ``element`` is not present in the structure.
"""
structure = Structure.from_file(cif_path)
a, b, c = structure.lattice.a, structure.lattice.b, structure.lattice.c
al, be, ga = (structure.lattice.alpha,
structure.lattice.beta,
structure.lattice.gamma)
cell = [a, b, c, al, be, ga]
lat_mat = lattice_matrix(a, b, c, al, be, ga)
lat_vecs = np.array([lat_mat[:, 0], lat_mat[:, 1], lat_mat[:, 2]])
element_order: List[str] = []
for site in structure:
sym = site.specie.symbol
if sym not in element_order:
element_order.append(sym)
groups_dict: Dict[str, list] = {el: [] for el in element_order}
local_idx = {el: 1 for el in element_order}
global_idx = 1
for site in structure:
sym = site.specie.symbol
label = f"{sym}{local_idx[sym]}"
frac = list(site.frac_coords)
groups_dict[sym].append([label] + frac + [0, global_idx])
local_idx[sym] += 1
global_idx += 1
if element not in element_order:
raise ValueError(
f"Element '{element}' not found in CIF.\nAvailable: {element_order}"
)
all_species = [groups_dict[el] for el in element_order]
nb_species = [len(g) for g in all_species]
target_species = groups_dict[element]
return CrystalData(cell, all_species, target_species,
nb_species, global_idx, element_order,
lat_mat, lat_vecs)
[docs]
def get_mm_distances(cif_path: str, element: str, cutoff: float) -> list:
"""Return sorted list of unique M–M pairs within ``cutoff`` Å.
Each entry is a dict with keys ``site_i, site_j, label_i, label_j, image, distance``.
"""
structure = Structure.from_file(cif_path)
m_indices = [i for i, s in enumerate(structure) if s.specie.symbol == element]
if not m_indices:
raise ValueError(
f"No '{element}' sites found.\n"
f"Available: {sorted({s.specie.symbol for s in structure})}"
)
results, seen = [], set()
for idx_i in m_indices:
for nb in structure.get_neighbors(structure[idx_i], cutoff):
if nb.specie.symbol != element:
continue
idx_j = nb.index
image = tuple(int(x) for x in nb.image)
if idx_i == idx_j and image == (0, 0, 0):
continue
# Canonicalise so the smaller atom index comes first; the
# opposite enumeration (atom j sees atom i at -image) maps to
# the same key after this transformation.
if idx_i < idx_j:
key = (idx_i, idx_j, image)
elif idx_i > idx_j:
key = (idx_j, idx_i, tuple(-x for x in image))
else:
# Self-pair (i == j): +image and -image are the same
# physical bond by translation — pick a canonical sign so
# we count it once. Order images lexicographically and
# keep the lexicographically-greater of (image, -image).
neg_image = tuple(-x for x in image)
key = (idx_i, idx_j, max(image, neg_image))
if key in seen:
continue
seen.add(key)
dist = round(structure.get_distance(idx_i, idx_j, jimage=image), 6)
results.append(dict(
site_i=idx_i, site_j=idx_j,
label_i=f"{element}{idx_i}", label_j=f"{element}{idx_j}",
image=image, distance=dist,
))
results.sort(key=lambda d: d["distance"])
return results
[docs]
def group_inequivalent(results: list, tol: float = DIST_EQUIV_TOL) -> list:
"""Cluster M–M pairs into distance shells (J1, J2, …).
Returns a list of dicts ``{label, distance, dist_min, dist_max, multiplicity, pairs}``.
"""
if not results:
return []
shells, current = [], [results[0]]
for pair in results[1:]:
if abs(pair["distance"] - current[-1]["distance"]) <= tol:
current.append(pair)
else:
shells.append(current)
current = [pair]
shells.append(current)
grouped = []
for k, shell in enumerate(shells, 1):
dists = [p["distance"] for p in shell]
grouped.append(dict(
label = f"J{k}",
distance = sum(dists) / len(dists),
dist_min = min(dists),
dist_max = max(dists),
multiplicity = len(shell),
pairs = shell,
))
return grouped
[docs]
def get_distance_shells(cif_path: str, element: str,
cutoff: float,
tol: float = DIST_EQUIV_TOL) -> List[Dict]:
"""Lightweight version of :func:`group_inequivalent` for the magnon pipeline.
Returns ``[{"label": "J1", "distance": d_mean, "n_bonds": k}, …]``.
"""
structure = Structure.from_file(cif_path)
mag_idx = [i for i, s in enumerate(structure) if s.specie.symbol == element]
if not mag_idx:
raise ValueError(f"No '{element}' sites found in CIF.")
raw, seen = [], set()
for i in mag_idx:
for nb in structure.get_neighbors(structure[i], cutoff):
if nb.specie.symbol != element:
continue
j = nb.index
image = tuple(int(x) for x in nb.image)
if i == j and image == (0, 0, 0):
continue
key = (i, j, image)
if key in seen:
continue
seen.add(key)
raw.append(round(structure.get_distance(i, j, jimage=image), 6))
raw.sort()
shells, cur = [], [raw[0]]
for d in raw[1:]:
if abs(d - cur[-1]) <= tol:
cur.append(d)
else:
shells.append(cur)
cur = [d]
shells.append(cur)
return [{"label": f"J{k}",
"distance": sum(sh) / len(sh),
"n_bonds": len(sh)}
for k, sh in enumerate(shells, 1)]