QSpaceLanczos Module

tdscha.QSpaceLanczos

Q-Space Lanczos Module

This module implements the Lanczos algorithm in q-space (Bloch basis) to exploit momentum conservation and block structure from Bloch's theorem. This gives a speedup of ~N_cell over the real-space implementation.

Key differences from the real-space DynamicalLanczos.Lanczos: - Psi vector is complex128 (Hermitian Lanczos with sesquilinear inner product) - Two-phonon sector uses (q1, q2) pairs constrained by q1+q2 = q_pert + G - Symmetries are point-group only (translations handled by Fourier transform) - Requires Julia extension (tdscha_qspace.jl)

References: Implementation plan: implementation_plan.md Parent class: DynamicalLanczos.py

QSpaceLanczos(ensemble, lo_to_split=None, **kwargs)

Bases: Lanczos

Q-space Lanczos for spectral calculations exploiting Bloch's theorem.

This class works in the q-space mode basis to exploit momentum conservation, reducing the psi vector size by ~N_cell and the anharmonic computation by ~N_cell.

Only Wigner formalism is supported. Requires Julia extension.

Initialize the Q-Space Lanczos.

Parameters:

Name	Type	Description	Default
ensemble	Ensemble	The SSCHA ensemble.	required
lo_to_split	string, ndarray, or None	LO-TO splitting mode. If None (default), LO-TO splitting correction is neglected. If "random", a random direction is used. If an ndarray, it specifies the q-direction for the LO-TO splitting correction.	None
**kwargs		Additional keyword arguments passed to the parent Lanczos class.	{}

Source code in tdscha/QSpaceLanczos.py

def __init__(self, ensemble, lo_to_split=None, **kwargs):
    """Initialize the Q-Space Lanczos.

    Parameters
    ----------
    ensemble : sscha.Ensemble.Ensemble
        The SSCHA ensemble.
    lo_to_split : string, ndarray, or None
        LO-TO splitting mode. If None (default), LO-TO splitting correction
        is neglected. If "random", a random direction is used. If an ndarray,
        it specifies the q-direction for the LO-TO splitting correction.
    **kwargs
        Additional keyword arguments passed to the parent Lanczos class.
    """
    # Force Wigner and Julia mode
    kwargs['mode'] = DL.MODE_FAST_JULIA
    kwargs['use_wigner'] = True

    self.ensemble = ensemble
    super().__init__(ensemble, unwrap_symmetries=False, lo_to_split=lo_to_split, **kwargs)

    if not __JULIA_EXT__:
        raise ImportError(
            "QSpaceLanczos requires Julia. Install with: pip install julia"
        )
    self.use_wigner = True

    # -- Add the q-space attributes --
    qspace_attrs = [
        'q_points', 'n_q', 'n_bands', 'w_q', 'pols_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',
        # Distributed mode attributes
        '_distributed', '_N_global', '_N_eff_global', '_N_local',
    ]
    self.__total_attributes__.extend(qspace_attrs)

    # If ensemble is None, perform a bare initialization like the parent
    if ensemble is None:
        return

    # == 1. Get q-space eigenmodes ==
    ws_sc, pols_sc, w_q, pols_q = self.dyn.DiagonalizeSupercell(
        return_qmodes=True, lo_to_split=lo_to_split)

    self.q_points = np.array(self.dyn.q_tot)  # (n_q, 3)
    self.n_q = len(self.q_points)
    self.n_bands = 3 * self.uci_structure.N_atoms  # uniform band count
    self.w_q = w_q        # (n_bands, n_q) from DiagonalizeSupercell
    self.pols_q = pols_q  # (3*n_at, n_bands, n_q) complex eigenvectors

    # The masses needs to be restricted to the primitive cell only
    self.m = self.m[:self.n_bands]

    # 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()

    # Q-space specific state (set by build_q_pair_map)
    self.iq_pert = None
    self.q_pair_map = None
    self.unique_pairs = None
    self._psi_size = None
    self._block_offsets_a = None
    self._block_offsets_b = None
    self._block_sizes = None

    # Symmetry data for q-space
    self._qspace_sym_data = None
    self._qspace_sym_q_map = None
    self.n_syms_qspace = 0

    # Distributed mode attributes (for MPI configuration distribution)
    self._distributed = False
    self._N_global = self.N
    self._N_eff_global = self.N_eff
    self._N_local = self.N

build_q_pair_map(iq_pert)

Find all (iq1, iq2) pairs satisfying q1 + q2 = q_pert + G.

Parameters:

Name	Type	Description	Default
iq_pert	int	Index of the perturbation q-point.	required

Source code in tdscha/QSpaceLanczos.py

def build_q_pair_map(self, iq_pert):
    """Find all (iq1, iq2) pairs satisfying q1 + q2 = q_pert + G.

    Parameters
    ----------
    iq_pert : int
        Index of the perturbation q-point.
    """
    bg = self.uci_structure.get_reciprocal_vectors() / (2 * np.pi)
    q_pert = self.q_points[iq_pert]

    self.iq_pert = iq_pert
    self.q_pair_map = np.zeros(self.n_q, dtype=np.int32)

    for iq1 in range(self.n_q):
        q_target = q_pert - self.q_points[iq1]
        found = False
        for iq2 in range(self.n_q):
            if CC.Methods.get_min_dist_into_cell(bg, q_target, self.q_points[iq2]) < 1e-6:
                self.q_pair_map[iq1] = iq2
                found = True
                break
        if not found:
            raise ValueError(
                "Could not find partner for q1={} with q_pert={}".format(
                    self.q_points[iq1], q_pert))

    # Unique pairs: iq1 <= iq2 (avoids double-counting)
    self.unique_pairs = []
    for iq1 in range(self.n_q):
        iq2 = self.q_pair_map[iq1]
        if iq1 <= iq2:
            self.unique_pairs.append((iq1, iq2))

    # Pre-compute block layout
    self._compute_block_layout()

get_psi_size()

Return the total size of the psi vector.

Source code in tdscha/QSpaceLanczos.py

def get_psi_size(self):
    """Return the total size of the psi vector."""
    if self._psi_size is None:
        raise ValueError("Must call build_q_pair_map first")
    return self._psi_size

get_static_psi_size()

Return psi size for the static layout: [R, one W sector].

This equals the end of the a' sector, i.e. the start of b'.

Source code in tdscha/QSpaceLanczos.py

def get_static_psi_size(self):
    """Return psi size for the static layout: [R, one W sector].

    This equals the end of the a' sector, i.e. the start of b'.
    """
    return self._block_offsets_b[0]

get_block_offset(pair_idx, sector='a')

Get the offset into psi for a given pair index.

Parameters:

Name	Type	Description	Default
pair_idx	int	Index into self.unique_pairs.	required
sector	str	'a' for a' sector, 'b' for b' sector.	'a'

Source code in tdscha/QSpaceLanczos.py

def get_block_offset(self, pair_idx, sector='a'):
    """Get the offset into psi for a given pair index.

    Parameters
    ----------
    pair_idx : int
        Index into self.unique_pairs.
    sector : str
        'a' for a' sector, 'b' for b' sector.
    """
    if sector == 'a':
        return self._block_offsets_a[pair_idx]
    else:
        return self._block_offsets_b[pair_idx]

get_block_size(pair_idx)

Get the number of entries for this pair.

Source code in tdscha/QSpaceLanczos.py

def get_block_size(self, pair_idx):
    """Get the number of entries for this pair."""
    return self._block_sizes[pair_idx]

get_R1_q()

Extract R^(1) from psi (n_bands complex entries at q_pert).

Source code in tdscha/QSpaceLanczos.py

def get_R1_q(self):
    """Extract R^(1) from psi (n_bands complex entries at q_pert)."""
    return self.psi[:self.n_bands].copy()

get_block(pair_idx, sector='a', source=None)

Reconstruct full (n_bands, n_bands) matrix from psi storage.

Parameters:

Name	Type	Description	Default
pair_idx	int	Index into self.unique_pairs.	required
sector	str	'a' or 'b'.	'a'
source	ndarray or None	If provided, read from this array instead of self.psi.	None

Returns:

Type	Description
(ndarray(n_bands, n_bands), complex128)

Source code in tdscha/QSpaceLanczos.py

def get_block(self, pair_idx, sector='a', source=None):
    """Reconstruct full (n_bands, n_bands) matrix from psi storage.

    Parameters
    ----------
    pair_idx : int
        Index into self.unique_pairs.
    sector : str
        'a' or 'b'.
    source : ndarray or None
        If provided, read from this array instead of self.psi.

    Returns
    -------
    ndarray(n_bands, n_bands), complex128
    """
    nb = self.n_bands
    iq1, iq2 = self.unique_pairs[pair_idx]
    offset = self.get_block_offset(pair_idx, sector)
    size = self.get_block_size(pair_idx)

    data = self.psi if source is None else source
    raw = data[offset:offset + size]

    if iq1 < iq2:
        # Full block, row-major
        return raw.reshape(nb, nb).copy()
    else:
        # Upper triangle storage
        return self._unpack_upper_triangle(raw, nb)

get_a1_block(pair_idx)

Get the a'(1) block for pair_idx.

Source code in tdscha/QSpaceLanczos.py

def get_a1_block(self, pair_idx):
    """Get the a'(1) block for pair_idx."""
    return self.get_block(pair_idx, 'a')

get_b1_block(pair_idx)

Get the b'(1) block for pair_idx.

Source code in tdscha/QSpaceLanczos.py

def get_b1_block(self, pair_idx):
    """Get the b'(1) block for pair_idx."""
    return self.get_block(pair_idx, 'b')

set_block_in_psi(pair_idx, matrix, sector, target_psi)

Write a (n_bands, n_bands) block into the target psi vector.

Parameters:

Name	Type	Description	Default
pair_idx	int		required
matrix	ndarray(n_bands, n_bands)		required
sector	str(a or b)		required
target_psi	ndarray — the psi vector to write into		required

Source code in tdscha/QSpaceLanczos.py

def set_block_in_psi(self, pair_idx, matrix, sector, target_psi):
    """Write a (n_bands, n_bands) block into the target psi vector.

    Parameters
    ----------
    pair_idx : int
    matrix : ndarray(n_bands, n_bands)
    sector : str ('a' or 'b')
    target_psi : ndarray — the psi vector to write into
    """
    nb = self.n_bands
    iq1, iq2 = self.unique_pairs[pair_idx]
    offset = self.get_block_offset(pair_idx, sector)

    if iq1 < iq2:
        # Full block
        target_psi[offset:offset + nb * nb] = matrix.ravel()
    else:
        # Upper triangle
        target_psi[offset:offset + self.get_block_size(pair_idx)] = \
            self._pack_upper_triangle(matrix, nb)

mask_dot_wigner(debug=False)

Build the mask for Hermitian inner product with upper-triangle storage.

For full blocks (iq1 < iq2): factor 2 for the conjugate block (iq2, iq1). For diagonal blocks (iq1 == iq2): off-diagonal factor 2, diagonal factor 1.

Returns:

Type	Description
(ndarray(psi_size), float64)

Source code in tdscha/QSpaceLanczos.py

def mask_dot_wigner(self, debug=False):
    """Build the mask for Hermitian inner product with upper-triangle storage.

    For full blocks (iq1 < iq2): factor 2 for the conjugate block (iq2, iq1).
    For diagonal blocks (iq1 == iq2): off-diagonal factor 2, diagonal factor 1.

    Returns
    -------
    ndarray(psi_size,), float64
    """
    mask = np.ones(self.get_psi_size(), dtype=np.float64)

    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        offset_a = self.get_block_offset(pair_idx, 'a')
        offset_b = self.get_block_offset(pair_idx, 'b')
        size = self.get_block_size(pair_idx)

        if iq1 < iq2:
            # Full block: factor 2 for the conjugate block
            mask[offset_a:offset_a + size] = 2
            mask[offset_b:offset_b + size] = 2
        else:  # iq1 == iq2, upper triangle storage
            nb = self.n_bands
            idx_a = offset_a
            idx_b = offset_b
            for i in range(nb):
                # Diagonal element: factor 1
                mask[idx_a] = 1
                mask[idx_b] = 1
                idx_a += 1
                idx_b += 1
                # Off-diagonal elements: factor 2
                for j in range(i + 1, nb):
                    mask[idx_a] = 2
                    mask[idx_b] = 2
                    idx_a += 1
                    idx_b += 1

    return mask

apply_L1_FT(transpose=False)

Apply the harmonic part of L in q-space (Wigner formalism).

