Skip to content

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 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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',
    ]
    self.__total_attributes__.extend(qspace_attrs)

    # == 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

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 
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
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 
344
345
346
347
348
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 
350
351
352
353
354
355
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 
357
358
359
360
361
362
363
364
365
366
367
368
369
370
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 
372
373
374
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 
376
377
378
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 
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
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 
439
440
441
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 
443
444
445
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 
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
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 
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
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 
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
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 
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
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 
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
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 
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
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 
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
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)

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 
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
    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):
        """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
            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 
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
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).
    """
    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 
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
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 
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
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 
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
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 
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
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 
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
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 
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
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(np.int64(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 
1498
1499
1500
1501
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 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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))

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 
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
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',
    ]
    self.__total_attributes__.extend(qspace_attrs)

    # == 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

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 
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
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 
344
345
346
347
348
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 
350
351
352
353
354
355
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 
357
358
359
360
361
362
363
364
365
366
367
368
369
370
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 
372
373
374
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 
376
377
378
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 
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
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 
439
440
441
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 
443
444
445
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 
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
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 
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
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 
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
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 
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
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 
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
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 
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
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 
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
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 
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
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)

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 
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
    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):
        """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
            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 
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
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).
    """
    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 
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
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 
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
1189
1190
1191
1192
1193
1194
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 
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208
1209
1210
1211
1212
1213
1214
1215
1216
1217
1218
1219
1220
1221
1222
1223
1224
1225
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
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 
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288
1289
1290
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 
1292
1293
1294
1295
1296
1297
1298
1299
1300
1301
1302
1303
1304
1305
1306
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 
1311
1312
1313
1314
1315
1316
1317
1318
1319
1320
1321
1322
1323
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
1358
1359
1360
1361
1362
1363
1364
1365
1366
1367
1368
1369
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(np.int64(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 
1498
1499
1500
1501
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 
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
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))