Refactor acoustic mode masking in Q-space modules

Replace hardcoded acoustic_eps threshold with a unified valid_modes_q mask
that combines translation-based detection at Gamma with optional frequency-
based masking when ignore_small_w is True.

Changes:
- QSpaceLanczos.py: Build valid_modes_q mask using translations at Gamma
  plus small-frequency masking at all q-points when ignore_small_w=True.
  Pass mask to Julia extension.

- QSpaceHessian.py: Replace all acoustic_eps usages with valid_modes_q
  mask lookups for consistent mode exclusion.

- tdscha_qspace.jl: Accept valid_modes_q mask from Python and use it
  instead of hardcoded 1e-6 threshold for zeroing f_Y/f_psi.

All 13 q-space tests pass.

Author: mesonepigreco

Modules/QSpaceHessian.py

@@ -215,7 +215,6 @@ def _find_degenerate_blocks(self, iq):
             E.g. [[3], [4, 5, 6], [7, 8, 9]] for a singlet + two triplets.
         """
         w_qp = self.w_q[:, iq]
-        eps = self.qlanc.acoustic_eps
 
         deg_lists = CC.symmetries.get_degeneracies(w_qp)
 
@@ -224,12 +223,12 @@ def _find_degenerate_blocks(self, iq):
         for nu in range(self.n_bands):
             if nu in seen:
                 continue
-            if w_qp[nu] < eps:
+            if not self.qlanc.valid_modes_q[nu, iq]:
                 seen.add(nu)
                 continue
             block = sorted(deg_lists[nu].tolist())
-            # Filter out any acoustic modes that crept in
-            block = [b for b in block if w_qp[b] >= eps]
+            # Filter out any masked modes that crept in
+            block = [b for b in block if self.qlanc.valid_modes_q[b, iq]]
             for b in block:
                 seen.add(b)
             if len(block) > 0:
@@ -281,38 +280,35 @@ def _static_mask(self):
         """Build mask for static psi inner product.
 
         Same structure as FT mask but only one W sector (not a' + b').
-        Acoustic modes (w < eps) are masked out (set to 0) to keep
+        Masked modes (acoustic/small-w) are set to 0 to keep
         the null space purely in the W sector and avoid inconsistency.
         """
         psi_size = self._get_static_psi_size()
         mask = np.ones(psi_size, dtype=np.float64)
         nb = self.n_bands
-        eps = self.qlanc.acoustic_eps
 
-        # R sector: mask acoustic modes at the perturbation q-point
-        w_qp = self.w_q[:, self.qlanc.iq_pert]
-        mask[:nb] = np.where(w_qp < eps, 0.0, 1.0)
+        # R sector: mask modes at the perturbation q-point
+        valid_pert = self.qlanc.valid_modes_q[:, self.qlanc.iq_pert]
+        mask[:nb] = np.where(valid_pert, 1.0, 0.0)
 
         for pair_idx, (iq1, iq2) in enumerate(self.qlanc.unique_pairs):
             offset = self.qlanc.get_block_offset(pair_idx, 'a')
-            w1 = self.w_q[:, iq1]
-            w2 = self.w_q[:, iq2]
-            ac1 = w1 < eps
-            ac2 = w2 < eps
-            # acoustic_any[i,j] = True if either mode i or j is acoustic
-            acoustic_any = ac1[:, None] | ac2[None, :]
+            valid1 = self.qlanc.valid_modes_q[:, iq1]
+            valid2 = self.qlanc.valid_modes_q[:, iq2]
+            # invalid_any[i,j] = True if either mode i or j is masked
+            invalid_any = ~valid1[:, None] | ~valid2[None, :]
 
             if iq1 < iq2:
-                # Full block: 0 if acoustic, 2 if not
-                block_mask = np.where(acoustic_any, 0.0, 2.0)
+                # Full block: 0 if masked, 2 if not
+                block_mask = np.where(invalid_any, 0.0, 2.0)
                 mask[offset:offset + nb * nb] = block_mask.ravel()
             else:
                 # Upper triangle storage: pack diag (weight=1) + off-diag (weight=2)
                 size = self.qlanc.get_block_size(pair_idx)
-                # Build full matrix of weights: diag=1, off-diag=2, acoustic=0
-                full_weights = np.where(acoustic_any, 0.0, 2.0)
+                # Build full matrix of weights: diag=1, off-diag=2, masked=0
+                full_weights = np.where(invalid_any, 0.0, 2.0)
                 np.fill_diagonal(full_weights, np.where(
-                    ac1 | ac2, 0.0, 1.0))
+                    valid1 & valid2, 1.0, 0.0))
                 # Pack upper triangle in row-major order
                 tri = np.zeros(size, dtype=np.float64)
                 idx = 0
@@ -350,16 +346,16 @@ def _compute_lambda_q(self, iq1, iq2):
         n2 = np.zeros_like(w2)
         if T > __EPSILON__:
             beta = __RyToK__ / T
-            valid1 = w1 > self.qlanc.acoustic_eps
-            valid2 = w2 > self.qlanc.acoustic_eps
+            valid1 = self.qlanc.valid_modes_q[:, iq1]
+            valid2 = self.qlanc.valid_modes_q[:, iq2]
             n1[valid1] = 1.0 / (np.exp(w1[valid1] * beta) - 1.0)
             n2[valid2] = 1.0 / (np.exp(w2[valid2] * beta) - 1.0)
 
         n1_mat = np.tile(n1, (self.n_bands, 1)).T
         n2_mat = np.tile(n2, (self.n_bands, 1))
 
-        valid_mask = np.outer(w1 > self.qlanc.acoustic_eps,
-                              w2 > self.qlanc.acoustic_eps)
+        valid_mask = np.outer(self.qlanc.valid_modes_q[:, iq1],
+                              self.qlanc.valid_modes_q[:, iq2])
 
         # (n1 - n2) / (w1 - w2), regularized for w1 ≈ w2
         diff_n = np.zeros_like(w1_mat)
@@ -395,18 +391,18 @@ def _compute_lambda_q(self, iq1, iq2):
     def _compute_Y_w(self, iq):
         """Compute Y_w = 2*w/(2*n+1) for all bands at q-point iq.
 