L_harm is block-diagonal: R sector: -(w_q_pert[nu])^2 * R[nu] a' sector: -(w1 - w2)^2 * a' b' sector: -(w1 + w2)^2 * b'

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_L1_FT(self, transpose=False):
    """Apply the harmonic part of L in q-space (Wigner formalism).

    L_harm is block-diagonal:
      R sector: -(w_q_pert[nu])^2 * R[nu]
      a' sector: -(w1 - w2)^2 * a'
      b' sector: -(w1 + w2)^2 * b'

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    out = np.zeros(self.get_psi_size(), dtype=np.complex128)

    if self.ignore_harmonic:
        return out

    # R sector — mask acoustic modes
    w_qp = self.w_q[:, self.iq_pert]  # (n_bands,)
    v_pert = self.valid_modes_q[:, self.iq_pert]
    out[:self.n_bands] = -(w_qp ** 2) * self.psi[:self.n_bands]
    out[:self.n_bands] *= v_pert  # zero out acoustic

    # a' and b' sectors — mask pairs involving acoustic modes
    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        w1 = self.w_q[:, iq1]  # (n_bands,)
        w2 = self.w_q[:, iq2]  # (n_bands,)
        v1 = self.valid_modes_q[:, iq1]
        v2 = self.valid_modes_q[:, iq2]
        valid_mask = np.outer(v1, v2)

        a_block = self.get_a1_block(pair_idx)
        b_block = self.get_b1_block(pair_idx)

        w_minus2 = np.subtract.outer(w1, w2) ** 2  # (n_bands, n_bands)
        w_plus2 = np.add.outer(w1, w2) ** 2

        self.set_block_in_psi(pair_idx, -w_minus2 * a_block * valid_mask, 'a', out)
        self.set_block_in_psi(pair_idx, -w_plus2 * b_block * valid_mask, 'b', out)

    return out

get_chi_minus_q()

Get chi^- for each unique pair as a list of (n_bands, n_bands) matrices.

chi^-_{nu1, nu2} = (w1 - w2)(n1 - n2) / (2 * w1 * w2) Entries involving acoustic modes (w < acoustic_eps) are set to 0.

Source code in tdscha/QSpaceLanczos.py

def get_chi_minus_q(self):
    """Get chi^- for each unique pair as a list of (n_bands, n_bands) matrices.

    chi^-_{nu1, nu2} = (w1 - w2)(n1 - n2) / (2 * w1 * w2)
    Entries involving acoustic modes (w < acoustic_eps) are set to 0.
    """
    chi_list = []
    for iq1, iq2 in self.unique_pairs:
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        # Outer mask: both bands must be non-acoustic
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        # Safe division: use np.where to avoid 0/0
        denom = 2.0 * w1_mat * w2_mat
        chi = np.where(valid_mask,
                       (w1_mat - w2_mat) * (n1_mat - n2_mat) / np.where(valid_mask, denom, 1.0),
                       0.0)
        chi_list.append(chi)
    return chi_list

get_chi_plus_q()

Get chi^+ for each unique pair as a list of (n_bands, n_bands) matrices.

chi^+_{nu1, nu2} = (w1 + w2)(1 + n1 + n2) / (2 * w1 * w2) Entries involving acoustic modes (w < acoustic_eps) are set to 0.

Source code in tdscha/QSpaceLanczos.py

def get_chi_plus_q(self):
    """Get chi^+ for each unique pair as a list of (n_bands, n_bands) matrices.

    chi^+_{nu1, nu2} = (w1 + w2)(1 + n1 + n2) / (2 * w1 * w2)
    Entries involving acoustic modes (w < acoustic_eps) are set to 0.
    """
    chi_list = []
    for iq1, iq2 in self.unique_pairs:
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        denom = 2.0 * w1_mat * w2_mat
        chi = np.where(valid_mask,
                       (w1_mat + w2_mat) * (1 + n1_mat + n2_mat) / np.where(valid_mask, denom, 1.0),
                       0.0)
        chi_list.append(chi)
    return chi_list

get_alpha1_beta1_wigner_q(get_alpha=True)

Get the perturbation on alpha (Upsilon) from the q-space psi.

Transforms a'/b' blocks back to the alpha1 perturbation that the Julia code needs.

alpha1[iq1, iq2] = (w1w2/X) * [sqrt(-0.5chi_minus)a' - sqrt(0.5chi_plus)*b']

Returns:

Type	Description
list of ndarray(n_bands, n_bands) — one per unique pair

Source code in tdscha/QSpaceLanczos.py

def get_alpha1_beta1_wigner_q(self, get_alpha=True):
    """Get the perturbation on alpha (Upsilon) from the q-space psi.

    Transforms a'/b' blocks back to the alpha1 perturbation that the
    Julia code needs.

    alpha1[iq1, iq2] = (w1*w2/X) * [sqrt(-0.5*chi_minus)*a' - sqrt(0.5*chi_plus)*b']

    Returns
    -------
    list of ndarray(n_bands, n_bands) — one per unique pair
    """
    chi_minus_list = self.get_chi_minus_q()
    chi_plus_list = self.get_chi_plus_q()

    alpha1_blocks = []
    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        X = (1 + 2 * n1_mat) * (1 + 2 * n2_mat) / 8
        # Safe division for w2_on_X: when acoustic, X→∞ and w→0, set to 0
        w2_on_X = np.where(valid_mask,
                           (w1_mat * w2_mat) / np.where(valid_mask, X, 1.0),
                           0.0)
        chi_minus = chi_minus_list[pair_idx]
        chi_plus = chi_plus_list[pair_idx]

        a_block = self.get_a1_block(pair_idx)
        b_block = self.get_b1_block(pair_idx)

        if get_alpha:
            new_a = w2_on_X * np.sqrt(-0.5 * chi_minus) * a_block
            new_b = w2_on_X * np.sqrt(+0.5 * chi_plus) * b_block
            alpha1 = new_a - new_b
        else:
            X_safe = np.where(valid_mask, X, 1.0)
            new_a = np.where(valid_mask,
                             (np.sqrt(-0.5 * chi_minus) / X_safe) * a_block,
                             0.0)
            new_b = np.where(valid_mask,
                             (np.sqrt(+0.5 * chi_plus) / X_safe) * b_block,
                             0.0)
            alpha1 = new_a + new_b

        alpha1_blocks.append(alpha1)
    return alpha1_blocks

apply_anharmonic_FT(transpose=False, **kwargs)

Apply the anharmonic part of L in q-space (Wigner formalism).

Calls the Julia q-space extension to compute the perturbed averages, then assembles the output psi vector.

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_anharmonic_FT(self, transpose=False, **kwargs):
    """Apply the anharmonic part of L in q-space (Wigner formalism).

    Calls the Julia q-space extension to compute the perturbed averages,
    then assembles the output psi vector.

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    # If both D3 and D4 are ignored, return zero
    if self.ignore_v3 and self.ignore_v4:
        return np.zeros(self.get_psi_size(), dtype=np.complex128)

    import julia.Main

    R1 = self.get_R1_q()
    # If D3 is ignored, zero out R1 so that D3 weight is zero
    if self.ignore_v3:
        R1 = np.zeros_like(R1)
    alpha1_blocks = self.get_alpha1_beta1_wigner_q(get_alpha=True)
    alpha1_flat = self._flatten_blocks(alpha1_blocks)

    # Call Julia
    f_pert, d2v_blocks = self._call_julia_qspace(R1, alpha1_flat)

    # Build output psi
    final_psi = np.zeros(self.get_psi_size(), dtype=np.complex128)

    # R sector
    final_psi[:self.n_bands] = f_pert

    # a'/b' sectors
    chi_minus_list = self.get_chi_minus_q()
    chi_plus_list = self.get_chi_plus_q()

    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        d2v_block = d2v_blocks[pair_idx]
        pert_a = np.sqrt(-0.5 * chi_minus_list[pair_idx]) * d2v_block
        pert_b = -np.sqrt(+0.5 * chi_plus_list[pair_idx]) * d2v_block
        self.set_block_in_psi(pair_idx, pert_a, 'a', final_psi)
        self.set_block_in_psi(pair_idx, pert_b, 'b', final_psi)

    return final_psi

apply_full_L(target=None, force_t_0=False, force_FT=True, transpose=False, fast_lanczos=True)

Apply the full L operator in q-space.

Parameters:

Name	Type	Description	Default
target	ndarray or None	If provided, copy into self.psi first.	None
transpose	bool	Not used for Hermitian Lanczos.	False

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_full_L(self, target=None, force_t_0=False, force_FT=True,
                 transpose=False, fast_lanczos=True):
    """Apply the full L operator in q-space.

    Parameters
    ----------
    target : ndarray or None
        If provided, copy into self.psi first.
    transpose : bool
        Not used for Hermitian Lanczos.

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    if target is not None:
        self.psi = target.copy()

    result = self.apply_L1_FT(transpose=transpose)
    result += self.apply_anharmonic_FT(transpose=transpose)

    return result

run_FT(n_iter, save_dir=None, save_each=5, verbose=True, n_rep_orth=0, n_ortho=10, flush_output=True, debug=False, prefix='LANCZOS', run_simm=None, optimized=False, reorthogonalize=True)

Run the Hermitian Lanczos algorithm for q-space.

This is the same structure as the parent run_FT but with: 1. Forced run_simm = True (Hermitian) 2. Hermitian dot products: psi.conj().dot(psi * mask).real 3. Complex128 psi 4. Real coefficients (guaranteed by Hermitian L)

Source code in tdscha/QSpaceLanczos.py

    def run_FT(self, n_iter, save_dir=None, save_each=5, verbose=True,
               n_rep_orth=0, n_ortho=10, flush_output=True, debug=False,
               prefix="LANCZOS", run_simm=None, optimized=False,
               reorthogonalize=True):
        """Run the Hermitian Lanczos algorithm for q-space.

        This is the same structure as the parent run_FT but with:
        1. Forced run_simm = True (Hermitian)
        2. Hermitian dot products: psi.conj().dot(psi * mask).real
        3. Complex128 psi
        4. Real coefficients (guaranteed by Hermitian L)
        """
        self.verbose = verbose

        if not self.initialized:
            if verbose:
                print('Not initialized. Now we symmetrize\n')
            self.prepare_symmetrization()

        ERROR_MSG = """
Error, you must initialize a perturbation to start the Lanczos.
Use prepare_mode_q or prepare_perturbation_q before calling run_FT.
"""
        if self.psi is None:
            raise ValueError(ERROR_MSG)

        mask_dot = self.mask_dot_wigner(debug)
        psi_norm = np.real(np.conj(self.psi).dot(self.psi * mask_dot))
        if np.isnan(psi_norm) or psi_norm == 0:
            raise ValueError(ERROR_MSG)

        # Force symmetric Lanczos
        run_simm = True

        if Parallel.am_i_the_master():
            if save_dir is not None:
                if not os.path.exists(save_dir):
                    os.makedirs(save_dir)

        if verbose:
            print('Running the Hermitian Lanczos in q-space')
            print()

        # Get current step
        i_step = len(self.a_coeffs)

        if verbose:
            header = """
<=====================================>
|                                     |
|     Q-SPACE LANCZOS ALGORITHM       |
|                                     |
<=====================================>

