"""Magnon dispersion via the eigenvalues of J(q) in the primitive cell."""
from __future__ import annotations
import copy
import logging
from typing import Dict, Tuple
import numpy as np
from scipy.optimize import minimize
from mag4.magnon.exchange import build_jq_matrices
log = logging.getLogger(__name__)
[docs]
def compute_magnon_dispersion(bond_data: Dict,
kpath: Dict,
S: float,
ordering: Dict) -> Tuple[np.ndarray, Dict]:
"""Compute the magnon dispersion in the primitive cell.
.. math:: ω_k(q) = 2S × \\bigl(λ_k(J(q)) − λ_{\\min}(q_0)\\bigr)
where ``λ_k(J(q))`` are the eigenvalues of the Hermitian exchange matrix
at each k-point along the path, sorted ascending. This gives exactly
``n_sub`` branches in the primitive BZ for any ordering wavevector q₀.
The Goldstone mode (ω = 0) is at ``k = q₀`` by construction of the shift,
so no supercell or spectral unfolding is needed.
"""
q0_c = ordering["q0_cart"]
q0_norm = float(np.linalg.norm(q0_c)) # noqa: F841 — kept for parity with original
def lam_min_fn(q_vec):
q = np.array(q_vec, dtype=float).reshape(1, 3)
Jq = build_jq_matrices(bond_data, q)
Jh = 0.5 * (Jq[0] + Jq[0].T.conj())
return float(np.linalg.eigvalsh(Jh).min())
result = minimize(lam_min_fn, q0_c, method="L-BFGS-B",
options={"ftol": 1e-14, "gtol": 1e-10, "maxiter": 500})
lam0_exact = float(result.fun)
log.info(" lambda_min(q0) refined = %.6f meV.S2", lam0_exact)
q_cart = kpath["kpoints_cart"]
dist = kpath["distances"]
distance_breaks = {i for i in range(1, len(dist))
if abs(dist[i] - dist[i - 1]) < 1e-10}
tick_reset = {int(np.argmin(np.abs(dist - tp)))
for tp in kpath["tick_positions"]}
all_resets = distance_breaks | tick_reset
kpath_out = copy.copy(kpath)
kpath_out["reset_indices"] = all_resets
N_q = len(q_cart)
n_sub = bond_data["n_sub"]
Jq_raw = build_jq_matrices(bond_data, q_cart)
Jq_herm = 0.5 * (Jq_raw + Jq_raw.swapaxes(1, 2).conj())
evals = np.zeros((N_q, n_sub))
for iq in range(N_q):
evals[iq] = np.linalg.eigvalsh(Jq_herm[iq])
magnon_bands = 2.0 * S * (evals - lam0_exact)
magnon_bands = np.clip(magnon_bands, 0.0, None)
log.info(" %d bands, range [%.4f, %.4f] meV",
n_sub, magnon_bands.min(), magnon_bands.max())
return magnon_bands, kpath_out