"""Inequivalent magnetic sites and bond table for the multi-sublattice exchange matrix."""
from __future__ import annotations
import logging
from typing import Dict, List, Tuple
import numpy as np
from pymatgen.core import Structure
from mag4.constants import DIST_EQUIV_TOL
log = logging.getLogger(__name__)
[docs]
def get_inequivalent_sites(cif_path: str, element: str) -> List[List[int]]:
"""Group magnetic sites into crystallographically inequivalent classes.
Uses :class:`pymatgen.symmetry.analyzer.SpacegroupAnalyzer` to derive
Wyckoff orbits and returns one list of global site indices per class.
"""
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
structure = Structure.from_file(cif_path)
spa = SpacegroupAnalyzer(structure)
sym_struct = spa.get_symmetrized_structure()
inequiv_groups = []
for idx_group in sym_struct.equivalent_indices:
if structure[idx_group[0]].specie.symbol != element:
continue
assert all(structure[i].specie.symbol == element for i in idx_group), \
f"Mixed-element symmetry group: {idx_group}"
inequiv_groups.append(list(idx_group))
if not inequiv_groups:
avail = sorted({s.specie.symbol for s in structure})
raise ValueError(
f"No inequivalent {element} sites found. Available: {avail}")
inequiv_groups.sort(key=lambda g: tuple(structure[g[0]].frac_coords))
log.info("Space group: %s", spa.get_space_group_symbol())
log.info("Found %d inequivalent %s site(s):", len(inequiv_groups), element)
for k, grp in enumerate(inequiv_groups):
frac = structure[grp[0]].frac_coords
log.info(" Sublattice %d: multiplicity=%d "
"ref frac=(%.5f, %.5f, %.5f) site indices=%s",
k, len(grp), frac[0], frac[1], frac[2], grp)
return inequiv_groups
[docs]
def get_sublattice_bonds(cif_path: str, element: str,
cutoff: float,
shells: List[Dict],
j_values: Dict[str, float],
tol: float = DIST_EQUIV_TOL) -> Dict:
"""Build the bond table for the full ``N_mag × N_mag`` primitive-cell exchange matrix.
Every magnetic atom is treated as its own sublattice. The reciprocal-space
exchange matrix is
.. math:: J_{αβ}(q) = Σ_{\\text{images}} J_n · \\exp(i\\,q · r_{αβ,\\text{image}})
Returns
-------
dict
Keys: ``n_sub, site_labels, bonds, inequiv_groups, global_to_wyck,
mag_idx, cif_path, orig_rec_matrix, prim_fracs, lat_matrix``.
"""
structure = Structure.from_file(cif_path)
mag_idx = [i for i, s in enumerate(structure) if s.specie.symbol == element]
n_sub = len(mag_idx)
if n_sub == 0:
raise ValueError(f"No {element} sites found in CIF.")
global_to_local = {gi: li for li, gi in enumerate(mag_idx)}
inequiv_groups = get_inequivalent_sites(cif_path, element)
global_to_wyck: Dict[int, int] = {}
for wyck_idx, grp in enumerate(inequiv_groups):
for gi in grp:
global_to_wyck[gi] = wyck_idx
wyck_count = [0] * len(inequiv_groups)
site_labels = []
for gi in mag_idx:
wyck_idx = global_to_wyck[gi]
letter = chr(ord('a') + wyck_count[wyck_idx])
wyck_count[wyck_idx] += 1
site_labels.append(f"{element}{wyck_idx+1}{letter}")
log.info("Full exchange matrix: %d x %d (%d magnetic atoms)",
n_sub, n_sub, n_sub)
for k, grp in enumerate(inequiv_groups):
ref = structure[grp[0]].frac_coords
log.info(" Wyckoff %d (%s): %d sites ref=(%.4f,%.4f,%.4f)",
k+1, element, len(grp), *ref)
dist_to_J: Dict[float, float] = {}
for sh in shells:
lbl = sh["label"]
if lbl in j_values:
dist_to_J[sh["distance"]] = j_values[lbl]
bonds: Dict[Tuple[int, int], list] = {}
for local_alpha, global_i in enumerate(mag_idx):
cart_i = structure[global_i].coords
for nb in structure.get_neighbors(structure[global_i], cutoff):
if nb.specie.symbol != element:
continue
global_j = nb.index
if global_j not in global_to_local:
continue
local_beta = global_to_local[global_j]
image = tuple(int(x) for x in nb.image)
dist = round(structure.get_distance(
global_i, global_j, jimage=image), 6)
J_val = None
for d_ref, j_ref in dist_to_J.items():
if abs(dist - d_ref) <= tol:
J_val = j_ref
break
if J_val is None:
continue
frac_j = structure[global_j].frac_coords + np.array(image)
cart_j = structure.lattice.get_cartesian_coords(frac_j)
vec = cart_j - cart_i
key = (local_alpha, local_beta)
if key not in bonds:
bonds[key] = []
bonds[key].append((vec, J_val))
log.info("Bond pairs (alpha,beta): %d (out of %d possible)",
len(bonds), n_sub * n_sub)
log.info("Total bond entries: %d", sum(len(v) for v in bonds.values()))
return {"n_sub": n_sub,
"site_labels": site_labels,
"bonds": bonds,
"inequiv_groups": inequiv_groups,
"global_to_wyck": global_to_wyck,
"mag_idx": mag_idx,
"cif_path": cif_path,
"orig_rec_matrix": structure.lattice.reciprocal_lattice.matrix,
"prim_fracs": np.array([structure[i].frac_coords
for i in mag_idx]),
"lat_matrix": structure.lattice.matrix}