-        Acoustic modes (w < eps) get Y_w = 0.
+        Masked modes (acoustic/small-w) get Y_w = 0.
         """
         w = self.w_q[:, iq]
         T = self.qlanc.T
 
         n = np.zeros_like(w)
         if T > __EPSILON__:
-            valid = w > self.qlanc.acoustic_eps
+            valid = self.qlanc.valid_modes_q[:, iq]
             n[valid] = 1.0 / (np.exp(w[valid] * __RyToK__ / T) - 1.0)
 
         Y_w = np.zeros_like(w)
-        valid = w > self.qlanc.acoustic_eps
+        valid = self.qlanc.valid_modes_q[:, iq]
         Y_w[valid] = 2.0 * w[valid] / (2.0 * n[valid] + 1.0)
         return Y_w
 
@@ -421,13 +417,13 @@ def _precompute_static_quantities(self):
           - _cached_yw_outer[pair_idx]: -2*outer(Y_w1,Y_w2) for anharmonic
         """
         nb = self.n_bands
-        eps = self.qlanc.acoustic_eps
 
         # R sector
         w_qp = self.w_q[:, self.qlanc.iq_pert]
         self._cached_w_qp_sq = w_qp ** 2
+        valid_pert = self.qlanc.valid_modes_q[:, self.qlanc.iq_pert]
         self._cached_inv_w_qp_sq = np.where(
-            w_qp > eps, 1.0 / np.where(w_qp > eps, w_qp, 1.0) ** 2, 0.0)
+            valid_pert, 1.0 / np.where(valid_pert, w_qp, 1.0) ** 2, 0.0)
 
         # W sector: per-pair Lambda and Y_w caches
         n_pairs = len(self.qlanc.unique_pairs)
@@ -643,7 +639,7 @@ def apply_M_tilde(x_tilde):
         # 5. Identify non-acoustic bands and degenerate blocks
         w_qp = self.w_q[:, iq]
         non_acoustic = [nu for nu in range(nb)
-                        if w_qp[nu] > self.qlanc.acoustic_eps]
+                        if self.qlanc.valid_modes_q[nu, iq]]
 
         # Build solve schedule: list of (band_to_solve, block_members)
         if use_mode_symmetry:

