"""q-point grid generation for the BZ scan."""
from __future__ import annotations
import logging
from typing import Tuple
import numpy as np
from pymatgen.core import Structure
log = logging.getLogger(__name__)
[docs]
def build_qmesh(structure: Structure,
na: int, nb: int, nc: int) -> Tuple[np.ndarray, np.ndarray]:
"""Γ-centred uniform q-mesh over the full primitive BZ.
Returns
-------
q_cart, q_frac : ndarray
Cartesian (Å⁻¹) and fractional reciprocal-lattice coordinates,
each of shape (na·nb·nc, 3).
"""
rec = structure.lattice.reciprocal_lattice
def _axis(n):
x = np.arange(n) / n
x[x >= 0.5] -= 1.0
return x
QA, QB, QC = np.meshgrid(_axis(na), _axis(nb), _axis(nc), indexing='ij')
q_frac = np.stack([QA.ravel(), QB.ravel(), QC.ravel()], axis=1)
return q_frac @ rec.matrix, q_frac
[docs]
def auto_qmesh(structure: Structure,
target_points: int = 125000) -> Tuple[int, int, int]:
"""Anisotropy-adapted q-mesh targeting ~``target_points`` q-points."""
lattice = structure.lattice
rec = lattice.reciprocal_lattice
real_lengths = np.array([lattice.a, lattice.b, lattice.c])
rec_lengths = np.array([rec.a, rec.b, rec.c ])
axis_names = ["a", "b", "c"]
shortest_idx = int(np.argmin(real_lengths))
log.info("Real-space cell: a=%.4f b=%.4f c=%.4f A", *real_lengths)
log.info("Reciprocal cell: a*=%.4f b*=%.4f c*=%.4f A-1", *rec_lengths)
log.info("Shortest real-space axis: %s = %.4f A",
axis_names[shortest_idx], real_lengths[shortest_idx])
rec_max = rec_lengths.max()
weights = rec_lengths / rec_max
scale = (target_points / np.prod(weights)) ** (1.0 / 3.0)
ns = np.maximum(1, np.round(weights * scale)).astype(int)
log.info("Auto q-mesh -> %d x %d x %d = %d points (target %d)",
ns[0], ns[1], ns[2], int(np.prod(ns)), target_points)
return int(ns[0]), int(ns[1]), int(ns[2])