Starting the algorithm.
Starting from step %d
""" % i_step
            print(header)

        # Initialize
        if i_step == 0:
            self.basis_Q = []
            self.basis_P = []
            self.s_norm = []
            norm = np.sqrt(np.real(np.conj(self.psi).dot(self.psi * mask_dot)))
            first_vector = self.psi / norm
            self.basis_Q.append(first_vector)
            self.basis_P.append(first_vector)
            self.s_norm.append(1)
        else:
            if verbose:
                print('Restarting the Lanczos')
            self.basis_Q = list(self.basis_Q)
            self.basis_P = list(self.basis_P)
            self.s_norm = list(self.s_norm)
            self.a_coeffs = list(self.a_coeffs)
            self.b_coeffs = list(self.b_coeffs)
            self.c_coeffs = list(self.c_coeffs)

        psi_q = self.basis_Q[-1]
        psi_p = self.basis_P[-1]

        next_converged = False
        converged = False

        for i in range(i_step, i_step + n_iter):
            if verbose:
                print("\n ===== NEW STEP %d =====\n" % (i + 1))
                if flush_output:
                    sys.stdout.flush()

            # Apply L (Hermitian => p_L = L_q)
            t1 = time.time()
            L_q = self.apply_full_L(psi_q)
            p_L = np.copy(L_q)
            t2 = time.time()

            # p normalization
            c_old = 1
            if len(self.c_coeffs) > 0:
                c_old = self.c_coeffs[-1]
            p_norm = self.s_norm[-1] / c_old

            # a coefficient (real for Hermitian L)
            a_coeff = np.real(np.conj(psi_p).dot(L_q * mask_dot)) * p_norm

            if np.isnan(a_coeff):
                raise ValueError("Invalid value in Lanczos. Check frequencies/initialization.")

            # Residuals
            rk = L_q - a_coeff * psi_q
            if len(self.basis_Q) > 1:
                rk -= self.c_coeffs[-1] * self.basis_Q[-2]

            sk = p_L - a_coeff * psi_p
            if len(self.basis_P) > 1:
                old_p_norm = self.s_norm[-2]
                if len(self.c_coeffs) >= 2:
                    old_p_norm = self.s_norm[-2] / self.c_coeffs[-2]
                sk -= self.b_coeffs[-1] * self.basis_P[-2] * (old_p_norm / p_norm)

            # s_norm
            s_norm = np.sqrt(np.real(np.conj(sk).dot(sk * mask_dot)))
            sk_tilde = sk / s_norm
            s_norm *= p_norm

            # b and c coefficients (real, should be equal for Hermitian L)
            b_coeff = np.sqrt(np.real(np.conj(rk).dot(rk * mask_dot)))
            c_coeff = np.real(np.conj(sk_tilde).dot((rk / b_coeff) * mask_dot)) * s_norm

            self.a_coeffs.append(a_coeff)

            if np.abs(b_coeff) < __EPSILON__ or next_converged:
                if verbose:
                    print("Converged (b = {})".format(b_coeff))
                converged = True
                break
            if np.abs(c_coeff) < __EPSILON__:
                if verbose:
                    print("Converged (c = {})".format(c_coeff))
                converged = True
                break

            psi_q = rk / b_coeff
            psi_p = sk_tilde.copy()

            # Gram-Schmidt reorthogonalization
            if reorthogonalize:
                # Correct Hermitian MGS: orthogonalize against all Q vectors
                # (which are unit-normalized in the mask inner product)
                new_q = psi_q.copy()

                for j in range(len(self.basis_Q)):
                    coeff = np.real(np.conj(self.basis_Q[j]).dot(new_q * mask_dot))
                    new_q -= coeff * self.basis_Q[j]

                normq = np.sqrt(np.real(np.conj(new_q).dot(new_q * mask_dot)))
                if normq < __EPSILON__:
                    next_converged = True
                new_q /= normq

                # Hermitian L: P = Q, s_norm = c_coeff (since <P,Q> = <Q,Q> = 1)
                new_p = new_q.copy()
                s_norm = c_coeff
            else:
                # Existing biconjugate GS (for backward compatibility)
                new_q = psi_q.copy()
                new_p = psi_p.copy()

                for k_orth in range(n_rep_orth):
                    start = max(0, len(self.basis_P) - (n_ortho or len(self.basis_P)))

                    for j in range(start, len(self.basis_P)):
                        coeff1 = np.real(np.conj(self.basis_P[j]).dot(new_q * mask_dot))
                        coeff2 = np.real(np.conj(self.basis_Q[j]).dot(new_p * mask_dot))
                        new_q -= coeff1 * self.basis_P[j]
                        new_p -= coeff2 * self.basis_Q[j]

                    normq = np.sqrt(np.real(np.conj(new_q).dot(new_q * mask_dot)))
                    if normq < __EPSILON__:
                        next_converged = True
                    new_q /= normq

                    normp = np.real(np.conj(new_p).dot(new_p * mask_dot))
                    if np.abs(normp) < __EPSILON__:
                        next_converged = True
                    new_p /= normp

                    s_norm = c_coeff / np.real(np.conj(new_p).dot(new_q * mask_dot))

            if not converged:
                self.basis_Q.append(new_q)
                self.basis_P.append(new_p)
                psi_q = new_q.copy()
                psi_p = new_p.copy()

                self.b_coeffs.append(b_coeff)
                self.c_coeffs.append(c_coeff)
                self.s_norm.append(s_norm)

            if verbose:
                print("Time for L application: %d s" % (t2 - t1))
                print("a_%d = %.8e" % (i, self.a_coeffs[-1]))
                print("b_%d = %.8e" % (i, self.b_coeffs[-1]))
                print("c_%d = %.8e" % (i, self.c_coeffs[-1]))
                print("|b-c| = %.8e" % np.abs(self.b_coeffs[-1] - self.c_coeffs[-1]))

            if save_dir is not None:
                if (i + 1) % save_each == 0:
                    self.save_status("%s/%s_STEP%d" % (save_dir, prefix, i + 1))

            if verbose:
                print("Lanczos step %d completed." % (i + 1))

        if converged and verbose:
            print("   last a coeff = {}".format(self.a_coeffs[-1]))

prepare_mode_q(iq, band_index)

Prepare perturbation for mode (q, nu).

Parameters:

Name	Type	Description	Default
iq	int	Index of the q-point.	required
band_index	int	Band index (0-based).	required

Source code in tdscha/QSpaceLanczos.py

def prepare_mode_q(self, iq, band_index):
    """Prepare perturbation for mode (q, nu).

    Parameters
    ----------
    iq : int
        Index of the q-point.
    band_index : int
        Band index (0-based).
    """
    if band_index < 0 or band_index >= self.n_bands:
        raise ValueError("Invalid band index for perturbation: {}".format(band_index))

    self.build_q_pair_map(iq)
    self.reset_q()
    self.psi[band_index] = 1.0 + 0j
    self.perturbation_modulus = 1.0

prepare_ir(effective_charges=None, pol_vec=np.array([1.0, 0.0, 0.0]))

PREPARE LANCZOS FOR INFRARED SPECTRUM COMPUTATION

In this subroutine we prepare the lanczos algorithm for the computation of the infrared spectrum signal.

Source code in tdscha/QSpaceLanczos.py

def prepare_ir(self, effective_charges = None, pol_vec = np.array([1.0, 0.0, 0.0])):
    """
    PREPARE LANCZOS FOR INFRARED SPECTRUM COMPUTATION
    =================================================

    In this subroutine we prepare the lanczos algorithm for the computation of the
    infrared spectrum signal.

    Parameters
    ----------
        effective_charges : ndarray(size = (n_atoms, 3, 3), dtype = np.double)
            The effective charges. Indices are: Number of atoms in the unit cell,
            electric field component, atomic coordinate. If None, the effective charges
            contained in the dynamical matrix will be considered.
        pol_vec : ndarray(size = 3)
            The polarization vector of the light.
    """

    ec = self.dyn.effective_charges
    if not effective_charges is None:
        ec = effective_charges

    assert not ec is None, "Error, no effective charge found. Cannot initialize IR responce"

    z_eff = np.einsum("abc, b", ec, pol_vec)

    # Get the gamma effective charge
    # FIX: remove double mass scaling and add supercell factor
    n_cell = np.prod(self.dyn.GetSupercell())
    new_zeff = z_eff.ravel() * np.sqrt(n_cell)

    # This is a Gamma perturbation
    self.prepare_perturbation_q(0, new_zeff)

prepare_raman(pol_vec_in=np.array([1.0, 0.0, 0.0]), pol_vec_out=np.array([1.0, 0.0, 0.0]), mixed=False, pol_in_2=None, pol_out_2=None, unpolarized=None)

PREPARE LANCZOS FOR RAMAN SPECTRUM COMPUTATION

In this subroutine we prepare the lanczos algorithm for the computation of the Raman spectrum signal.

Source code in tdscha/QSpaceLanczos.py

def prepare_raman(self, pol_vec_in=np.array([1.0, 0.0, 0.0]), pol_vec_out=np.array([1.0, 0.0, 0.0]), 
                 mixed=False, pol_in_2=None, pol_out_2=None, unpolarized=None):
    """
    PREPARE LANCZOS FOR RAMAN SPECTRUM COMPUTATION
    ==============================================

    In this subroutine we prepare the lanczos algorithm for the computation of the
    Raman spectrum signal.

    Parameters
    ----------
        pol_vec_in : ndarray(size = 3)
            The polarization vector of the incoming light
        pol_vec_out : ndarray(size = 3)
            The polarization vector for the outcoming light
        mixed : bool
            If True, add another component of the Raman tensor
        pol_in_2 : ndarray(size = 3) or None
            Second incoming polarization if mixed=True
        pol_out_2 : ndarray(size = 3) or None  
            Second outcoming polarization if mixed=True
        unpolarized : int or None
            The perturbation for unpolarized raman (if different from None, overrides the behaviour
            of pol_vec_in and pol_vec_out). Indices goes from 0 to 6 (included).
            0 is alpha^2
            1 + 2 + 3 + 4 + 5 + 6 are beta^2
            alpha_0 = (xx + yy + zz)^2/9
            beta_1 = (xx -yy)^2 / 2
            beta_2 = (xx -zz)^2 / 2
            beta_3 = (yy -zz)^2 / 2
            beta_4 = 3xy^2
            beta_5 = 3xz^2
            beta_6 = 3yz^2

            The total unpolarized raman intensity is 45 alpha^2 + 7 beta^2
    """
    # Check if the raman tensor is present
    assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize Raman response"

    n_cell = np.prod(self.dyn.GetSupercell())

    if unpolarized is None:
        # Get the raman vector (apply the ASR and contract the raman tensor with the polarization vectors)
        raman_v = self.dyn.GetRamanVector(pol_vec_in, pol_vec_out)

        if mixed:
            print('Prepare Raman')
            print('Adding other component of the Raman tensor')
            raman_v += self.dyn.GetRamanVector(pol_in_2, pol_out_2)

        # Scale for Γ-point constant perturbation
        new_raman_v = raman_v.ravel() * np.sqrt(n_cell)

        # Convert in the polarization basis
        self.prepare_perturbation_q(0, new_raman_v)
    else:
        px = np.array([1, 0, 0])
        py = np.array([0, 1, 0])
        pz = np.array([0, 0, 1])

        if unpolarized == 0:
            # Alpha = (xx + yy + zz)^2/9
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v, add=True)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 1:
            # (xx - yy)^2 / 2
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 2:
            # (xx - zz)^2 / 2
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 3:
            # (yy - zz)^2 / 2
            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 4:
            # 3 xy^2
            raman_v = self.dyn.GetRamanVector(px, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        elif unpolarized == 5:
            # 3 yz^2
            raman_v = self.dyn.GetRamanVector(py, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        elif unpolarized == 6:
            # 3 xz^2
            raman_v = self.dyn.GetRamanVector(px, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        else:
            raise ValueError(f"Error, unpolarized must be between [0, ..., 6] got invalid {unpolarized}.")

prepare_unpolarized_raman(index=0, debug=False)

PREPARE UNPOLARIZED RAMAN SIGNAL

The raman tensor is read from the dynamical matrix provided by the original ensemble.

The perturbations are prepared according to the formula (see https://doi.org/10.1021/jp5125266)

..math:

I_unpol = 45/9 (xx + yy + zz)^2
          + 7/2 [(xx-yy)^2 + (xx-zz)^2 + (yy-zz)^2]
          + 7 * 3 [(xy)^2 + (yz)^2 + (xz)^2]

Note: This method prepares the raw components WITHOUT prefactors. Use get_prefactors_unpolarized_raman() to get the correct prefactors.

Source code in tdscha/QSpaceLanczos.py

def prepare_unpolarized_raman(self, index=0, debug=False):
    """
    PREPARE UNPOLARIZED RAMAN SIGNAL
    ================================

    The raman tensor is read from the dynamical matrix provided by the original ensemble.

    The perturbations are prepared according to the formula (see https://doi.org/10.1021/jp5125266)

    ..math:

        I_unpol = 45/9 (xx + yy + zz)^2
                  + 7/2 [(xx-yy)^2 + (xx-zz)^2 + (yy-zz)^2]
                  + 7 * 3 [(xy)^2 + (yz)^2 + (xz)^2]

    Note: This method prepares the raw components WITHOUT prefactors.
    Use get_prefactors_unpolarized_raman() to get the correct prefactors.
    """
    # Check if the raman tensor is present
    assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize the Raman response"

    labels = [i for i in range(7)]
    if index not in labels:
        raise ValueError(f'{index} should be in {labels}')

    epols = {'x': np.array([1, 0, 0]),
             'y': np.array([0, 1, 0]),
             'z': np.array([0, 0, 1])}

    n_cell = np.prod(self.dyn.GetSupercell())

    # (xx + yy + zz)^2
    if index == 0:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v += self.dyn.GetRamanVector(epols['y'], epols['y'])
        raman_v += self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (xx - yy)^2    
    elif index == 1:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v -= self.dyn.GetRamanVector(epols['y'], epols['y'])
    # (xx - zz)^2       
    elif index == 2:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (yy - zz)^2   
    elif index == 3:
        raman_v = self.dyn.GetRamanVector(epols['y'], epols['y'])
        raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (xy)^2
    elif index == 4:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['y'])
    # (xz)^2
    elif index == 5:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['z'])
    # (yz)^2
    elif index == 6:
        raman_v = self.dyn.GetRamanVector(epols['y'], epols['z'])

    if debug:
        np.save(f'raman_v_{index}', raman_v)

    # Scale for Γ-point constant perturbation
    new_raman_v = raman_v.ravel() * np.sqrt(n_cell)

    # Convert in the polarization basis
    self.prepare_perturbation_q(0, new_raman_v)

    if debug:
        print(f'[NEW] Perturbation modulus with eq Raman tensors = {self.perturbation_modulus}')
    print()

prepare_perturbation_q(iq, vector, add=False)

Prepare perturbation at q from a real-space vector (3*n_at_uc,).

Projects the vector onto q-space eigenmodes at iq.

Parameters:

Name	Type	Description	Default
iq	int	Index of the q-point.	required
vector	ndarray(3 * n_at_uc)	Perturbation vector in Cartesian real space.	required
add	bool	If true, the perturbation is added on top of the one already setup. Calling add does not cause a reset of the Lanczos.	False

Source code in tdscha/QSpaceLanczos.py

def prepare_perturbation_q(self, iq, vector, add=False):
    """Prepare perturbation at q from a real-space vector (3*n_at_uc,).

    Projects the vector onto q-space eigenmodes at iq.

    Parameters
    ----------
    iq : int
        Index of the q-point.
    vector : ndarray(3*n_at_uc,)
        Perturbation vector in Cartesian real space.
    add : bool
        If true, the perturbation is added on top of the one already setup.
        Calling add does not cause a reset of the Lanczos.
    """
    if not add:
        self.build_q_pair_map(iq)
        self.reset_q()

    m = np.tile(self.uci_structure.get_masses_array(), (3, 1)).T.ravel()
    v_scaled = vector / np.sqrt(m)
    R1 = np.conj(self.pols_q[:, :, iq]).T @ v_scaled  # (n_bands,) complex
    self.psi[:self.n_bands] += R1
    self.perturbation_modulus = np.real(np.conj(R1) @ R1)

reset_q()

Reset the Lanczos state for q-space.

Source code in tdscha/QSpaceLanczos.py

def reset_q(self):
    """Reset the Lanczos state for q-space."""
    n = self.get_psi_size()
    self.psi = np.zeros(n, dtype=np.complex128)

    self.eigvals = None
    self.eigvects = None

    self.a_coeffs = []
    self.b_coeffs = []
    self.c_coeffs = []
    self.krilov_basis = []
    self.basis_P = []
    self.basis_Q = []
    self.s_norm = []

prepare_symmetrization(no_sym=False, verbose=True, symmetries=None)

Build q-space symmetry matrices and cache them in Julia.

Overrides the parent to build sparse complex symmetry matrices in the q-space mode basis.

Uses spglib on the unit cell (not supercell) to get correct fractional-coordinate rotations and translations, then converts to Cartesian for the representation matrices.

Source code in tdscha/QSpaceLanczos.py

def prepare_symmetrization(self, no_sym=False, verbose=True, symmetries=None):
    """Build q-space symmetry matrices and cache them in Julia.

    Overrides the parent to build sparse complex symmetry matrices
    in the q-space mode basis.

    Uses spglib on the unit cell (not supercell) to get correct
    fractional-coordinate rotations and translations, then converts
    to Cartesian for the representation matrices.
    """
    self.initialized = True

    if no_sym:
        # Identity only
        self.n_syms_qspace = 1
        n_total = self.n_q * self.n_bands
        # Build identity sparse matrix
        import julia.Main
        julia.Main.eval("""
        function init_identity_qspace(n_total::Int64)
            I_sparse = SparseArrays.sparse(
                Int32.(1:n_total), Int32.(1:n_total),
                ComplexF64.(ones(n_total)), n_total, n_total)
            _cached_qspace_symmetries[] = [I_sparse]
            return nothing
        end
        """)
        julia.Main.init_identity_qspace(int(n_total))
        return

    if not __SPGLIB__:
        raise ImportError("spglib required for symmetrization")

    # Get symmetries from the UNIT CELL directly via spglib.
    # spglib returns rotations and translations in fractional
    # (crystal) coordinates. We convert rotations to Cartesian via
    # R_cart = M @ R_frac @ M^{-1} where M = unit_cell.T.
    spg_data = spglib.get_symmetry(self.uci_structure.get_spglib_cell())
    rot_frac_all = spg_data['rotations']   # (n_sym, 3, 3) integer
    trans_frac_all = spg_data['translations']  # (n_sym, 3) fractional

    M = self.uci_structure.unit_cell.T       # columns = lattice vectors
    Minv = np.linalg.inv(M)

    # Extract unique point-group rotations (keep first occurrence)
    unique_pg = {}
    for i in range(len(rot_frac_all)):
        key = rot_frac_all[i].tobytes()
        if key not in unique_pg:
            unique_pg[key] = i
    pg_indices = list(unique_pg.values())

    if verbose:
        print("Q-space: {} PG symmetries from {} total unit cell symmetries".format(
            len(pg_indices), len(rot_frac_all)))

    self._build_qspace_symmetries(
        rot_frac_all, trans_frac_all, pg_indices, M, Minv,
        verbose=verbose)

init(use_symmetries=True)

Initialize the q-space Lanczos calculation.

Source code in tdscha/QSpaceLanczos.py

def init(self, use_symmetries=True):
    """Initialize the q-space Lanczos calculation."""
    self.prepare_symmetrization(no_sym=not use_symmetries)
    self.initialized = True

find_q_index(q_target, q_points, bg, tol=1e-06)

Find the index of q_target in q_points up to a reciprocal lattice vector.

Parameters:

Name	Type	Description	Default
q_target	ndarray(3)	The q-point to find (Cartesian coordinates).	required
q_points	ndarray(n_q, 3)	Array of q-points.	required
bg	ndarray(3, 3)	Reciprocal lattice vectors / (2*pi), rows are vectors.	required
tol	float	Tolerance for matching.	1e-06

Returns:

Type	Description
int	Index of the matching q-point.

Source code in tdscha/QSpaceLanczos.py

def find_q_index(q_target, q_points, bg, tol=1e-6):
    """Find the index of q_target in q_points up to a reciprocal lattice vector.

    Parameters
    ----------
    q_target : ndarray(3,)
        The q-point to find (Cartesian coordinates).
    q_points : ndarray(n_q, 3)
        Array of q-points.
    bg : ndarray(3, 3)
        Reciprocal lattice vectors / (2*pi), rows are vectors.
    tol : float
        Tolerance for matching.

    Returns
    -------
    int
        Index of the matching q-point.
    """
    for iq, q in enumerate(q_points):
        dist = CC.Methods.get_min_dist_into_cell(bg, q_target, q)
        if dist < tol:
            return iq
    raise ValueError("Could not find q-point {} in the q-point list".format(q_target))

load_distributed_tdscha(data_dir, population_id, dyn, T, lo_to_split=None, use_symmetries=True, n_configs=None, final_dyn=None, final_T=None, **kwargs)

Load QSpaceLanczos with distributed configurations across MPI ranks.

Loads the ensemble on master rank only, then distributes configuration data across all ranks. Each process stores only N/n_procs configs instead of all N. This avoids having the full ensemble replicated in memory on all ranks.

Parameters:

Name	Type	Description	Default
data_dir	str	Directory containing the ensemble data files.	required
population_id	int	Population ID of the ensemble to load.	required
dyn	Phonons	Dynamical matrix object (used to create the ensemble).	required
T	float	Temperature in Kelvin.	required
lo_to_split	str, ndarray, or None	LO-TO splitting mode.	None
use_symmetries	bool	If True, use q-space symmetries.	True
n_configs	int or None	Number of configs to load. If None, loads all available.	None
final_dyn	Phonons	Final dynamical matrix from the SSCHA calculation. If provided, the ensemble weights are updated using update_weights(final_dyn, final_T). This is recommended for production calculations where the ensemble was generated with a preliminary dynamical matrix.	None
final_T	float	Temperature for weight updates. Defaults to T if not specified. Use this if the final temperature differs from the ensemble temperature.	None
**kwargs		Additional arguments passed to QSpaceLanczos.	{}

Returns:

Type	Description
QSpaceLanczos	QSpaceLanczos with distributed configs. Each rank has N_local = N/n_procs configs. With n_procs=1, rank 0 holds all configs (N_local = N).

Usage

mpirun -np 8 python your_script.py

Flow: - Rank 0: loads ensemble, optionally updates weights, creates QSpaceLanczos, broadcasts metadata, sends slices - Ranks 1..n-1: receive metadata, build bare QSpaceLanczos, receive slices - All: have only their local slice, _distributed=True

Source code in tdscha/QSpaceLanczos.py

def load_distributed_tdscha(data_dir, population_id, dyn, T, lo_to_split=None,
                           use_symmetries=True, n_configs=None,
                           final_dyn=None, final_T=None, **kwargs):
    """Load QSpaceLanczos with distributed configurations across MPI ranks.

    Loads the ensemble on master rank only, then distributes configuration data
    across all ranks. Each process stores only N/n_procs configs instead of all N.
    This avoids having the full ensemble replicated in memory on all ranks.

    Parameters
    ----------
    data_dir : str
        Directory containing the ensemble data files.
    population_id : int
        Population ID of the ensemble to load.
    dyn : CC.Phonons.Phonons
        Dynamical matrix object (used to create the ensemble).
    T : float
        Temperature in Kelvin.
    lo_to_split : str, ndarray, or None
        LO-TO splitting mode.
    use_symmetries : bool
        If True, use q-space symmetries.
    n_configs : int or None
        Number of configs to load. If None, loads all available.
    final_dyn : CC.Phonons.Phonons, optional
        Final dynamical matrix from the SSCHA calculation. If provided,
        the ensemble weights are updated using update_weights(final_dyn, final_T).
        This is recommended for production calculations where the ensemble
        was generated with a preliminary dynamical matrix.
    final_T : float, optional
        Temperature for weight updates. Defaults to T if not specified.
        Use this if the final temperature differs from the ensemble temperature.
    **kwargs
        Additional arguments passed to QSpaceLanczos.

    Returns
    -------
    QSpaceLanczos
        QSpaceLanczos with distributed configs. Each rank has N_local = N/n_procs configs.
        With n_procs=1, rank 0 holds all configs (N_local = N).

    Usage
    -----
    mpirun -np 8 python your_script.py

    Flow:
    - Rank 0: loads ensemble, optionally updates weights, creates QSpaceLanczos,
              broadcasts metadata, sends slices
    - Ranks 1..n-1: receive metadata, build bare QSpaceLanczos, receive slices
    - All: have only their local slice, _distributed=True
    """
    comm = mpi4py.MPI.COMM_WORLD
    rank = comm.Get_rank()
    n_procs = Parallel.GetNProc()

    if Parallel.am_i_the_master():
        # ========== MASTER (RANK 0) ==========
        ensemble = sscha.Ensemble.Ensemble(dyn, T)
        if n_configs is not None:
            ensemble.load_bin(data_dir, population_id, n_configs=n_configs)
        else:
            ensemble.load_bin(data_dir, population_id)

        # Update weights if final_dyn is provided
        if final_dyn is not None:
            T_for_update = final_T if final_T is not None else T
            ensemble.update_weights(final_dyn, T_for_update)

        qlanc = QSpaceLanczos(ensemble, lo_to_split=lo_to_split, **kwargs)
        qlanc.init(use_symmetries=use_symmetries)

        # Free ensemble - we only need QSpaceLanczos arrays
        del ensemble

        # Extract global info
        N_global = qlanc.N
        N_eff_global = qlanc.N_eff

        # Broadcast metadata (structure arrays, NO config data)
        metadata = {
            'T': qlanc.T, 'dyn': qlanc.dyn,
            'uci_structure': qlanc.uci_structure,
            'super_structure': qlanc.super_structure,
            'q_points': qlanc.q_points, 'n_q': qlanc.n_q,
            'n_bands': qlanc.n_bands, 'w_q': qlanc.w_q,
            'pols_q': qlanc.pols_q, 'valid_modes_q': qlanc.valid_modes_q,
            'm': qlanc.m,
            'ignore_v3': qlanc.ignore_v3, 'ignore_v4': qlanc.ignore_v4,
            'N_global': N_global, 'N_eff_global': N_eff_global,
            'n_syms_qspace': qlanc.n_syms_qspace,
            '_qspace_sym_data': qlanc._qspace_sym_data,
            '_qspace_sym_q_map': qlanc._qspace_sym_q_map,
        }
        comm.bcast(metadata, root=0)

        # Barrier to ensure all ranks have received metadata before we start sending slices
        comm.barrier()

        # Distribute configs to other ranks
        configs_per_proc = N_global // n_procs
        remainder = N_global % n_procs

        for target in range(1, n_procs):
            if target < remainder:
                start = target * (configs_per_proc + 1)
                end = start + configs_per_proc + 1
            else:
                start = target * configs_per_proc + remainder
                end = start + configs_per_proc

            N_local = end - start
            comm.Send(np.array([N_local], dtype=np.int64), dest=target, tag=0)
            comm.Send(qlanc.X_q[:, start:end, :].copy(), dest=target, tag=1)
            comm.Send(qlanc.Y_q[:, start:end, :].copy(), dest=target, tag=2)
            comm.Send(qlanc.rho[start:end].copy(), dest=target, tag=3)

        # Master keeps its slice (rank 0)
        if remainder > 0:
            start, end = 0, configs_per_proc + 1
        else:
            start, end = 0, configs_per_proc

        N_local = end - start
        qlanc.X_q = qlanc.X_q[:, start:end, :].copy()
        qlanc.Y_q = qlanc.Y_q[:, start:end, :].copy()
        qlanc.rho = qlanc.rho[start:end].copy()

        # Set distributed state
        qlanc._distributed = True
        qlanc._N_global = N_global
        qlanc._N_eff_global = N_eff_global
        qlanc._N_local = N_local
        qlanc.N = N_local
        qlanc.N_eff = int(np.sum(qlanc.rho))

        # Free unused arrays
        if hasattr(qlanc, 'X') and qlanc.X is not None:
            qlanc.X = None
        if hasattr(qlanc, 'Y') and qlanc.Y is not None:
            qlanc.Y = None

        return qlanc

    else:
        # ========== OTHER RANKS ==========
        metadata = comm.bcast(None, root=0)

        # Barrier to ensure all ranks have received metadata before slices are sent
        comm.barrier()

        # Create bare QSpaceLanczos and populate from metadata
        qlanc = QSpaceLanczos(ensemble=None, lo_to_split=lo_to_split, **kwargs)
        qlanc.T = metadata['T']
        qlanc.dyn = metadata['dyn']
        qlanc.uci_structure = metadata['uci_structure']
        qlanc.super_structure = metadata['super_structure']
        qlanc.q_points = metadata['q_points']
        qlanc.n_q = metadata['n_q']
        qlanc.n_bands = metadata['n_bands']
        qlanc.w_q = metadata['w_q']
        qlanc.pols_q = metadata['pols_q']
        qlanc.valid_modes_q = metadata['valid_modes_q']
        qlanc.m = metadata['m']
        qlanc.ignore_v3 = metadata['ignore_v3']
        qlanc.ignore_v4 = metadata['ignore_v4']
        qlanc.n_syms_qspace = metadata['n_syms_qspace']
        qlanc._qspace_sym_data = metadata['_qspace_sym_data']
        qlanc._qspace_sym_q_map = metadata['_qspace_sym_q_map']

        # Receive local config slice
        N_local_arr = np.array([0], dtype=np.int64)
        comm.Recv(N_local_arr, source=0, tag=0)
        N_local = int(N_local_arr[0])

        qlanc.X_q = np.zeros((qlanc.n_q, N_local, qlanc.n_bands), dtype=np.complex128)
        qlanc.Y_q = np.zeros((qlanc.n_q, N_local, qlanc.n_bands), dtype=np.complex128)
        qlanc.rho = np.zeros(N_local, dtype=np.float64)

        comm.Recv(qlanc.X_q, source=0, tag=1)
        comm.Recv(qlanc.Y_q, source=0, tag=2)
        comm.Recv(qlanc.rho, source=0, tag=3)

        # Set distributed state
        qlanc._distributed = True
        qlanc._N_global = metadata['N_global']
        qlanc._N_eff_global = metadata['N_eff_global']
        qlanc._N_local = N_local
        qlanc.N = N_local
        qlanc.N_eff = int(np.sum(qlanc.rho))

        # Build Julia symmetry cache
        qlanc.prepare_symmetrization(no_sym=not use_symmetries)
        qlanc.initialized = True

        return qlanc

QSpaceLanczos Class

tdscha.QSpaceLanczos.QSpaceLanczos(ensemble, lo_to_split=None, **kwargs)

Bases: Lanczos

Q-space Lanczos for spectral calculations exploiting Bloch's theorem.

This class works in the q-space mode basis to exploit momentum conservation, reducing the psi vector size by ~N_cell and the anharmonic computation by ~N_cell.

Only Wigner formalism is supported. Requires Julia extension.

Initialize the Q-Space Lanczos.

Parameters:

Name	Type	Description	Default
ensemble	Ensemble	The SSCHA ensemble.	required
lo_to_split	string, ndarray, or None	LO-TO splitting mode. If None (default), LO-TO splitting correction is neglected. If "random", a random direction is used. If an ndarray, it specifies the q-direction for the LO-TO splitting correction.	None
**kwargs		Additional keyword arguments passed to the parent Lanczos class.	{}

Source code in tdscha/QSpaceLanczos.py

def __init__(self, ensemble, lo_to_split=None, **kwargs):
    """Initialize the Q-Space Lanczos.

    Parameters
    ----------
    ensemble : sscha.Ensemble.Ensemble
        The SSCHA ensemble.
    lo_to_split : string, ndarray, or None
        LO-TO splitting mode. If None (default), LO-TO splitting correction
        is neglected. If "random", a random direction is used. If an ndarray,
        it specifies the q-direction for the LO-TO splitting correction.
    **kwargs
        Additional keyword arguments passed to the parent Lanczos class.
    """
    # Force Wigner and Julia mode
    kwargs['mode'] = DL.MODE_FAST_JULIA
    kwargs['use_wigner'] = True

    self.ensemble = ensemble
    super().__init__(ensemble, unwrap_symmetries=False, lo_to_split=lo_to_split, **kwargs)

    if not __JULIA_EXT__:
        raise ImportError(
            "QSpaceLanczos requires Julia. Install with: pip install julia"
        )
    self.use_wigner = True

    # -- Add the q-space attributes --
    qspace_attrs = [
        'q_points', 'n_q', 'n_bands', 'w_q', 'pols_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',
        # Distributed mode attributes
        '_distributed', '_N_global', '_N_eff_global', '_N_local',
    ]
    self.__total_attributes__.extend(qspace_attrs)

    # If ensemble is None, perform a bare initialization like the parent
    if ensemble is None:
        return

    # == 1. Get q-space eigenmodes ==
    ws_sc, pols_sc, w_q, pols_q = self.dyn.DiagonalizeSupercell(
        return_qmodes=True, lo_to_split=lo_to_split)

    self.q_points = np.array(self.dyn.q_tot)  # (n_q, 3)
    self.n_q = len(self.q_points)
    self.n_bands = 3 * self.uci_structure.N_atoms  # uniform band count
    self.w_q = w_q        # (n_bands, n_q) from DiagonalizeSupercell
    self.pols_q = pols_q  # (3*n_at, n_bands, n_q) complex eigenvectors

    # The masses needs to be restricted to the primitive cell only
    self.m = self.m[:self.n_bands]

    # 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()

    # Q-space specific state (set by build_q_pair_map)
    self.iq_pert = None
    self.q_pair_map = None
    self.unique_pairs = None
    self._psi_size = None
    self._block_offsets_a = None
    self._block_offsets_b = None
    self._block_sizes = None

    # Symmetry data for q-space
    self._qspace_sym_data = None
    self._qspace_sym_q_map = None
    self.n_syms_qspace = 0

    # Distributed mode attributes (for MPI configuration distribution)
    self._distributed = False
    self._N_global = self.N
    self._N_eff_global = self.N_eff
    self._N_local = self.N

build_q_pair_map(iq_pert)

Find all (iq1, iq2) pairs satisfying q1 + q2 = q_pert + G.

Parameters:

Name	Type	Description	Default
iq_pert	int	Index of the perturbation q-point.	required

Source code in tdscha/QSpaceLanczos.py

def build_q_pair_map(self, iq_pert):
    """Find all (iq1, iq2) pairs satisfying q1 + q2 = q_pert + G.

    Parameters
    ----------
    iq_pert : int
        Index of the perturbation q-point.
    """
    bg = self.uci_structure.get_reciprocal_vectors() / (2 * np.pi)
    q_pert = self.q_points[iq_pert]

    self.iq_pert = iq_pert
    self.q_pair_map = np.zeros(self.n_q, dtype=np.int32)

    for iq1 in range(self.n_q):
        q_target = q_pert - self.q_points[iq1]
        found = False
        for iq2 in range(self.n_q):
            if CC.Methods.get_min_dist_into_cell(bg, q_target, self.q_points[iq2]) < 1e-6:
                self.q_pair_map[iq1] = iq2
                found = True
                break
        if not found:
            raise ValueError(
                "Could not find partner for q1={} with q_pert={}".format(
                    self.q_points[iq1], q_pert))

    # Unique pairs: iq1 <= iq2 (avoids double-counting)
    self.unique_pairs = []
    for iq1 in range(self.n_q):
        iq2 = self.q_pair_map[iq1]
        if iq1 <= iq2:
            self.unique_pairs.append((iq1, iq2))

    # Pre-compute block layout
    self._compute_block_layout()

get_psi_size()

Return the total size of the psi vector.

Source code in tdscha/QSpaceLanczos.py

def get_psi_size(self):
    """Return the total size of the psi vector."""
    if self._psi_size is None:
        raise ValueError("Must call build_q_pair_map first")
    return self._psi_size

get_static_psi_size()

Return psi size for the static layout: [R, one W sector].

This equals the end of the a' sector, i.e. the start of b'.

Source code in tdscha/QSpaceLanczos.py

def get_static_psi_size(self):
    """Return psi size for the static layout: [R, one W sector].

    This equals the end of the a' sector, i.e. the start of b'.
    """
    return self._block_offsets_b[0]

get_block_offset(pair_idx, sector='a')

Get the offset into psi for a given pair index.

Parameters:

Name	Type	Description	Default
pair_idx	int	Index into self.unique_pairs.	required
sector	str	'a' for a' sector, 'b' for b' sector.	'a'

Source code in tdscha/QSpaceLanczos.py

def get_block_offset(self, pair_idx, sector='a'):
    """Get the offset into psi for a given pair index.

    Parameters
    ----------
    pair_idx : int
        Index into self.unique_pairs.
    sector : str
        'a' for a' sector, 'b' for b' sector.
    """
    if sector == 'a':
        return self._block_offsets_a[pair_idx]
    else:
        return self._block_offsets_b[pair_idx]

get_block_size(pair_idx)

Get the number of entries for this pair.

Source code in tdscha/QSpaceLanczos.py

def get_block_size(self, pair_idx):
    """Get the number of entries for this pair."""
    return self._block_sizes[pair_idx]

get_R1_q()

Extract R^(1) from psi (n_bands complex entries at q_pert).

Source code in tdscha/QSpaceLanczos.py

def get_R1_q(self):
    """Extract R^(1) from psi (n_bands complex entries at q_pert)."""
    return self.psi[:self.n_bands].copy()

get_block(pair_idx, sector='a', source=None)

Reconstruct full (n_bands, n_bands) matrix from psi storage.

Parameters:

Name	Type	Description	Default
pair_idx	int	Index into self.unique_pairs.	required
sector	str	'a' or 'b'.	'a'
source	ndarray or None	If provided, read from this array instead of self.psi.	None

Returns:

Type	Description
(ndarray(n_bands, n_bands), complex128)

Source code in tdscha/QSpaceLanczos.py

def get_block(self, pair_idx, sector='a', source=None):
    """Reconstruct full (n_bands, n_bands) matrix from psi storage.

    Parameters
    ----------
    pair_idx : int
        Index into self.unique_pairs.
    sector : str
        'a' or 'b'.
    source : ndarray or None
        If provided, read from this array instead of self.psi.

    Returns
    -------
    ndarray(n_bands, n_bands), complex128
    """
    nb = self.n_bands
    iq1, iq2 = self.unique_pairs[pair_idx]
    offset = self.get_block_offset(pair_idx, sector)
    size = self.get_block_size(pair_idx)

    data = self.psi if source is None else source
    raw = data[offset:offset + size]

    if iq1 < iq2:
        # Full block, row-major
        return raw.reshape(nb, nb).copy()
    else:
        # Upper triangle storage
        return self._unpack_upper_triangle(raw, nb)

get_a1_block(pair_idx)

Get the a'(1) block for pair_idx.

Source code in tdscha/QSpaceLanczos.py

def get_a1_block(self, pair_idx):
    """Get the a'(1) block for pair_idx."""
    return self.get_block(pair_idx, 'a')

get_b1_block(pair_idx)

Get the b'(1) block for pair_idx.

Source code in tdscha/QSpaceLanczos.py

def get_b1_block(self, pair_idx):
    """Get the b'(1) block for pair_idx."""
    return self.get_block(pair_idx, 'b')

set_block_in_psi(pair_idx, matrix, sector, target_psi)

Write a (n_bands, n_bands) block into the target psi vector.

Parameters:

Name	Type	Description	Default
pair_idx	int		required
matrix	ndarray(n_bands, n_bands)		required
sector	str(a or b)		required
target_psi	ndarray — the psi vector to write into		required

Source code in tdscha/QSpaceLanczos.py

def set_block_in_psi(self, pair_idx, matrix, sector, target_psi):
    """Write a (n_bands, n_bands) block into the target psi vector.

    Parameters
    ----------
    pair_idx : int
    matrix : ndarray(n_bands, n_bands)
    sector : str ('a' or 'b')
    target_psi : ndarray — the psi vector to write into
    """
    nb = self.n_bands
    iq1, iq2 = self.unique_pairs[pair_idx]
    offset = self.get_block_offset(pair_idx, sector)

    if iq1 < iq2:
        # Full block
        target_psi[offset:offset + nb * nb] = matrix.ravel()
    else:
        # Upper triangle
        target_psi[offset:offset + self.get_block_size(pair_idx)] = \
            self._pack_upper_triangle(matrix, nb)

mask_dot_wigner(debug=False)

Build the mask for Hermitian inner product with upper-triangle storage.

For full blocks (iq1 < iq2): factor 2 for the conjugate block (iq2, iq1). For diagonal blocks (iq1 == iq2): off-diagonal factor 2, diagonal factor 1.

Returns:

Type	Description
(ndarray(psi_size), float64)

Source code in tdscha/QSpaceLanczos.py

def mask_dot_wigner(self, debug=False):
    """Build the mask for Hermitian inner product with upper-triangle storage.

    For full blocks (iq1 < iq2): factor 2 for the conjugate block (iq2, iq1).
    For diagonal blocks (iq1 == iq2): off-diagonal factor 2, diagonal factor 1.

    Returns
    -------
    ndarray(psi_size,), float64
    """
    mask = np.ones(self.get_psi_size(), dtype=np.float64)

    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        offset_a = self.get_block_offset(pair_idx, 'a')
        offset_b = self.get_block_offset(pair_idx, 'b')
        size = self.get_block_size(pair_idx)

        if iq1 < iq2:
            # Full block: factor 2 for the conjugate block
            mask[offset_a:offset_a + size] = 2
            mask[offset_b:offset_b + size] = 2
        else:  # iq1 == iq2, upper triangle storage
            nb = self.n_bands
            idx_a = offset_a
            idx_b = offset_b
            for i in range(nb):
                # Diagonal element: factor 1
                mask[idx_a] = 1
                mask[idx_b] = 1
                idx_a += 1
                idx_b += 1
                # Off-diagonal elements: factor 2
                for j in range(i + 1, nb):
                    mask[idx_a] = 2
                    mask[idx_b] = 2
                    idx_a += 1
                    idx_b += 1

    return mask

apply_L1_FT(transpose=False)

Apply the harmonic part of L in q-space (Wigner formalism).

L_harm is block-diagonal: R sector: -(w_q_pert[nu])^2 * R[nu] a' sector: -(w1 - w2)^2 * a' b' sector: -(w1 + w2)^2 * b'

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_L1_FT(self, transpose=False):
    """Apply the harmonic part of L in q-space (Wigner formalism).

    L_harm is block-diagonal:
      R sector: -(w_q_pert[nu])^2 * R[nu]
      a' sector: -(w1 - w2)^2 * a'
      b' sector: -(w1 + w2)^2 * b'

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    out = np.zeros(self.get_psi_size(), dtype=np.complex128)

    if self.ignore_harmonic:
        return out

    # R sector — mask acoustic modes
    w_qp = self.w_q[:, self.iq_pert]  # (n_bands,)
    v_pert = self.valid_modes_q[:, self.iq_pert]
    out[:self.n_bands] = -(w_qp ** 2) * self.psi[:self.n_bands]
    out[:self.n_bands] *= v_pert  # zero out acoustic

    # a' and b' sectors — mask pairs involving acoustic modes
    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        w1 = self.w_q[:, iq1]  # (n_bands,)
        w2 = self.w_q[:, iq2]  # (n_bands,)
        v1 = self.valid_modes_q[:, iq1]
        v2 = self.valid_modes_q[:, iq2]
        valid_mask = np.outer(v1, v2)

        a_block = self.get_a1_block(pair_idx)
        b_block = self.get_b1_block(pair_idx)

        w_minus2 = np.subtract.outer(w1, w2) ** 2  # (n_bands, n_bands)
        w_plus2 = np.add.outer(w1, w2) ** 2

        self.set_block_in_psi(pair_idx, -w_minus2 * a_block * valid_mask, 'a', out)
        self.set_block_in_psi(pair_idx, -w_plus2 * b_block * valid_mask, 'b', out)

    return out

get_chi_minus_q()

Get chi^- for each unique pair as a list of (n_bands, n_bands) matrices.

chi^-_{nu1, nu2} = (w1 - w2)(n1 - n2) / (2 * w1 * w2) Entries involving acoustic modes (w < acoustic_eps) are set to 0.

Source code in tdscha/QSpaceLanczos.py

def get_chi_minus_q(self):
    """Get chi^- for each unique pair as a list of (n_bands, n_bands) matrices.

    chi^-_{nu1, nu2} = (w1 - w2)(n1 - n2) / (2 * w1 * w2)
    Entries involving acoustic modes (w < acoustic_eps) are set to 0.
    """
    chi_list = []
    for iq1, iq2 in self.unique_pairs:
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        # Outer mask: both bands must be non-acoustic
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        # Safe division: use np.where to avoid 0/0
        denom = 2.0 * w1_mat * w2_mat
        chi = np.where(valid_mask,
                       (w1_mat - w2_mat) * (n1_mat - n2_mat) / np.where(valid_mask, denom, 1.0),
                       0.0)
        chi_list.append(chi)
    return chi_list

get_chi_plus_q()

Get chi^+ for each unique pair as a list of (n_bands, n_bands) matrices.

chi^+_{nu1, nu2} = (w1 + w2)(1 + n1 + n2) / (2 * w1 * w2) Entries involving acoustic modes (w < acoustic_eps) are set to 0.

Source code in tdscha/QSpaceLanczos.py

def get_chi_plus_q(self):
    """Get chi^+ for each unique pair as a list of (n_bands, n_bands) matrices.

    chi^+_{nu1, nu2} = (w1 + w2)(1 + n1 + n2) / (2 * w1 * w2)
    Entries involving acoustic modes (w < acoustic_eps) are set to 0.
    """
    chi_list = []
    for iq1, iq2 in self.unique_pairs:
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        denom = 2.0 * w1_mat * w2_mat
        chi = np.where(valid_mask,
                       (w1_mat + w2_mat) * (1 + n1_mat + n2_mat) / np.where(valid_mask, denom, 1.0),
                       0.0)
        chi_list.append(chi)
    return chi_list

get_alpha1_beta1_wigner_q(get_alpha=True)

Get the perturbation on alpha (Upsilon) from the q-space psi.

Transforms a'/b' blocks back to the alpha1 perturbation that the Julia code needs.

alpha1[iq1, iq2] = (w1w2/X) * [sqrt(-0.5chi_minus)a' - sqrt(0.5chi_plus)*b']

Returns:

Type	Description
list of ndarray(n_bands, n_bands) — one per unique pair

Source code in tdscha/QSpaceLanczos.py

def get_alpha1_beta1_wigner_q(self, get_alpha=True):
    """Get the perturbation on alpha (Upsilon) from the q-space psi.

    Transforms a'/b' blocks back to the alpha1 perturbation that the
    Julia code needs.

    alpha1[iq1, iq2] = (w1*w2/X) * [sqrt(-0.5*chi_minus)*a' - sqrt(0.5*chi_plus)*b']

    Returns
    -------
    list of ndarray(n_bands, n_bands) — one per unique pair
    """
    chi_minus_list = self.get_chi_minus_q()
    chi_plus_list = self.get_chi_plus_q()

    alpha1_blocks = []
    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        w1 = self.w_q[:, iq1]
        w2 = self.w_q[:, iq2]
        n1, v1 = self._safe_bose_and_mask(iq1)
        n2, v2 = self._safe_bose_and_mask(iq2)
        valid_mask = np.outer(v1, v2)

        w1_mat = np.tile(w1, (self.n_bands, 1)).T
        w2_mat = np.tile(w2, (self.n_bands, 1))
        n1_mat = np.tile(n1, (self.n_bands, 1)).T
        n2_mat = np.tile(n2, (self.n_bands, 1))

        X = (1 + 2 * n1_mat) * (1 + 2 * n2_mat) / 8
        # Safe division for w2_on_X: when acoustic, X→∞ and w→0, set to 0
        w2_on_X = np.where(valid_mask,
                           (w1_mat * w2_mat) / np.where(valid_mask, X, 1.0),
                           0.0)
        chi_minus = chi_minus_list[pair_idx]
        chi_plus = chi_plus_list[pair_idx]

        a_block = self.get_a1_block(pair_idx)
        b_block = self.get_b1_block(pair_idx)

        if get_alpha:
            new_a = w2_on_X * np.sqrt(-0.5 * chi_minus) * a_block
            new_b = w2_on_X * np.sqrt(+0.5 * chi_plus) * b_block
            alpha1 = new_a - new_b
        else:
            X_safe = np.where(valid_mask, X, 1.0)
            new_a = np.where(valid_mask,
                             (np.sqrt(-0.5 * chi_minus) / X_safe) * a_block,
                             0.0)
            new_b = np.where(valid_mask,
                             (np.sqrt(+0.5 * chi_plus) / X_safe) * b_block,
                             0.0)
            alpha1 = new_a + new_b

        alpha1_blocks.append(alpha1)
    return alpha1_blocks

apply_anharmonic_FT(transpose=False, **kwargs)

Apply the anharmonic part of L in q-space (Wigner formalism).

Calls the Julia q-space extension to compute the perturbed averages, then assembles the output psi vector.

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_anharmonic_FT(self, transpose=False, **kwargs):
    """Apply the anharmonic part of L in q-space (Wigner formalism).

    Calls the Julia q-space extension to compute the perturbed averages,
    then assembles the output psi vector.

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    # If both D3 and D4 are ignored, return zero
    if self.ignore_v3 and self.ignore_v4:
        return np.zeros(self.get_psi_size(), dtype=np.complex128)

    import julia.Main

    R1 = self.get_R1_q()
    # If D3 is ignored, zero out R1 so that D3 weight is zero
    if self.ignore_v3:
        R1 = np.zeros_like(R1)
    alpha1_blocks = self.get_alpha1_beta1_wigner_q(get_alpha=True)
    alpha1_flat = self._flatten_blocks(alpha1_blocks)

    # Call Julia
    f_pert, d2v_blocks = self._call_julia_qspace(R1, alpha1_flat)

    # Build output psi
    final_psi = np.zeros(self.get_psi_size(), dtype=np.complex128)

    # R sector
    final_psi[:self.n_bands] = f_pert

    # a'/b' sectors
    chi_minus_list = self.get_chi_minus_q()
    chi_plus_list = self.get_chi_plus_q()

    for pair_idx, (iq1, iq2) in enumerate(self.unique_pairs):
        d2v_block = d2v_blocks[pair_idx]
        pert_a = np.sqrt(-0.5 * chi_minus_list[pair_idx]) * d2v_block
        pert_b = -np.sqrt(+0.5 * chi_plus_list[pair_idx]) * d2v_block
        self.set_block_in_psi(pair_idx, pert_a, 'a', final_psi)
        self.set_block_in_psi(pair_idx, pert_b, 'b', final_psi)

    return final_psi

apply_full_L(target=None, force_t_0=False, force_FT=True, transpose=False, fast_lanczos=True)

Apply the full L operator in q-space.

Parameters:

Name	Type	Description	Default
target	ndarray or None	If provided, copy into self.psi first.	None
transpose	bool	Not used for Hermitian Lanczos.	False

Returns:

Type	Description
(ndarray(psi_size), complex128)

Source code in tdscha/QSpaceLanczos.py

def apply_full_L(self, target=None, force_t_0=False, force_FT=True,
                 transpose=False, fast_lanczos=True):
    """Apply the full L operator in q-space.

    Parameters
    ----------
    target : ndarray or None
        If provided, copy into self.psi first.
    transpose : bool
        Not used for Hermitian Lanczos.

    Returns
    -------
    ndarray(psi_size,), complex128
    """
    if target is not None:
        self.psi = target.copy()

    result = self.apply_L1_FT(transpose=transpose)
    result += self.apply_anharmonic_FT(transpose=transpose)

    return result

run_FT(n_iter, save_dir=None, save_each=5, verbose=True, n_rep_orth=0, n_ortho=10, flush_output=True, debug=False, prefix='LANCZOS', run_simm=None, optimized=False, reorthogonalize=True)

Run the Hermitian Lanczos algorithm for q-space.

This is the same structure as the parent run_FT but with: 1. Forced run_simm = True (Hermitian) 2. Hermitian dot products: psi.conj().dot(psi * mask).real 3. Complex128 psi 4. Real coefficients (guaranteed by Hermitian L)

Source code in tdscha/QSpaceLanczos.py

    def run_FT(self, n_iter, save_dir=None, save_each=5, verbose=True,
               n_rep_orth=0, n_ortho=10, flush_output=True, debug=False,
               prefix="LANCZOS", run_simm=None, optimized=False,
               reorthogonalize=True):
        """Run the Hermitian Lanczos algorithm for q-space.

        This is the same structure as the parent run_FT but with:
        1. Forced run_simm = True (Hermitian)
        2. Hermitian dot products: psi.conj().dot(psi * mask).real
        3. Complex128 psi
        4. Real coefficients (guaranteed by Hermitian L)
        """
        self.verbose = verbose

        if not self.initialized:
            if verbose:
                print('Not initialized. Now we symmetrize\n')
            self.prepare_symmetrization()

        ERROR_MSG = """
Error, you must initialize a perturbation to start the Lanczos.
Use prepare_mode_q or prepare_perturbation_q before calling run_FT.
"""
        if self.psi is None:
            raise ValueError(ERROR_MSG)

        mask_dot = self.mask_dot_wigner(debug)
        psi_norm = np.real(np.conj(self.psi).dot(self.psi * mask_dot))
        if np.isnan(psi_norm) or psi_norm == 0:
            raise ValueError(ERROR_MSG)

        # Force symmetric Lanczos
        run_simm = True

        if Parallel.am_i_the_master():
            if save_dir is not None:
                if not os.path.exists(save_dir):
                    os.makedirs(save_dir)

        if verbose:
            print('Running the Hermitian Lanczos in q-space')
            print()

        # Get current step
        i_step = len(self.a_coeffs)

        if verbose:
            header = """
<=====================================>
|                                     |
|     Q-SPACE LANCZOS ALGORITHM       |
|                                     |
<=====================================>

Starting the algorithm.
Starting from step %d
""" % i_step
            print(header)

        # Initialize
        if i_step == 0:
            self.basis_Q = []
            self.basis_P = []
            self.s_norm = []
            norm = np.sqrt(np.real(np.conj(self.psi).dot(self.psi * mask_dot)))
            first_vector = self.psi / norm
            self.basis_Q.append(first_vector)
            self.basis_P.append(first_vector)
            self.s_norm.append(1)
        else:
            if verbose:
                print('Restarting the Lanczos')
            self.basis_Q = list(self.basis_Q)
            self.basis_P = list(self.basis_P)
            self.s_norm = list(self.s_norm)
            self.a_coeffs = list(self.a_coeffs)
            self.b_coeffs = list(self.b_coeffs)
            self.c_coeffs = list(self.c_coeffs)

        psi_q = self.basis_Q[-1]
        psi_p = self.basis_P[-1]

        next_converged = False
        converged = False

        for i in range(i_step, i_step + n_iter):
            if verbose:
                print("\n ===== NEW STEP %d =====\n" % (i + 1))
                if flush_output:
                    sys.stdout.flush()

            # Apply L (Hermitian => p_L = L_q)
            t1 = time.time()
            L_q = self.apply_full_L(psi_q)
            p_L = np.copy(L_q)
            t2 = time.time()

            # p normalization
            c_old = 1
            if len(self.c_coeffs) > 0:
                c_old = self.c_coeffs[-1]
            p_norm = self.s_norm[-1] / c_old

            # a coefficient (real for Hermitian L)
            a_coeff = np.real(np.conj(psi_p).dot(L_q * mask_dot)) * p_norm

            if np.isnan(a_coeff):
                raise ValueError("Invalid value in Lanczos. Check frequencies/initialization.")

            # Residuals
            rk = L_q - a_coeff * psi_q
            if len(self.basis_Q) > 1:
                rk -= self.c_coeffs[-1] * self.basis_Q[-2]

            sk = p_L - a_coeff * psi_p
            if len(self.basis_P) > 1:
                old_p_norm = self.s_norm[-2]
                if len(self.c_coeffs) >= 2:
                    old_p_norm = self.s_norm[-2] / self.c_coeffs[-2]
                sk -= self.b_coeffs[-1] * self.basis_P[-2] * (old_p_norm / p_norm)

            # s_norm
            s_norm = np.sqrt(np.real(np.conj(sk).dot(sk * mask_dot)))
            sk_tilde = sk / s_norm
            s_norm *= p_norm

            # b and c coefficients (real, should be equal for Hermitian L)
            b_coeff = np.sqrt(np.real(np.conj(rk).dot(rk * mask_dot)))
            c_coeff = np.real(np.conj(sk_tilde).dot((rk / b_coeff) * mask_dot)) * s_norm

            self.a_coeffs.append(a_coeff)

            if np.abs(b_coeff) < __EPSILON__ or next_converged:
                if verbose:
                    print("Converged (b = {})".format(b_coeff))
                converged = True
                break
            if np.abs(c_coeff) < __EPSILON__:
                if verbose:
                    print("Converged (c = {})".format(c_coeff))
                converged = True
                break

            psi_q = rk / b_coeff
            psi_p = sk_tilde.copy()

            # Gram-Schmidt reorthogonalization
            if reorthogonalize:
                # Correct Hermitian MGS: orthogonalize against all Q vectors
                # (which are unit-normalized in the mask inner product)
                new_q = psi_q.copy()

                for j in range(len(self.basis_Q)):
                    coeff = np.real(np.conj(self.basis_Q[j]).dot(new_q * mask_dot))
                    new_q -= coeff * self.basis_Q[j]

                normq = np.sqrt(np.real(np.conj(new_q).dot(new_q * mask_dot)))
                if normq < __EPSILON__:
                    next_converged = True
                new_q /= normq

                # Hermitian L: P = Q, s_norm = c_coeff (since <P,Q> = <Q,Q> = 1)
                new_p = new_q.copy()
                s_norm = c_coeff
            else:
                # Existing biconjugate GS (for backward compatibility)
                new_q = psi_q.copy()
                new_p = psi_p.copy()

                for k_orth in range(n_rep_orth):
                    start = max(0, len(self.basis_P) - (n_ortho or len(self.basis_P)))

                    for j in range(start, len(self.basis_P)):
                        coeff1 = np.real(np.conj(self.basis_P[j]).dot(new_q * mask_dot))
                        coeff2 = np.real(np.conj(self.basis_Q[j]).dot(new_p * mask_dot))
                        new_q -= coeff1 * self.basis_P[j]
                        new_p -= coeff2 * self.basis_Q[j]

                    normq = np.sqrt(np.real(np.conj(new_q).dot(new_q * mask_dot)))
                    if normq < __EPSILON__:
                        next_converged = True
                    new_q /= normq

                    normp = np.real(np.conj(new_p).dot(new_p * mask_dot))
                    if np.abs(normp) < __EPSILON__:
                        next_converged = True
                    new_p /= normp

                    s_norm = c_coeff / np.real(np.conj(new_p).dot(new_q * mask_dot))

            if not converged:
                self.basis_Q.append(new_q)
                self.basis_P.append(new_p)
                psi_q = new_q.copy()
                psi_p = new_p.copy()

                self.b_coeffs.append(b_coeff)
                self.c_coeffs.append(c_coeff)
                self.s_norm.append(s_norm)

            if verbose:
                print("Time for L application: %d s" % (t2 - t1))
                print("a_%d = %.8e" % (i, self.a_coeffs[-1]))
                print("b_%d = %.8e" % (i, self.b_coeffs[-1]))
                print("c_%d = %.8e" % (i, self.c_coeffs[-1]))
                print("|b-c| = %.8e" % np.abs(self.b_coeffs[-1] - self.c_coeffs[-1]))

            if save_dir is not None:
                if (i + 1) % save_each == 0:
                    self.save_status("%s/%s_STEP%d" % (save_dir, prefix, i + 1))

            if verbose:
                print("Lanczos step %d completed." % (i + 1))

        if converged and verbose:
            print("   last a coeff = {}".format(self.a_coeffs[-1]))

prepare_mode_q(iq, band_index)

Prepare perturbation for mode (q, nu).

Parameters:

Name	Type	Description	Default
iq	int	Index of the q-point.	required
band_index	int	Band index (0-based).	required

Source code in tdscha/QSpaceLanczos.py

def prepare_mode_q(self, iq, band_index):
    """Prepare perturbation for mode (q, nu).

    Parameters
    ----------
    iq : int
        Index of the q-point.
    band_index : int
        Band index (0-based).
    """
    if band_index < 0 or band_index >= self.n_bands:
        raise ValueError("Invalid band index for perturbation: {}".format(band_index))

    self.build_q_pair_map(iq)
    self.reset_q()
    self.psi[band_index] = 1.0 + 0j
    self.perturbation_modulus = 1.0

prepare_ir(effective_charges=None, pol_vec=np.array([1.0, 0.0, 0.0]))

PREPARE LANCZOS FOR INFRARED SPECTRUM COMPUTATION

In this subroutine we prepare the lanczos algorithm for the computation of the infrared spectrum signal.

Source code in tdscha/QSpaceLanczos.py

def prepare_ir(self, effective_charges = None, pol_vec = np.array([1.0, 0.0, 0.0])):
    """
    PREPARE LANCZOS FOR INFRARED SPECTRUM COMPUTATION
    =================================================

    In this subroutine we prepare the lanczos algorithm for the computation of the
    infrared spectrum signal.

    Parameters
    ----------
        effective_charges : ndarray(size = (n_atoms, 3, 3), dtype = np.double)
            The effective charges. Indices are: Number of atoms in the unit cell,
            electric field component, atomic coordinate. If None, the effective charges
            contained in the dynamical matrix will be considered.
        pol_vec : ndarray(size = 3)
            The polarization vector of the light.
    """

    ec = self.dyn.effective_charges
    if not effective_charges is None:
        ec = effective_charges

    assert not ec is None, "Error, no effective charge found. Cannot initialize IR responce"

    z_eff = np.einsum("abc, b", ec, pol_vec)

    # Get the gamma effective charge
    # FIX: remove double mass scaling and add supercell factor
    n_cell = np.prod(self.dyn.GetSupercell())
    new_zeff = z_eff.ravel() * np.sqrt(n_cell)

    # This is a Gamma perturbation
    self.prepare_perturbation_q(0, new_zeff)

prepare_raman(pol_vec_in=np.array([1.0, 0.0, 0.0]), pol_vec_out=np.array([1.0, 0.0, 0.0]), mixed=False, pol_in_2=None, pol_out_2=None, unpolarized=None)

PREPARE LANCZOS FOR RAMAN SPECTRUM COMPUTATION

In this subroutine we prepare the lanczos algorithm for the computation of the Raman spectrum signal.

Source code in tdscha/QSpaceLanczos.py

def prepare_raman(self, pol_vec_in=np.array([1.0, 0.0, 0.0]), pol_vec_out=np.array([1.0, 0.0, 0.0]), 
                 mixed=False, pol_in_2=None, pol_out_2=None, unpolarized=None):
    """
    PREPARE LANCZOS FOR RAMAN SPECTRUM COMPUTATION
    ==============================================

    In this subroutine we prepare the lanczos algorithm for the computation of the
    Raman spectrum signal.

    Parameters
    ----------
        pol_vec_in : ndarray(size = 3)
            The polarization vector of the incoming light
        pol_vec_out : ndarray(size = 3)
            The polarization vector for the outcoming light
        mixed : bool
            If True, add another component of the Raman tensor
        pol_in_2 : ndarray(size = 3) or None
            Second incoming polarization if mixed=True
        pol_out_2 : ndarray(size = 3) or None  
            Second outcoming polarization if mixed=True
        unpolarized : int or None
            The perturbation for unpolarized raman (if different from None, overrides the behaviour
            of pol_vec_in and pol_vec_out). Indices goes from 0 to 6 (included).
            0 is alpha^2
            1 + 2 + 3 + 4 + 5 + 6 are beta^2
            alpha_0 = (xx + yy + zz)^2/9
            beta_1 = (xx -yy)^2 / 2
            beta_2 = (xx -zz)^2 / 2
            beta_3 = (yy -zz)^2 / 2
            beta_4 = 3xy^2
            beta_5 = 3xz^2
            beta_6 = 3yz^2

            The total unpolarized raman intensity is 45 alpha^2 + 7 beta^2
    """
    # Check if the raman tensor is present
    assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize Raman response"

    n_cell = np.prod(self.dyn.GetSupercell())

    if unpolarized is None:
        # Get the raman vector (apply the ASR and contract the raman tensor with the polarization vectors)
        raman_v = self.dyn.GetRamanVector(pol_vec_in, pol_vec_out)

        if mixed:
            print('Prepare Raman')
            print('Adding other component of the Raman tensor')
            raman_v += self.dyn.GetRamanVector(pol_in_2, pol_out_2)

        # Scale for Γ-point constant perturbation
        new_raman_v = raman_v.ravel() * np.sqrt(n_cell)

        # Convert in the polarization basis
        self.prepare_perturbation_q(0, new_raman_v)
    else:
        px = np.array([1, 0, 0])
        py = np.array([0, 1, 0])
        pz = np.array([0, 0, 1])

        if unpolarized == 0:
            # Alpha = (xx + yy + zz)^2/9
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v, add=True)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / 3
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 1:
            # (xx - yy)^2 / 2
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 2:
            # (xx - zz)^2 / 2
            raman_v = self.dyn.GetRamanVector(px, px)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 3:
            # (yy - zz)^2 / 2
            raman_v = self.dyn.GetRamanVector(py, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v)

            raman_v = self.dyn.GetRamanVector(pz, pz)
            new_raman_v = -raman_v.ravel() * np.sqrt(n_cell) / np.sqrt(2)
            self.prepare_perturbation_q(0, new_raman_v, add=True)
        elif unpolarized == 4:
            # 3 xy^2
            raman_v = self.dyn.GetRamanVector(px, py)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        elif unpolarized == 5:
            # 3 yz^2
            raman_v = self.dyn.GetRamanVector(py, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        elif unpolarized == 6:
            # 3 xz^2
            raman_v = self.dyn.GetRamanVector(px, pz)
            new_raman_v = raman_v.ravel() * np.sqrt(n_cell) * np.sqrt(3)
            self.prepare_perturbation_q(0, new_raman_v)
        else:
            raise ValueError(f"Error, unpolarized must be between [0, ..., 6] got invalid {unpolarized}.")

prepare_unpolarized_raman(index=0, debug=False)

PREPARE UNPOLARIZED RAMAN SIGNAL

The raman tensor is read from the dynamical matrix provided by the original ensemble.

The perturbations are prepared according to the formula (see https://doi.org/10.1021/jp5125266)

..math:

I_unpol = 45/9 (xx + yy + zz)^2
          + 7/2 [(xx-yy)^2 + (xx-zz)^2 + (yy-zz)^2]
          + 7 * 3 [(xy)^2 + (yz)^2 + (xz)^2]

Note: This method prepares the raw components WITHOUT prefactors. Use get_prefactors_unpolarized_raman() to get the correct prefactors.

Source code in tdscha/QSpaceLanczos.py

def prepare_unpolarized_raman(self, index=0, debug=False):
    """
    PREPARE UNPOLARIZED RAMAN SIGNAL
    ================================

    The raman tensor is read from the dynamical matrix provided by the original ensemble.

    The perturbations are prepared according to the formula (see https://doi.org/10.1021/jp5125266)

    ..math:

        I_unpol = 45/9 (xx + yy + zz)^2
                  + 7/2 [(xx-yy)^2 + (xx-zz)^2 + (yy-zz)^2]
                  + 7 * 3 [(xy)^2 + (yz)^2 + (xz)^2]

    Note: This method prepares the raw components WITHOUT prefactors.
    Use get_prefactors_unpolarized_raman() to get the correct prefactors.
    """
    # Check if the raman tensor is present
    assert not self.dyn.raman_tensor is None, "Error, no Raman tensor found. Cannot initialize the Raman response"

    labels = [i for i in range(7)]
    if index not in labels:
        raise ValueError(f'{index} should be in {labels}')

    epols = {'x': np.array([1, 0, 0]),
             'y': np.array([0, 1, 0]),
             'z': np.array([0, 0, 1])}

    n_cell = np.prod(self.dyn.GetSupercell())

    # (xx + yy + zz)^2
    if index == 0:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v += self.dyn.GetRamanVector(epols['y'], epols['y'])
        raman_v += self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (xx - yy)^2    
    elif index == 1:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v -= self.dyn.GetRamanVector(epols['y'], epols['y'])
    # (xx - zz)^2       
    elif index == 2:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['x'])
        raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (yy - zz)^2   
    elif index == 3:
        raman_v = self.dyn.GetRamanVector(epols['y'], epols['y'])
        raman_v -= self.dyn.GetRamanVector(epols['z'], epols['z'])
    # (xy)^2
    elif index == 4:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['y'])
    # (xz)^2
    elif index == 5:
        raman_v = self.dyn.GetRamanVector(epols['x'], epols['z'])
    # (yz)^2
    elif index == 6:
        raman_v = self.dyn.GetRamanVector(epols['y'], epols['z'])

    if debug:
        np.save(f'raman_v_{index}', raman_v)

    # Scale for Γ-point constant perturbation
    new_raman_v = raman_v.ravel() * np.sqrt(n_cell)

    # Convert in the polarization basis
    self.prepare_perturbation_q(0, new_raman_v)

    if debug:
        print(f'[NEW] Perturbation modulus with eq Raman tensors = {self.perturbation_modulus}')
    print()

