Skip unnecessary QR decompositions in make_unique_kpoint_irreps

Author: bpuchala

python/libcasm/configuration/_ReciprocalSupercell.py

@@ -7,7 +7,7 @@
     get_all_subgroup_orbits,
     get_cyclic_subgroup_orbits,
 )
-from libcasm.irreps import IrrepDecomposition, IrrepInfo, make_invariant_subspace
+from libcasm.irreps import IrrepDecomposition, IrrepInfo
 from libcasm.xtal import (
     UnitCellIndexConverter,
     min_periodic_displacement,
@@ -222,6 +222,14 @@ def __init__(self, supercell: Supercell, tol: float = 1e-5):
         K-points are ordered as generated by :py:attr:`index_converter`.
         """
 
+        self.neg_kpoint_index = []
+        """list[int | None]: For each k-point index i, the index of −k in
+        `self.coordinates` (mod primitive reciprocal lattice), or None if not found.
+
+        For TRIM points where k ≡ −k (mod primitive reciprocal lattice),
+        ``neg_kpoint_index[i] == i``.
+        """
+
         self.orbits = []
         """list[list[int]]: Symmetry orbits of k-point indices under the supercell
         factor group.
@@ -232,8 +240,14 @@ def __init__(self, supercell: Supercell, tol: float = 1e-5):
         """
 
         self.equivalence_map = []
-        """list[list[libcasm.xtal.SymOp]: List of symmetry operations mapping the first
-        k-point in each orbit to the other k-points in the orbit.
+        """list[list[libcasm.configuration.SupercellSymOp]]: List of SupercellSymOp
+        operations mapping the first k-point in each orbit to the other k-points in
+        the orbit.
+
+        Each element ``equivalence_map[i_orbit][j]`` is a factor group
+        SupercellSymOp (``translation_index() == 0``) that maps
+        ``coordinates[:, orbits[i_orbit][0]]`` to
+        ``coordinates[:, orbits[i_orbit][j]]``.
         """
 
         self.little_groups = []
@@ -244,6 +258,16 @@ def __init__(self, supercell: Supercell, tol: float = 1e-5):
         supercell factor group that leave the k-point invariant.
         """
 
+        self.orbit_independent_indices = []
+        """list[list[int]]: For each orbit, the positions j within the orbit whose
+        real-valued plane wave subspace is independent of all earlier members.
+
+        Position 0 (the prototype) is always included. A later position j is excluded
+        when ``coordinates[:, orbits[i_orbit][j]]`` is the negative (mod primitive
+        reciprocal lattice) of an already-independent k-point in the same orbit,
+        because the real-valued plane wave basis for k and −k spans the same subspace.
+        """
+
         # build data
         self._build_kpoints()
         self._build_kpoint_orbits()
@@ -272,6 +296,15 @@ def _build_kpoints(self):
 
         # Compute k-points by applying S_recip to the integer index columns
         self.coordinates = pretty(S_recip @ self.indices)
+
+        # For each k-point, find the index of its negative (mod primitive reciprocal
+        # lattice), or None if not found
+        n_kpts = self.coordinates.shape[1]
+        neg_kpoint_index = []
+        for i in range(n_kpts):
+            idx = self._get_kpoint_index(-self.coordinates[:, i])
+            neg_kpoint_index.append(idx if idx != -1 else None)
+        self.neg_kpoint_index = neg_kpoint_index
         return
 
     def _get_kpoint_index(self, kpt: np.array) -> int:
@@ -293,14 +326,12 @@ def _build_kpoint_orbits(self):
 
         found = [False] * self.coordinates.shape[1]
         kpoint_orbits = []
-        kpoint_equivalence_map = []
         elements = self.supercell.factor_group.elements
         for i in range(self.coordinates.shape[1]):
             if found[i]:
                 continue
             kpt = self.coordinates[:, i]
             orbit = [i]
-            equivalence_map = [elements[0].copy()]
             found[i] = True
             for op in elements:
                 matrix = op.matrix()
@@ -312,46 +343,88 @@ def _build_kpoint_orbits(self):
                     )
                 if not found[kpt_index]:
                     orbit.append(kpt_index)
-                    equivalence_map.append(op.copy())
                     found[kpt_index] = True
-            # Sort orbit and equivalence_map together by kpoint index
-            sorted_pairs = sorted(zip(orbit, equivalence_map), key=lambda p: p[0])
-            orbit, equivalence_map = map(list, zip(*sorted_pairs))
+            orbit.sort()
             kpoint_orbits.append(orbit)
-            kpoint_equivalence_map.append(equivalence_map)
         self.orbits = kpoint_orbits
-        self.equivalence_map = kpoint_equivalence_map
         return
 
     def _build_little_groups(self):
-        """Populate self.little_groups using symmetry operations of the supercell."""
+        """Populate self.little_groups and self.equivalence_map using symmetry
+        operations of the supercell."""
         if self.coordinates is None:
             raise RuntimeError("coordinates must be built before computing orbits")
 
-        little_groups = []
-        for i in range(self.coordinates.shape[1]):
-            little_group = []
-
-            kpt = self.coordinates[:, i]
-
-            it = SupercellSymOp.begin(self.supercell)
-            end = SupercellSymOp.end(self.supercell)
-            while it != end:
-                op = it.copy()
-                matrix = op.to_symop().matrix()
-                kpt_transformed = matrix @ kpt
-                kpt_index = self._get_kpoint_index(kpt_transformed)
-
-                if kpt_index == -1:
+        n_kpts = self.coordinates.shape[1]
+
+        # Build lookup tables from k-point index to orbit position
+        orbit_index_of = {}
+        position_in_orbit = {}
+        for i_orbit, orbit in enumerate(self.orbits):
+            for j, kpt_idx in enumerate(orbit):
+                orbit_index_of[kpt_idx] = i_orbit
+                position_in_orbit[kpt_idx] = j
+
+        prototype_indices = [orbit[0] for orbit in self.orbits]
+
+        little_groups = [[] for _ in range(n_kpts)]
+        equivalence_map = [[None] * len(orbit) for orbit in self.orbits]
+
+        it = SupercellSymOp.begin(self.supercell)
+        end = SupercellSymOp.end(self.supercell)
+        while it != end:
+            op = it.copy()
+            matrix = op.to_symop().matrix()
+            is_factor_group_op = op.translation_index() == 0
+
+            # Compute transformed k-point index for each initial k-point index
+            kpt_transformed_indices = []
+            for kpt_index_init in range(n_kpts):
+                kpt = self.coordinates[:, kpt_index_init]
+                kpt_index_final = self._get_kpoint_index(matrix @ kpt)
+                if kpt_index_final == -1:
                     raise ValueError(
                         "Transformed k-point not found in list of k-points"
                     )
-                if kpt_index == i:
-                    little_group.append(op)
-                it.next()
+                kpt_transformed_indices.append(kpt_index_final)
+                if kpt_index_final == kpt_index_init:
+                    little_groups[kpt_index_init].append(op)
+
+            # Use factor group ops to build equivalence_map (maps prototype k-point
+            # in each orbit to each other k-point in the orbit)
+            if is_factor_group_op:
+                for i_orbit, proto_idx in enumerate(prototype_indices):
+                    dest_idx = kpt_transformed_indices[proto_idx]
+                    j_dest = position_in_orbit.get(dest_idx)
+                    if j_dest is not None and equivalence_map[i_orbit][j_dest] is None:
+                        equivalence_map[i_orbit][j_dest] = op
+
+            it.next()
 
-            little_groups.append(little_group)
         self.little_groups = little_groups
+        self.equivalence_map = equivalence_map
+
+        # For each orbit, determine which positions j contribute an independent
+        # real-valued plane wave subspace. Position j is dependent if k_j is the
+        # negative of an already-independent k-point, because the real-valued cos/sin
+        # basis for k and −k spans the same subspace.
+        orbit_independent_indices = []
+        for orbit in self.orbits:
+            independent = [0]
+            covered = {orbit[0]}
+            neg = self.neg_kpoint_index[orbit[0]]
+            if neg is not None:
+                covered.add(neg)
+            for j in range(1, len(orbit)):
+                if orbit[j] in covered:
+                    continue
+                independent.append(j)
+                covered.add(orbit[j])
+                neg = self.neg_kpoint_index[orbit[j]]
+                if neg is not None:
+                    covered.add(neg)
+            orbit_independent_indices.append(independent)
+        self.orbit_independent_indices = orbit_independent_indices
 
 
 class SupercellDoF:
@@ -1141,14 +1214,10 @@ def make_unique_kpoint_irreps(
     in the supercell, and combines the symmetry-adapted bases from each irrep
     decomposition into a full symmetry-adapted DoF space for the supercell.
 
-    Note that an irrep subspaces is constructed for a single k-point in each orbit
-    and operations from the equivalence map are used to transform the subspace to the
-    subspaces for the other k-points in the orbit. Therefore, the axes of the combined
-    symmetry-adapted basis are not completely symmetrized as is obtained by
-    an irrep decomposition for DoFSpace.
-
-    these subspaces are not
-    fully symmetrized with respect to the symmetry of the full supercell.
+    Note that this method results in symmetry-adapted modes that are restricted to the
+    plane wave subspace for each k-point and do not mix with each other. This may
+    result in axes that are less symmetrized than if the symmetry-adapted modes were
+    constructed using the full supercell symmetry from the start.
 
     Parameters
     ----------
@@ -1177,36 +1246,23 @@ def make_unique_kpoint_irreps(
     # Get the full DoF space for the supercell
     full_dof_space = supercell_dof.dof_space
 
-    # Make matrix representation of supercell factor group on full DoF space
-    factor_group = []
-    it = SupercellSymOp.begin(supercell_kpoints.supercell)
-    end = SupercellSymOp.end(supercell_kpoints.supercell)
-    while it != end:
-        if it.translation_index() == 0:
-            factor_group.append(it.copy())
-        it.next()
-    matrix_rep = make_dof_space_rep(
-        group=factor_group,
-        dof_space=full_dof_space,
-    )
-
     ## Iterate over k-point orbits to build irreps and symmetry-adapted bases ##
 
     # Irreps generated from a single k-point in each orbit
     kpoint_irreps: list[IrrepInfo] = []
 
-    # Bases generated by applying the supercell factor group to expand the symmetry
-    # adapted basis from a single k-point in each orbit to the full orbit. This
-    # will be combined at the end to form the full symmetry-adapted DoF space.
+    # Bases generated by applying the equivalence map to expand the symmetry-adapted
+    # basis from the prototype k-point in each orbit to the other k-points in the
+    # orbit. Combined at the end to form the full symmetry-adapted DoF space.
     orbit_adapted_bases: list[np.ndarray] = []
 
     # Information about each axis in the final symmetry-adapted basis
     axis_irrep_info = []
     next_irrep_char_index = 0
 
     # For each k-point orbit, build irreps for the plane wave basis from the first
-    # k-point in the orbit, then apply the supercell factor group to extend the
-    # symmetry-adapted basis to the space for the full orbit.
+    # k-point in the orbit, then use the equivalence map to extend the
+    # symmetry-adapted basis to the other k-points in the orbit.
     for i_orbit, orbit in enumerate(supercell_kpoints.orbits):
         # Get irreps for first k-point in orbit
         kpoint_index = orbit[0]
@@ -1220,19 +1276,49 @@ def make_unique_kpoint_irreps(
         B = irrep_decomp.symmetry_adapted_subspace
         B = pretty(B)
 
+        # Matrix reps for independent orbit members (skipping −k partners)
+        independent_indices = supercell_kpoints.orbit_independent_indices[i_orbit]
+        independent_ops = [
+            supercell_kpoints.equivalence_map[i_orbit][j] for j in independent_indices
+        ]
+        equivalence_map_reps = make_dof_space_rep(
+            group=independent_ops,
+            dof_space=full_dof_space,
+        )
+
         # For each irrep in the decomposition, extend basis to full orbit and
         # store axis irrep info
         i_axis = 0
         for i_irrep, irrep in enumerate(irrep_decomp.irreps):
 
             d = irrep.irrep_dim
 
-            # Apply supercell factor group to adapt basis to full symmetry
+            # Apply equivalence map ops for independent orbit members to generate
+            # orthogonal subspaces. −k partners are skipped because the real-valued
+            # plane wave basis for k and −k spans the same subspace.
             kpoint_adapted_basis = B[:, i_axis : i_axis + d]
-            orbit_adapted_basis = make_invariant_subspace(
-                matrix_rep=matrix_rep,
-                init_subspace=kpoint_adapted_basis,
-            )
+            # Sanity check: columns of kpoint_adapted_basis should be orthonormal
+            inner = kpoint_adapted_basis.T @ kpoint_adapted_basis
+            if not np.allclose(inner, np.eye(d), atol=1e-8):
+                raise ValueError(
+                    "Error in make_unique_kpoint_irreps: kpoint_adapted_basis "
+                    "columns are not orthonormal."
+                )
+            orbit_cols = [kpoint_adapted_basis]
+            for D in equivalence_map_reps[1:]:
+                transformed = D @ kpoint_adapted_basis
+                # Sanity check: transformed subspace should be orthogonal to prototype
+                overlap = kpoint_adapted_basis.T @ transformed
+                if not np.allclose(overlap, np.zeros_like(overlap), atol=1e-8):
+                    raise ValueError(
+                        "Error in make_unique_kpoint_irreps: equivalence map "
+                        "transformed subspace is not orthogonal to prototype subspace."
+                    )
+                orbit_cols.append(transformed)
+            orbit_adapted_basis = np.hstack(orbit_cols)
+            # q, r = np.linalg.qr(orbit_adapted_basis)
+            # rank = np.linalg.matrix_rank(r)
+            # orbit_adapted_basis = pretty(q[:, :rank])
 
             orbit_adapted_bases.append(orbit_adapted_basis)
 
@@ -1258,11 +1344,20 @@ def make_unique_kpoint_irreps(
             kpoint_irreps.append(irrep)
             i_axis += d
 
-    # Combine all symmetry-adapted bases from each irrep decomposition
-    symmetry_adapted_basis = np.hstack(orbit_adapted_bases)
-    q, r = np.linalg.qr(symmetry_adapted_basis)
-    rank = np.linalg.matrix_rank(r)
-    symmetry_adapted_basis = pretty(q[:, :rank])
+    # Combine all symmetry-adapted bases from each irrep decomposition. Each
+    # orbit-adapted basis should be orthogonal to all previously accumulated bases.
+    symmetry_adapted_basis = orbit_adapted_bases[0]
+    for basis in orbit_adapted_bases[1:]:
+        overlap = symmetry_adapted_basis.T @ basis
+        if not np.allclose(overlap, np.zeros_like(overlap), atol=1e-8):
+            raise ValueError(
+                "Error in make_unique_kpoint_irreps: orbit-adapted basis is not "
+                "orthogonal to previously accumulated symmetry-adapted basis."
+            )
+        symmetry_adapted_basis = np.hstack([symmetry_adapted_basis, basis])
+    # q, r = np.linalg.qr(symmetry_adapted_basis)
+    # rank = np.linalg.matrix_rank(r)
+    # symmetry_adapted_basis = pretty(q[:, :rank])
 
     dof_space = DoFSpace(
         dof_key=supercell_dof.dof_key,

python/tests/configuration/test_SupercellKpoints.py

@@ -146,7 +146,7 @@ def test_make_unique_kpoint_irreps(self, sc_kpts, sc_dof):
 
     def test_axis_irrep_info(self, sc_kpts, sc_dof):
         """Each axis_irrep_info entry has valid orbit_index and kpoint_irreps_index."""
-        _, kpoint_irreps, axis_irrep_info = casmconfig.make_unique_kpoint_irreps(
+        kpoint_irreps, _, axis_irrep_info = casmconfig.make_unique_kpoint_irreps(
             supercell_kpoints=sc_kpts,
             supercell_dof=sc_dof,
         )