"""Bridging-ligand bond-angle analysis (M – L – M)."""
from __future__ import annotations
import logging
from typing import List
import numpy as np
from mag4.constants import DIST_EQUIV_TOL
log = logging.getLogger(__name__)
[docs]
def get_mm_pairs(structure, element: str, cutoff: float):
"""Return sorted list of unique M–M pairs within ``cutoff`` Å.
Mirrors :func:`mag4.io_cif.get_mm_distances` but returns lighter
dicts (no ``label_*`` fields) suited to the bond-angle workflow.
"""
m_indices = [i for i, s in enumerate(structure)
if s.specie.symbol == element]
if not m_indices:
avail = sorted({s.specie.symbol for s in structure})
raise ValueError(f"No '{element}' sites found. Available: {avail}")
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
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: ±image is the same bond by translation.
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,
image=image, distance=dist,
))
results.sort(key=lambda d: d["distance"])
return results
[docs]
def group_into_shells(results, tol: float = DIST_EQUIV_TOL):
"""Cluster M–M pairs into distance shells (J1, J2, …)."""
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),
multiplicity=len(shell),
pairs=shell,
))
return grouped
[docs]
def find_bridging_angles(structure, site_i, site_j, jimage,
ligand_elements, bond_cutoff):
"""Find every ligand L bridging ``site_i`` and ``site_j`` (at ``jimage``).
A ligand L bridges the pair when ``dist(site_i, L) < bond_cutoff`` AND
``dist(site_j + jimage, L) < bond_cutoff``.
Returns a list of dicts sorted by angle.
"""
lattice = structure.lattice
cart_i = structure[site_i].coords
cart_j = structure[site_j].coords + lattice.get_cartesian_coords(jimage)
neighbors_i = structure.get_neighbors(structure[site_i], bond_cutoff)
ligand_nbs = [nb for nb in neighbors_i
if nb.specie.symbol in ligand_elements]
angles = []
for nb in ligand_nbs:
cart_L = nb.coords
dist_Lj = float(np.linalg.norm(cart_L - cart_j))
if dist_Lj > bond_cutoff:
continue
dist_Li = float(np.linalg.norm(cart_L - cart_i))
vec_Li = cart_i - cart_L
vec_Lj = cart_j - cart_L
cos_a = np.dot(vec_Li, vec_Lj) / (np.linalg.norm(vec_Li)
* np.linalg.norm(vec_Lj))
angle = float(np.degrees(np.arccos(np.clip(cos_a, -1.0, 1.0))))
angles.append(dict(
ligand_symbol=nb.specie.symbol,
ligand_index=nb.index,
ligand_image=tuple(int(x) for x in nb.image),
dist_iL=round(dist_Li, 4),
dist_Lj=round(dist_Lj, 4),
angle=round(angle, 2),
))
unique, seen_keys = [], set()
for a in angles:
key = (a["ligand_index"], a["ligand_image"])
if key not in seen_keys:
seen_keys.add(key)
unique.append(a)
unique.sort(key=lambda a: a["angle"])
return unique
[docs]
def find_dihedral_paths(structure, site_i, site_j, jimage,
ligand_elements, bond_cutoff, ll_cutoff):
"""Find every M1–L1–L2–M2 super-super-exchange path between
``site_i`` (home) and ``site_j`` (at ``jimage``).
A path exists when:
* ``L1`` is bonded to ``M1`` (distance ≤ ``bond_cutoff``);
* ``L2`` is bonded to ``M2 + jimage`` (distance ≤ ``bond_cutoff``);
* ``L1`` and ``L2`` are bonded to each other
(distance ≤ ``ll_cutoff``);
* ``L1`` and ``L2`` are **distinct** ligand atoms (otherwise the path
is a direct M–L–M bridge already covered by
:func:`find_bridging_angles`).
The returned dihedral is the **signed** torsion angle of the
four-atom path in degrees, in ``[-180°, +180°]`` following the IUPAC
convention.
Returns a list of dicts, deduplicated on ``(L1 identity, L2 identity)``
and sorted by signed dihedral.
"""
lattice = structure.lattice
cart_i = structure[site_i].coords
cart_j = structure[site_j].coords + lattice.get_cartesian_coords(jimage)
def _candidates(centre_coord, anchor_image):
"""Return list of (ligand_index, absolute_image_tuple, cart_coords)
for ligand sites within ``bond_cutoff`` of ``centre_coord``. The
absolute image is given relative to lattice origin so we can detect
whether two candidates from the two M centres refer to the same
physical atom."""
out = []
# We use Structure.get_neighbors on the anchor site to find ligand
# neighbours. The pymatgen API guarantees this respects PBC.
anchor_idx_image = anchor_image
# We want neighbours relative to the centre coord — but
# Structure.get_neighbors works with site indices, not arbitrary
# points. Use site_i and site_j directly; shift by jimage as needed.
return out
# Ligand neighbours of M1 in its home position.
ligands_i = []
for nb in structure.get_neighbors(structure[site_i], bond_cutoff):
if nb.specie.symbol in ligand_elements:
ligands_i.append(dict(
index=nb.index,
image=tuple(int(x) for x in nb.image),
coords=np.asarray(nb.coords, dtype=float),
))
# Ligand neighbours of M2 (relative to M2's home), shifted by jimage
# so coordinates live in the same Cartesian frame as cart_j.
j_offset = lattice.get_cartesian_coords(jimage)
ligands_j = []
for nb in structure.get_neighbors(structure[site_j], bond_cutoff):
if nb.specie.symbol in ligand_elements:
ligands_j.append(dict(
index=nb.index,
# Effective image = nb.image + jimage (in lattice units).
image=tuple(int(nb.image[k]) + int(jimage[k]) for k in range(3)),
coords=np.asarray(nb.coords, dtype=float) + j_offset,
))
paths = []
seen = set()
for L1 in ligands_i:
for L2 in ligands_j:
# Reject direct-bridge case (same physical ligand atom).
if L1["index"] == L2["index"] and L1["image"] == L2["image"]:
continue
d_LL = float(np.linalg.norm(L1["coords"] - L2["coords"]))
if d_LL > ll_cutoff:
continue
# Signed dihedral M1–L1–L2–M2 (IUPAC convention).
b1 = L1["coords"] - cart_i
b2 = L2["coords"] - L1["coords"] # axis
b3 = cart_j - L2["coords"]
n1 = np.cross(b1, b2)
n2 = np.cross(b2, b3)
if np.linalg.norm(n1) < 1e-9 or np.linalg.norm(n2) < 1e-9:
continue # degenerate (colinear)
b2_n = b2 / np.linalg.norm(b2)
m1 = np.cross(n1, b2_n)
x = float(np.dot(n1, n2))
y = float(np.dot(m1, n2))
dihedral = float(np.degrees(np.arctan2(y, x)))
# arctan2 returns in (-180°, 180°] but tiny y noise can yield
# exactly -180°; that aliases to +180° (same physical angle).
if dihedral <= -179.999:
dihedral = 180.0
# Dedup on (L1 identity, L2 identity) so the same physical
# path isn't counted twice.
key = (L1["index"], L1["image"], L2["index"], L2["image"])
if key in seen:
continue
seen.add(key)
paths.append(dict(
ligand1_symbol=structure[L1["index"]].specie.symbol,
ligand1_index=L1["index"],
ligand1_image=L1["image"],
ligand2_symbol=structure[L2["index"]].specie.symbol,
ligand2_index=L2["index"],
ligand2_image=L2["image"],
dist_iL1=float(np.linalg.norm(L1["coords"] - cart_i)),
dist_L1L2=d_LL,
dist_L2j=float(np.linalg.norm(L2["coords"] - cart_j)),
dihedral=round(dihedral, 2),
))
paths.sort(key=lambda p: p["dihedral"])
return paths
[docs]
def compress_angles(ang_vals: List[float], tol: float = 0.01) -> str:
"""Group similar angles: ``[159.07, 159.07, 159.07]`` → ``'159.07° (x3)'``."""
if not ang_vals:
return ""
parts, i = [], 0
while i < len(ang_vals):
v = ang_vals[i]
count = 1
while i + count < len(ang_vals) and abs(ang_vals[i + count] - v) <= tol:
count += 1
if count > 1:
parts.append(f"{v:.2f}° (x{count})")
else:
parts.append(f"{v:.2f}°")
i += count
return " ".join(parts)
[docs]
def print_shell_summary(groups, element):
"""Brief summary of the J shells."""
W = 72
print("=" * W)
print(f" {element}-{element} coupling shells")
print("=" * W)
print(f" {'Coupling':<10} {'d_mean (Å)':>10} {'Multiplicity':>12}")
print("-" * W)
for g in groups:
print(f" {g['label']:<10} {g['distance']:>10.4f} {g['multiplicity']:>12}")
print("=" * W)
print()
[docs]
def print_angles_for_shell(group, structure, element, ligand_elements,
bond_cutoff, lines_out):
"""Bond-angle table for one J shell. Appends each printed line to ``lines_out``."""
W = 72
label = group["label"]
d_avg = group["distance"]
def out(msg=""):
print(msg)
lines_out.append(msg)
out("=" * W)
out(f" {label} (d = {d_avg:.4f} Å) — "
f"{element}-L-{element} bridging angles")
out("=" * W)
for pair in group["pairs"]:
si, sj = pair["site_i"], pair["site_j"]
img = pair["image"]
dist = pair["distance"]
label_i = f"{element}{si}"
img_str = f"({img[0]:+d} {img[1]:+d} {img[2]:+d})"
label_j = f"{element}{sj}{img_str}"
out(f"\n Pair: {label_i} — {label_j} d = {dist:.4f} Å")
out(f" {'Ligand':<16} {'Image':<16} "
f"{'d(M1-L) Å':>10} {'d(M2-L) Å':>10} {'Angle (°)':>10}")
out(" " + "-" * (W - 4))
angles = find_bridging_angles(
structure, si, sj, img, ligand_elements, bond_cutoff)
if not angles:
out(f" (no bridging {'/'.join(ligand_elements)} atom found "
f"within {bond_cutoff:.2f} Å of both sites)")
else:
for a in angles:
lig_lbl = f"{a['ligand_symbol']}{a['ligand_index']}"
lig_img = "({:+d} {:+d} {:+d})".format(*a["ligand_image"])
out(f" {lig_lbl:<16} {lig_img:<16} "
f"{a['dist_iL']:>10.4f} {a['dist_Lj']:>10.4f} "
f"{a['angle']:>10.2f}")
ang_vals = sorted(a["angle"] for a in angles)
out(f"\n {len(angles)} bridging path(s): "
+ compress_angles(ang_vals))
out()
[docs]
def write_report(lines_out, output_file):
"""Write accumulated report lines to ``output_file``."""
with open(output_file, "w") as f:
f.write("\n".join(lines_out) + "\n")
print(f"Report written to {output_file}")