prepare_perturbation_q(iq, vector, add=False)

Prepare perturbation at q from a real-space vector (3*n_at_uc,).

Projects the vector onto q-space eigenmodes at iq.

Parameters:

Name	Type	Description	Default
iq	int	Index of the q-point.	required
vector	ndarray(3 * n_at_uc)	Perturbation vector in Cartesian real space.	required
add	bool	If true, the perturbation is added on top of the one already setup. Calling add does not cause a reset of the Lanczos.	False

Source code in tdscha/QSpaceLanczos.py

def prepare_perturbation_q(self, iq, vector, add=False):
    """Prepare perturbation at q from a real-space vector (3*n_at_uc,).

    Projects the vector onto q-space eigenmodes at iq.

    Parameters
    ----------
    iq : int
        Index of the q-point.
    vector : ndarray(3*n_at_uc,)
        Perturbation vector in Cartesian real space.
    add : bool
        If true, the perturbation is added on top of the one already setup.
        Calling add does not cause a reset of the Lanczos.
    """
    if not add:
        self.build_q_pair_map(iq)
        self.reset_q()

    m = np.tile(self.uci_structure.get_masses_array(), (3, 1)).T.ravel()
    v_scaled = vector / np.sqrt(m)
    R1 = np.conj(self.pols_q[:, :, iq]).T @ v_scaled  # (n_bands,) complex
    self.psi[:self.n_bands] += R1
    self.perturbation_modulus = np.real(np.conj(R1) @ R1)