Modules/QSpaceLanczos.py

@@ -142,7 +142,7 @@ def __init__(self, ensemble, **kwargs):
         # -- Add the q-space attributes --
         qspace_attrs = [
             'q_points', 'n_q', 'n_bands', 'w_q', 'pols_q',
-            'acoustic_eps', 'valid_modes_q', 'X_q', 'Y_q',
+            'valid_modes_q', 'X_q', 'Y_q',
             'iq_pert', 'q_pair_map', 'unique_pairs',
             '_psi_size', '_block_offsets_a', '_block_offsets_b', '_block_sizes',
             '_qspace_sym_data', '_qspace_sym_q_map', 'n_syms_qspace',
@@ -161,17 +161,19 @@ def __init__(self, ensemble, **kwargs):
         # The masses needs to be restricted to the primitive cell only
         self.m = self.m[:self.n_bands]
 
-        # Small frequency threshold for acoustic mode masking
-        self.acoustic_eps = 1e-6
-
-        # Build translation-based valid mask for acoustic modes
-        # At Gamma (iq=0), identify translations from polarization vectors
-        # At q != 0, all modes are valid (acoustic modes only have zero freq at Gamma)
+        # Build valid_modes_q mask for acoustic mode exclusion
+        # Always apply translation-based mask at Gamma (iq=0)
         masses_uc = self.dyn.structure.get_masses_array()
         self.valid_modes_q = np.ones((self.n_bands, self.n_q), dtype=bool)
         trans_mask = CC.Methods.get_translations(np.real(self.pols_q[:, :, 0]), masses_uc)
         self.valid_modes_q[:, 0] = ~trans_mask
 
+        # If ignore_small_w is True, also mask small frequencies at ALL q-points
+        if ensemble.ignore_small_w:
+            for iq in range(self.n_q):
+                small_freq_mask = np.abs(self.w_q[:, iq]) < CC.Phonons.__EPSILON_W__
+                self.valid_modes_q[:, iq] &= ~small_freq_mask
+
         # == 2. Bloch transform ensemble data ==
         self._bloch_transform_ensemble()
 
@@ -780,7 +782,8 @@ def get_combined(start_end):
                 np.int64(self.iq_pert + 1),
                 self.q_pair_map + 1,  # 1-indexed
                 unique_pairs_arr,
-                int(start_end[0]), int(start_end[1])
+                int(start_end[0]), int(start_end[1]),
+                self.valid_modes_q  # Pass mask to Julia
             )
 
         combined = Parallel.GoParallel(get_combined, indices, "+")

Modules/tdscha_qspace.jl

@@ -649,7 +649,8 @@ function get_perturb_averages_qspace(
     q_pair_map::Vector{Int32},
     unique_pairs::Matrix{Int32},
     start_index::Int64,
-    end_index::Int64
+    end_index::Int64,
+    valid_modes_q::Matrix{Bool}  # Mask from Python: false for acoustic/small-w modes
 )
     n_q = size(X_q, 1)
     n_bands = size(X_q, 3)
@@ -662,19 +663,18 @@ function get_perturb_averages_qspace(
     end
 
     # Precompute occupation numbers and scaling factors
-    # Acoustic modes (w < threshold) get f_Y=0, f_psi=0 to avoid NaN/Inf
+    # Masked modes (valid_modes_q == false) get f_Y=0, f_psi=0 to avoid NaN/Inf
     f_Y = zeros(Float64, n_bands, n_q)
     f_psi = zeros(Float64, n_bands, n_q)
-    acoustic_eps = 1e-6
 
     for iq in 1:n_q
         for nu in 1:n_bands
-            w = w_q[nu, iq]
-            if w < acoustic_eps
+            if !valid_modes_q[nu, iq]  # Masked mode -> zero out
                 f_Y[nu, iq] = 0.0
                 f_psi[nu, iq] = 0.0
                 continue
             end
+            w = w_q[nu, iq]
             if temperature > 0
                 nw = 1.0 / (exp(w * RY_TO_K_Q / temperature) - 1.0)
             else