Add DFTB-SCC support to DeePTB

docs/adr/0003-kpoint-coordinate-convention.md

@@ -0,0 +1,45 @@
+# ADR 0003: Use First-Brillouin-Zone Coordinates For K-Point Meshes
+
+## Status
+
+Accepted
+
+## Context
+
+K-point mesh generation, path construction, and symmetry reduction are used by band, DOS, SCC, transport, and export workflows. Historically, related helpers lived in `dptb.utils.make_kpoints` and `dptb.utils.ksampling`, with overlapping responsibilities and different assumptions about how fractional k-points should be represented.
+
+Gamma-centered meshes previously exposed points in a `[0, 1)` fractional-coordinate convention in some utility paths. Symmetry reduction, especially time-reversal reduction, needs to compare `k` and `-k` while handling Brillouin-zone boundary points such as `0.5` and `-0.5`. Keeping multiple conventions in active implementation code makes boundary handling easy to get wrong.
+
+## Decision
+
+Core k-point implementation lives in `dptb.kpoints`, split by responsibility:
+
+- `dptb.kpoints.mesh` for mesh generation;
+- `dptb.kpoints.path` for high-symmetry paths;
+- `dptb.kpoints.reduction` for symmetry and time-reversal reduction;
+- `dptb.kpoints.sampling` for higher-level sampling.
+
+Mesh k-points should be wrapped to the first-Brillouin-zone fractional-coordinate convention `[-0.5, 0.5)`. For example, a Gamma-centered `[4, 1, 1]` mesh is represented as:
+
+```text
+0.00, 0.25, -0.50, -0.25
+```
+
+rather than:
+
+```text
+0.00, 0.25, 0.50, 0.75
+```
+
+These are equivalent modulo reciprocal lattice vectors (`0.75 == -0.25`, `0.50 == -0.50`), but `[-0.5, 0.5)` is the canonical representation inside DeePTB k-point utilities.
+
+`dptb.utils.make_kpoints` and `dptb.utils.ksampling` remain compatibility wrappers. New code should import from `dptb.kpoints` directly.
+
+## Consequences
+
+- Symmetry reduction has one canonical coordinate convention for comparing equivalent k-points.
+- Time-reversal reduction can treat Brillouin-zone boundary points consistently.
+- Physical calculations that consume k-points through periodic Hamiltonian evaluation should be unchanged by this representation choice.
+- User code that assumes Gamma-centered mesh coordinates are non-negative or lie in `[0, 1)` may observe different raw k-point arrays and should wrap or compare k-points modulo reciprocal lattice vectors.
+- Regression tests should cover both `dptb.kpoints` core functions and the compatibility wrappers under `dptb.utils`.
+- AI-assisted refactors should not reintroduce duplicated mesh or reduction logic in downstream modules.

docs/adr/index.md

@@ -7,5 +7,5 @@ Architecture Decision Records explain decisions that future reviews should not r
 
 0001-reduced-matrix-elements
 0002-orbitalmapper-as-index-authority