reset_q()

Reset the Lanczos state for q-space.

Source code in tdscha/QSpaceLanczos.py

def reset_q(self):
    """Reset the Lanczos state for q-space."""
    n = self.get_psi_size()
    self.psi = np.zeros(n, dtype=np.complex128)

    self.eigvals = None
    self.eigvects = None

    self.a_coeffs = []
    self.b_coeffs = []
    self.c_coeffs = []
    self.krilov_basis = []
    self.basis_P = []
    self.basis_Q = []
    self.s_norm = []

prepare_symmetrization(no_sym=False, verbose=True, symmetries=None)

Build q-space symmetry matrices and cache them in Julia.

Overrides the parent to build sparse complex symmetry matrices in the q-space mode basis.

Uses spglib on the unit cell (not supercell) to get correct fractional-coordinate rotations and translations, then converts to Cartesian for the representation matrices.

Source code in tdscha/QSpaceLanczos.py

def prepare_symmetrization(self, no_sym=False, verbose=True, symmetries=None):
    """Build q-space symmetry matrices and cache them in Julia.

    Overrides the parent to build sparse complex symmetry matrices
    in the q-space mode basis.

    Uses spglib on the unit cell (not supercell) to get correct
    fractional-coordinate rotations and translations, then converts
    to Cartesian for the representation matrices.
    """
    self.initialized = True

    if no_sym:
        # Identity only
        self.n_syms_qspace = 1
        n_total = self.n_q * self.n_bands
        # Build identity sparse matrix
        import julia.Main
        julia.Main.eval("""
        function init_identity_qspace(n_total::Int64)
            I_sparse = SparseArrays.sparse(
                Int32.(1:n_total), Int32.(1:n_total),
                ComplexF64.(ones(n_total)), n_total, n_total)
            _cached_qspace_symmetries[] = [I_sparse]
            return nothing
        end
        """)
        julia.Main.init_identity_qspace(int(n_total))
        return

    if not __SPGLIB__:
        raise ImportError("spglib required for symmetrization")

    # Get symmetries from the UNIT CELL directly via spglib.
    # spglib returns rotations and translations in fractional
    # (crystal) coordinates. We convert rotations to Cartesian via
    # R_cart = M @ R_frac @ M^{-1} where M = unit_cell.T.
    spg_data = spglib.get_symmetry(self.uci_structure.get_spglib_cell())
    rot_frac_all = spg_data['rotations']   # (n_sym, 3, 3) integer
    trans_frac_all = spg_data['translations']  # (n_sym, 3) fractional

    M = self.uci_structure.unit_cell.T       # columns = lattice vectors
    Minv = np.linalg.inv(M)

    # Extract unique point-group rotations (keep first occurrence)
    unique_pg = {}
    for i in range(len(rot_frac_all)):
        key = rot_frac_all[i].tobytes()
        if key not in unique_pg:
            unique_pg[key] = i
    pg_indices = list(unique_pg.values())

    if verbose:
        print("Q-space: {} PG symmetries from {} total unit cell symmetries".format(
            len(pg_indices), len(rot_frac_all)))

    self._build_qspace_symmetries(
        rot_frac_all, trans_frac_all, pg_indices, M, Minv,
        verbose=verbose)

init(use_symmetries=True)

Initialize the q-space Lanczos calculation.

Source code in tdscha/QSpaceLanczos.py

def init(self, use_symmetries=True):
    """Initialize the q-space Lanczos calculation."""
    self.prepare_symmetrization(no_sym=not use_symmetries)
    self.initialized = True

Key Differences from Lanczos

Feature	Lanczos (real-space)	QSpaceLanczos (q-space)
Psi vector	Real (float64)	Complex (complex128)
Inner product	Standard dot product	Hermitian: $\langle p
Lanczos type	Bi-conjugate or symmetric	Hermitian symmetric (\(b = c\), real coefficients)
Two-phonon pairs	All supercell mode pairs	\((q_1, q_2)\) pairs with \(q_1 + q_2 = q_\text{pert}\)
Symmetries	Full space group	Point group only (translations analytic)
Backend	C / MPI / Julia	Julia only

Utility Functions

tdscha.QSpaceLanczos.find_q_index(q_target, q_points, bg, tol=1e-06)

Find the index of q_target in q_points up to a reciprocal lattice vector.

Parameters:

Name	Type	Description	Default
q_target	ndarray(3)	The q-point to find (Cartesian coordinates).	required
q_points	ndarray(n_q, 3)	Array of q-points.	required
bg	ndarray(3, 3)	Reciprocal lattice vectors / (2*pi), rows are vectors.	required
tol	float	Tolerance for matching.	1e-06

Returns:

Type	Description
int	Index of the matching q-point.

Source code in tdscha/QSpaceLanczos.py

def find_q_index(q_target, q_points, bg, tol=1e-6):
    """Find the index of q_target in q_points up to a reciprocal lattice vector.

    Parameters
    ----------
    q_target : ndarray(3,)
        The q-point to find (Cartesian coordinates).
    q_points : ndarray(n_q, 3)
        Array of q-points.
    bg : ndarray(3, 3)
        Reciprocal lattice vectors / (2*pi), rows are vectors.
    tol : float
        Tolerance for matching.

    Returns
    -------
    int
        Index of the matching q-point.
    """
    for iq, q in enumerate(q_points):
        dist = CC.Methods.get_min_dist_into_cell(bg, q_target, q)
        if dist < tol:
            return iq
    raise ValueError("Could not find q-point {} in the q-point list".format(q_target))