+0003-kpoint-coordinate-convention
 ```
-

dptb/kpoints/__init__.py

@@ -0,0 +1,57 @@
+"""K-point mesh generation, path construction, and symmetry reduction."""
+
+from dptb.kpoints.geometry import (
+    calculate_reciprocal_vectors,
+    get_symm_ops,
+    is_integer,
+    rot_revlatt_2D,
+)
+from dptb.kpoints.mesh import (
+    build_kmeshgrid,
+    gamma_center,
+    kgrid_spacing,
+    kmesh_fs,
+    kmesh_sampling,
+    kmesh_sampling_negf,
+    monkhorst_pack,
+    monkhorst_pack_sampling,
+    mp,
+    time_symmetry_reduce,
+)
+from dptb.kpoints.path import (
+    abacus_kpath,
+    ase_kpath,
+    vasp_kpath,
+)
+from dptb.kpoints.reduction import (
+    reduce,
+    reduce_rotation,
+    reduce_time_inversion,
+)
+from dptb.kpoints.sampling import (
+    sample,
+)
+
+__all__ = [
+    "abacus_kpath",
+    "ase_kpath",
+    "build_kmeshgrid",
+    "calculate_reciprocal_vectors",
+    "gamma_center",
+    "get_symm_ops",
+    "is_integer",
+    "kgrid_spacing",
+    "kmesh_fs",
+    "kmesh_sampling",
+    "kmesh_sampling_negf",
+    "monkhorst_pack",
+    "monkhorst_pack_sampling",
+    "mp",
+    "reduce",
+    "reduce_rotation",
+    "reduce_time_inversion",
+    "rot_revlatt_2D",
+    "sample",
+    "time_symmetry_reduce",
+    "vasp_kpath",
+]

dptb/kpoints/geometry.py

@@ -0,0 +1,74 @@
+import logging
+from typing import List, Tuple
+
+import numpy as np
+from ase.atoms import Atoms
+
+log = logging.getLogger(__name__)
+
+
+def is_integer(num):
+    return abs(round(num) - num) < 1e-10
+
+
+def calculate_reciprocal_vectors(
+    a: np.ndarray,
+    b: np.ndarray,
+    c: np.ndarray,
+) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Calculate reciprocal lattice vectors."""
+    vol = np.dot(a, np.cross(b, c))
+    if np.abs(vol) <= 1e-10:
+        raise ValueError("Lattice vectors are singular and cannot define reciprocal vectors.")
+    rmat = np.vstack([x.reshape((1, 3)) for x in [a, b, c]])
+    gmat = np.linalg.solve(rmat, np.eye(3)).T * 2 * np.pi
+    return gmat[0], gmat[1], gmat[2]
+
+
+def get_symm_ops(atoms: Atoms) -> List[np.ndarray]:
+    """Get rotational symmetry operations for an atomic structure."""
+    from spglib import get_symmetry as _spglib_get_symmetry
+
+    latt = atoms.get_cell()
+    coords = atoms.get_positions()
+    type_map = atoms.get_atomic_numbers()
+    symmetry = _spglib_get_symmetry((latt, coords, type_map))
+    if symmetry is None or "rotations" not in symmetry:
+        raise ValueError("spglib failed to determine rotational symmetry operations.")
+    rotations = symmetry["rotations"]
+    if len(rotations) == 0:
+        raise ValueError("spglib returned no rotational symmetry operations.")
+    return [op for op in rotations]
+
+
+def rot_revlatt_2D(rev_latt, index=None):  # 0, x; 1,y, 2,z
+    """Transform reciprocal lattice vectors to a 2D-oriented coordinate system."""
+    if index is None:
+        index = (0, 1)
+    if len(index) != 2 or len(set(index)) != 2 or not all(isinstance(i, int) and i in {0, 1, 2} for i in index):
+        raise ValueError("index must contain two distinct integer axes from {0, 1, 2}.")
+    rev_latt = np.asarray(rev_latt)
+    if rev_latt.shape != (3, 3):
+        log.error("Error! rev_latt must be a 3x3 array!")
+        raise ValueError("rev_latt must be a 3x3 array.")
+
+    index_left = [0, 1, 2]
+    for i in index:
+        index_left.remove(i)
+
+    vec1 = np.array(rev_latt[index[0]]).reshape(-1)
+    vec2 = np.array(rev_latt[index[1]]).reshape(-1)
+
+    avec1 = vec1 / np.linalg.norm(vec1)
+    avec3 = np.cross(avec1, vec2) / np.linalg.norm(np.cross(avec1, vec2))
+    avec2 = np.cross(avec3, avec1)
+    if np.dot(np.cross(avec1, avec2), avec3) < 0:
+        avec3 = -avec3
+
+    newcorr = np.zeros((3, 3))
+    newcorr[index[0]] = avec1
+    newcorr[index[1]] = avec2
+    newcorr[index_left[0]] = avec3
+
+    rev_latt_new = rev_latt @ np.linalg.inv(newcorr)
+    return rev_latt_new, newcorr

dptb/kpoints/mesh.py

@@ -0,0 +1,159 @@
+import logging
+from itertools import product as itprod
+from typing import List, Optional, Union
+
+import ase
+import numpy as np
+
+from dptb.kpoints.geometry import calculate_reciprocal_vectors, is_integer
+
+log = logging.getLogger(__name__)
+
+
+def _validate_meshgrid(meshgrid) -> List[int]:
+    if not isinstance(meshgrid, (list, tuple, np.ndarray)) or len(meshgrid) != 3:
+        raise ValueError("meshgrid must be a sequence of 3 positive integers.")
+    meshgrid = np.asarray(meshgrid)
+    if not np.all(np.isfinite(meshgrid)) or not np.all(np.equal(meshgrid, meshgrid.astype(int))):
+        raise ValueError("meshgrid must contain 3 positive integers.")
+    meshgrid = meshgrid.astype(int)
+    if not (meshgrid > 0).all():
+        raise ValueError("meshgrid must contain 3 positive integers.")
+    return meshgrid.tolist()
+
+
+def mp(
+    nk1: int,
+    nk2: int,
+    nk3: int,
+    b1: Optional[np.ndarray] = None,
+    b2: Optional[np.ndarray] = None,
+    b3: Optional[np.ndarray] = None,
+    gamma_centered: bool = True,
+    direct: bool = True,
+):
+    """Monkhorst-Pack sampling, including Gamma-centered and shifted grids."""
+    nk1, nk2, nk3 = _validate_meshgrid([nk1, nk2, nk3])
+
+    if gamma_centered:
+        def k_1d(n):
+            k = np.arange(n) / n
+            k[k >= 0.5] -= 1
+            return k
+    else:
+        k_1d = lambda n: (np.arange(n) + 0.5) / n - 0.5
+
+    k_taud = np.array(list(itprod(*[k_1d(n) for n in [nk1, nk2, nk3]])))
+    k_taud = k_taud.reshape(-1, 3)
+    if len(k_taud) != nk1 * nk2 * nk3:
+        raise RuntimeError("Generated k-point mesh size does not match the requested grid.")
+
+    if direct:
+        return k_taud
+
+    if not all(x is not None for x in [b1, b2, b3]):
+        raise ValueError("Reciprocal vectors b1, b2, and b3 are required when direct=False.")
+    return np.tensordot(k_taud, np.array([b1, b2, b3]), axes=(1, 0))
+
+
+def build_kmeshgrid(b1, b2, b3, kspac: Union[int, float, List[float]]) -> List[int]:
+    """Build a k-mesh grid from reciprocal vectors and target k-spacing."""
+    norms = [np.linalg.norm(x) for x in [b1, b2, b3]]
+    kspac = [kspac] * 3 if isinstance(kspac, (int, float)) else kspac
+    if not isinstance(kspac, (list, tuple)) or len(norms) != len(kspac):
+        raise ValueError(f"kspacing should be a list of 3 positive floats: {kspac}")
+    if not all(np.isfinite(x) and x > 0 for x in kspac):
+        raise ValueError("kspacing values must be positive and finite.")
+    norms = [int(norm / kspac) for norm, kspac in zip(norms, kspac)]
+    return list(map(lambda x: max(1, x + 1), norms))
+
+
+def monkhorst_pack_sampling(
+    nk1: int = 1,
+    nk2: int = 1,
+    nk3: int = 1,
+    cell: Optional[np.ndarray] = None,
+    kspac: Optional[Union[int, float, List[float]]] = None,
+    gamma_centered: bool = True,
+    direct: bool = True,
+) -> np.ndarray:
+    """Sample k-points by explicit mesh or by target k-spacing."""
+    b1, b2, b3 = None, None, None
+
+    if kspac is not None:
+        cell = np.asarray(cell)
+        if cell.shape not in [(3, 3), (9,)]:
+            raise ValueError("cell must have shape (3, 3) or (9,) when kspac is provided.")
+        cell = cell.reshape((3, 3))
+        b1, b2, b3 = calculate_reciprocal_vectors(cell[0], cell[1], cell[2])
+        nk1, nk2, nk3 = build_kmeshgrid(b1, b2, b3, kspac)
+
+    return mp(nk1, nk2, nk3, b1, b2, b3, gamma_centered=gamma_centered, direct=direct)
+
+
+def monkhorst_pack(meshgrid=[1, 1, 1]):
+    """Generate shifted Monkhorst-Pack k-points in [-0.5, 0.5)."""
+    meshgrid = _validate_meshgrid(meshgrid)
+    return mp(*meshgrid, gamma_centered=False)
+
+
+def gamma_center(meshgrid=[1, 1, 1]):
+    """Generate Gamma-centered k-points in [-0.5, 0.5)."""
+    meshgrid = _validate_meshgrid(meshgrid)
+    return mp(*meshgrid, gamma_centered=True)
+
+
+def kmesh_sampling(meshgrid=[1, 1, 1], is_gamma_center=True):
+    """Generate a Gamma-centered or shifted Monkhorst-Pack mesh."""
+    return gamma_center(meshgrid) if is_gamma_center else monkhorst_pack(meshgrid)
+
+
+def time_symmetry_reduce(meshgrid=[1, 1, 1], is_gamma_center=True):
+    """Generate a mesh and reduce it by time-reversal symmetry."""
+    from dptb.kpoints.reduction import reduce_time_inversion
+
+    kpoints = kmesh_sampling(meshgrid, is_gamma_center=is_gamma_center)
+    reduced, weights = reduce_time_inversion(kpoints)
+    weights = weights / len(kpoints)
+    if abs(weights.sum() - 1.0) >= 1e-5:
+        raise RuntimeError("The sum of k-point weights is not 1.0.")
+    return reduced, weights
+
+
+def kmesh_sampling_negf(meshgrid=[1, 1, 1], is_gamma_center=True, is_time_reversal=True):
+    """Generate k-points for NEGF calculations."""
+    if is_time_reversal:
+        return time_symmetry_reduce(meshgrid, is_gamma_center=is_gamma_center)
+
+    kpoints = kmesh_sampling(meshgrid, is_gamma_center=is_gamma_center)
+    wk = np.ones(len(kpoints)) / len(kpoints)
+    return kpoints, wk
+
+
+def kmesh_fs(meshgrid=[1, 1, 1]):
+    """Generate endpoint-inclusive k-points for Fermi-surface calculations."""
+    nx, ny, nz = meshgrid
+    lx, ly, lz = np.linspace(0, 1, nx), np.linspace(0, 1, ny), np.linspace(0, 1, nz)
+    xx, yy, zz = np.meshgrid(lx, ly, lz, indexing="ij")
+    kgrids = np.array([xx.reshape(-1), yy.reshape(-1), zz.reshape(-1)]).T
+    return (lx, ly, lz), kgrids
+
+
+def kgrid_spacing(structase, kspacing: float, sampling="MP"):
+    """Generate k-points from a target k-spacing."""
+    if not isinstance(structase, ase.Atoms):
+        raise TypeError("structase must be an ase.Atoms instance.")
+    if not np.isfinite(kspacing) or kspacing <= 0:
+        raise ValueError("kspacing must be positive and finite.")
+    rev_latt = 2 * np.pi * np.linalg.inv(np.array(structase.cell)).T
+    meshgrid = np.maximum(1, np.floor(np.linalg.norm(rev_latt, axis=1) / kspacing).astype(int) + 1)
+
+    if sampling == "MP":
+        kpoints = monkhorst_pack(meshgrid)
+    elif sampling == "Gamma":
+        kpoints = gamma_center(meshgrid)
+    else:
+        log.error("Error! sampling must be either 'MP' or 'Gamma'!")
+        raise ValueError("sampling must be either 'MP' or 'Gamma'.")
+
+    return kpoints