QSpaceHessian API Reference

tdscha.QSpaceHessian

Q-Space Free Energy Hessian

Computes the free energy Hessian d²F/dRdR in q-space (Bloch basis).

Instead of building the full D4 matrix explicitly (O(N⁴) memory, O(N⁶) time), solves L_static(q) x = e_i for each band at each irreducible q-point, where L_static is the static Liouvillian operator. The Hessian is then H(q) = inv(G(q)), where G(q) is the static susceptibility.

The static L operator is DIFFERENT from the spectral L used in Lanczos: - Static: R sector = +w², W sector = 1/Lambda (one 2-phonon sector) - Spectral: R sector = -w², a'/b' sectors = -(w1∓w2)² (two sectors)

Both share the same anharmonic core (ensemble averages of D3/D4).

The q-space block-diagonal structure gives a speedup of ~N_cell³ / N_q_irr over the real-space approach.

References: Monacelli & Mauri 2021 (Phys. Rev. B)

QSpaceHessian(ensemble, verbose=True, ignore_v3=False, ignore_v4=False, **kwargs)

Compute the free energy Hessian in q-space via iterative linear solves.

Uses the static Liouvillian operator L_static, which has the structure: - R sector: w² * R (positive, unlike spectral -w²) - W sector: (1/Lambda) * W (static 2-phonon propagator) - Anharmonic coupling via ensemble averages (same Julia core)

Parameters:

Name	Type	Description	Default
ensemble	Ensemble	The SSCHA ensemble.	required
verbose	bool	If True, print progress information.	True
**kwargs		Additional keyword arguments passed to QSpaceLanczos.	{}

Source code in tdscha/QSpaceHessian.py

def __init__(self, ensemble, verbose=True, ignore_v3=False, ignore_v4=False, **kwargs):
    from tdscha.QSpaceLanczos import QSpaceLanczos

    self.qlanc = QSpaceLanczos(ensemble, **kwargs)
    self.ensemble = ensemble
    self.verbose = verbose

    # Store flags locally AND set on QSpaceLanczos
    self.ignore_v3 = ignore_v3
    self.ignore_v4 = ignore_v4
    self.qlanc.ignore_v3 = ignore_v3
    self.qlanc.ignore_v4 = ignore_v4

    # Shortcuts
    self.n_q = self.qlanc.n_q
    self.n_bands = self.qlanc.n_bands
    self.q_points = self.qlanc.q_points
    self.w_q = self.qlanc.w_q
    self.pols_q = self.qlanc.pols_q

    # Results storage: iq -> H_q matrix (n_bands x n_bands)
    self.H_q_dict = {}

    # Irreducible q-point data (set by init)
    self.irr_qpoints = None
    self.q_star_map = None

    # Cached static quantities (set by _precompute_static_quantities)
    self._cached_w_qp_sq = None
    self._cached_inv_w_qp_sq = None
    self._cached_inv_lambda = None
    self._cached_lambda = None
    self._cached_yw_outer = None

    # Cached spglib symmetry data (set by _find_irreducible_qpoints)
    self._has_symmetry_data = False

init(use_symmetries=True)

Initialize the Lanczos engine and find irreducible q-points.

Parameters:

Name	Type	Description	Default
use_symmetries	bool	If True, use symmetries to reduce q-points.	True

Source code in tdscha/QSpaceHessian.py

def init(self, use_symmetries=True):
    """Initialize the Lanczos engine and find irreducible q-points.

    Parameters
    ----------
    use_symmetries : bool
        If True, use symmetries to reduce q-points.
    """
    self.qlanc.init(use_symmetries=use_symmetries)
    if use_symmetries and __SPGLIB__:
        self._find_irreducible_qpoints()
    else:
        self.irr_qpoints = list(range(self.n_q))
        self.q_star_map = {iq: [(iq, np.eye(3), np.zeros(3))]
                           for iq in range(self.n_q)}

compute_hessian_at_q(iq, tol=1e-06, max_iters=500, use_preconditioner=True, dense_fallback=False, use_mode_symmetry=True)

Compute the free energy Hessian at a single q-point.

Solves L_static(q) x_i = e_i for each non-acoustic band. G_q[j,i] = x_i[j] (R-sector), H_q = inv(G_q).

When use_mode_symmetry=True and degenerate modes are present, exploits Schur's lemma: L_static commutes with the little group of q, so G_q restricted to a d-dimensional irrep block is c*I_d. Only one solve per degenerate block is needed instead of d solves.

Parameters:

Name	Type	Description	Default
iq	int	Q-point index.	required
tol	float	Convergence tolerance for the iterative solver.	1e-06
max_iters	int	Maximum number of iterations.	500
use_preconditioner	bool	If True, use harmonic preconditioner.	True
dense_fallback	bool	If True, fall back to dense solve when iterative solvers fail. WARNING: this builds a psi_size x psi_size dense matrix, which can be very large for big supercells. Default is False.	False
use_mode_symmetry	bool	If True, exploit mode degeneracy to reduce the number of GMRES solves. Within each degenerate block, only one solve is performed and G_q is filled using Schur's lemma (G_block = c * I).	True

Returns:

Name	Type	Description
H_q	(ndarray(n_bands, n_bands), complex128)	The Hessian matrix in the mode basis at q-point iq.

Source code in tdscha/QSpaceHessian.py

def compute_hessian_at_q(self, iq, tol=1e-6, max_iters=500,
                         use_preconditioner=True,
                         dense_fallback=False,
                         use_mode_symmetry=True):
    """Compute the free energy Hessian at a single q-point.

    Solves L_static(q) x_i = e_i for each non-acoustic band.
    G_q[j,i] = x_i[j] (R-sector), H_q = inv(G_q).

    When use_mode_symmetry=True and degenerate modes are present,
    exploits Schur's lemma: L_static commutes with the little group
    of q, so G_q restricted to a d-dimensional irrep block is c*I_d.
    Only one solve per degenerate block is needed instead of d solves.

    Parameters
    ----------
    iq : int
        Q-point index.
    tol : float
        Convergence tolerance for the iterative solver.
    max_iters : int
        Maximum number of iterations.
    use_preconditioner : bool
        If True, use harmonic preconditioner.
    dense_fallback : bool
        If True, fall back to dense solve when iterative solvers fail.
        WARNING: this builds a psi_size x psi_size dense matrix, which
        can be very large for big supercells. Default is False.
    use_mode_symmetry : bool
        If True, exploit mode degeneracy to reduce the number of GMRES
        solves. Within each degenerate block, only one solve is performed
        and G_q is filled using Schur's lemma (G_block = c * I).

    Returns
    -------
    H_q : ndarray(n_bands, n_bands), complex128
        The Hessian matrix in the mode basis at q-point iq.
    """
    nb = self.n_bands

    # 1. Setup pair map and block layout
    self.qlanc.build_q_pair_map(iq)
    psi_size = self._get_static_psi_size()

    # 2. Precompute all static quantities (Lambda, Y_w, w²) once
    self._precompute_static_quantities()

    mask = self._static_mask()

    # Similarity transform arrays
    sqrt_mask = np.sqrt(mask)
    inv_sqrt_mask = np.zeros_like(sqrt_mask)
    nonzero = mask > 0
    inv_sqrt_mask[nonzero] = 1.0 / sqrt_mask[nonzero]

    # 3. Define transformed operator: L_tilde = D^{1/2} L_static D^{-1/2}
    def apply_L_tilde(x_tilde):
        x = x_tilde * inv_sqrt_mask
        Lx = self._apply_L_static_q(x)
        return Lx * sqrt_mask

    L_op = scipy.sparse.linalg.LinearOperator(
        (psi_size, psi_size), matvec=apply_L_tilde, dtype=np.complex128)

    # 4. Preconditioner
    M_op = None
    if use_preconditioner:
        def apply_M_tilde(x_tilde):
            return self._apply_harmonic_preconditioner(
                x_tilde, sqrt_mask, inv_sqrt_mask)
        M_op = scipy.sparse.linalg.LinearOperator(
            (psi_size, psi_size), matvec=apply_M_tilde,
            dtype=np.complex128)

    # 5. Identify non-acoustic bands and degenerate blocks
    w_qp = self.w_q[:, iq]
    non_acoustic = [nu for nu in range(nb)
                    if self.qlanc.valid_modes_q[nu, iq]]

    # Build solve schedule: list of (band_to_solve, block_members)
    if use_mode_symmetry:
        blocks = self._find_degenerate_blocks(iq)
        solve_schedule = [(block[0], block) for block in blocks]
        n_solves = len(solve_schedule)
        n_skipped = len(non_acoustic) - n_solves
    else:
        solve_schedule = [(nu, [nu]) for nu in non_acoustic]
        n_solves = len(solve_schedule)
        n_skipped = 0

    if self.verbose:
        msg = "  Solving at q={} ({} non-acoustic bands, {} solves".format(
            iq, len(non_acoustic), n_solves)
        if n_skipped > 0:
            msg += ", {} skipped via degeneracy".format(n_skipped)
        print(msg + ")")

    # 6. Solve L_static x_i = e_i (one per degenerate block)
    G_q = np.zeros((nb, nb), dtype=np.complex128)
    total_iters = 0
    L_dense = None  # Built lazily if iterative solvers fail

    for band_i, block in solve_schedule:
        rhs = np.zeros(psi_size, dtype=np.complex128)
        rhs[band_i] = 1.0
        rhs_tilde = rhs * sqrt_mask

        # Initial guess from preconditioner
        x0 = None
        if use_preconditioner:
            x0 = apply_M_tilde(rhs_tilde)

        n_iters = [0]
        def _count(xk):
            n_iters[0] += 1

        t1 = time.time()

        # Use GMRES since L_static can be indefinite for unstable systems
        x_tilde, info = scipy.sparse.linalg.gmres(
            L_op, rhs_tilde, x0=x0, rtol=tol, maxiter=max_iters,
            M=M_op, callback=_count, callback_type='legacy')

        if info != 0:
            if self.verbose:
                print("    GMRES did not converge for band {} (info={}), "
                      "trying BiCGSTAB...".format(band_i, info))
            x_tilde, info = scipy.sparse.linalg.bicgstab(
                L_op, rhs_tilde, x0=x_tilde, rtol=tol, maxiter=max_iters,
                M=M_op)
            if info != 0:
                if not dense_fallback:
                    raise RuntimeError(
                        "Iterative solvers (GMRES, BiCGSTAB) did not "
                        "converge for band {} at iq={} (info={}). "
                        "Set dense_fallback=True to use a dense solve "
                        "(warning: O(psi_size^2) memory)."
                        .format(band_i, iq, info))
                if self.verbose:
                    print("    BiCGSTAB also did not converge "
                          "for band {} (info={}), using dense solve"
                          .format(band_i, info))
                # Build dense L matrix once, reuse for all failing bands
                if L_dense is None:
                    if self.verbose:
                        print("    Building dense L matrix "
                              "(size {})...".format(psi_size))
                    L_dense = np.zeros((psi_size, psi_size),
                                       dtype=np.complex128)
                    for j in range(psi_size):
                        e_j = np.zeros(psi_size, dtype=np.complex128)
                        e_j[j] = 1.0
                        L_dense[:, j] = apply_L_tilde(e_j)
                    L_dense = (L_dense + L_dense.conj().T) / 2
                x_tilde, _, _, _ = np.linalg.lstsq(
                    L_dense, rhs_tilde, rcond=1e-10)

        t2 = time.time()
        total_iters += n_iters[0]

        # Un-transform
        x = x_tilde * inv_sqrt_mask

        # Fill G_q using Schur's lemma for degenerate blocks.
        # G commutes with the little group, so between two d-dim
        # copies of the same irrep, G = c_cross * I_d.
        if len(block) == 1:
            # Non-degenerate: full column from R-sector
            G_q[:, band_i] = x[:nb]
        else:
            d = len(block)
            # Representative column: store the full solved column
            G_q[:, band_i] = x[:nb]
            # Non-rep columns: derive from Schur identity structure.
            # Within the block: G[block[i], block[j]] = c * delta(i,j)
            # Between two d-dim blocks A,B of same irrep:
            #   G[B[k], A[j]] = G[B[0], A[0]] * delta(k, j)
            for j in range(1, d):
                non_rep = block[j]
                # Within-block diagonal
                G_q[non_rep, non_rep] = x[band_i]
                # Cross-coupling with other same-dimension blocks
                for _, other_block in solve_schedule:
                    if other_block[0] == band_i:
                        continue
                    if len(other_block) == d:
                        # Same irrep type: shifted diagonal
                        G_q[other_block[j], non_rep] = \
                            x[other_block[0]]
                # Entries with different-dimension blocks and singlets
                # are zero by Schur (different irreps), left as 0.

        if self.verbose:
            block_str = "{}".format(block) if len(block) > 1 else ""
            print("    Band {}{}: {} iters, {:.2f}s".format(
                band_i,
                " (block {})".format(block_str) if block_str else "",
                n_iters[0], t2 - t1))

    # 7. Symmetrize G_q (should be Hermitian)
    G_q = (G_q + G_q.conj().T) / 2

    # 8. Invert to get H_q
    H_q = np.zeros((nb, nb), dtype=np.complex128)
    if len(non_acoustic) > 0:
        na = np.array(non_acoustic)
        G_sub = G_q[np.ix_(na, na)]

        cond = np.linalg.cond(G_sub)
        if cond > 1e12 and self.verbose:
            print("    WARNING: G_q ill-conditioned (cond={:.2e})".format(
                cond))

        H_sub = np.linalg.inv(G_sub)
        H_q[np.ix_(na, na)] = H_sub

    H_q = (H_q + H_q.conj().T) / 2

    if self.verbose:
        eigvals = np.sort(np.real(np.linalg.eigvalsh(H_q)))
        non_zero = eigvals[np.abs(eigvals) > 1e-10]
        if len(non_zero) > 0:
            print("  H_q eigenvalue range: [{:.6e}, {:.6e}]".format(
                non_zero[0], non_zero[-1]))
        print("  Total L applications: {}".format(total_iters))

    self.H_q_dict[iq] = H_q
    return H_q

compute_full_hessian(tol=1e-06, max_iters=500, use_preconditioner=True, dense_fallback=False, use_mode_symmetry=True)

Compute the Hessian at all q-points and return as CC.Phonons.

Parameters:

Name	Type	Description	Default
tol	float	Convergence tolerance for iterative solver.	1e-06
max_iters	int	Maximum iterations per linear solve.	500
use_preconditioner	bool	If True, use harmonic preconditioner.	True
dense_fallback	bool	If True, fall back to dense solve when iterative solvers fail. WARNING: this builds a psi_size x psi_size dense matrix, which can be very large for big supercells. Default is False.	False
use_mode_symmetry	bool	If True, exploit mode degeneracy to reduce GMRES solves.	True

Returns:

Name	Type	Description
hessian	Phonons	The free energy Hessian as a Phonons object (Ry/bohr²).

Source code in tdscha/QSpaceHessian.py

def compute_full_hessian(self, tol=1e-6, max_iters=500,
                         use_preconditioner=True,
                         dense_fallback=False,
                         use_mode_symmetry=True):
    """Compute the Hessian at all q-points and return as CC.Phonons.

    Parameters
    ----------
    tol : float
        Convergence tolerance for iterative solver.
    max_iters : int
        Maximum iterations per linear solve.
    use_preconditioner : bool
        If True, use harmonic preconditioner.
    dense_fallback : bool
        If True, fall back to dense solve when iterative solvers fail.
        WARNING: this builds a psi_size x psi_size dense matrix, which
        can be very large for big supercells. Default is False.
    use_mode_symmetry : bool
        If True, exploit mode degeneracy to reduce GMRES solves.

    Returns
    -------
    hessian : CC.Phonons.Phonons
        The free energy Hessian as a Phonons object (Ry/bohr²).
    """
    if self.verbose:
        print()
        print("=" * 50)
        print("  Q-SPACE FREE ENERGY HESSIAN")
        print("=" * 50)
        print()

    t_start = time.time()

    for iq_irr in self.irr_qpoints:
        if self.verbose:
            print("Irreducible q-point {} / {}".format(
                self.irr_qpoints.index(iq_irr) + 1,
                len(self.irr_qpoints)))
        self.compute_hessian_at_q(
            iq_irr, tol=tol, max_iters=max_iters,
            use_preconditioner=use_preconditioner,
            dense_fallback=dense_fallback,
            use_mode_symmetry=use_mode_symmetry)

    # Unfold to full BZ
    for iq_irr in self.irr_qpoints:
        H_irr = self.H_q_dict[iq_irr]
        for iq_rot, R_cart, t_cart in self.q_star_map[iq_irr]:
            if iq_rot == iq_irr:
                continue
            D = self._build_D_matrix(R_cart, t_cart, iq_irr, iq_rot)
            self.H_q_dict[iq_rot] = D @ H_irr @ D.conj().T

    t_end = time.time()
    if self.verbose:
        print()
        print("Total time: {:.1f}s".format(t_end - t_start))

    return self._build_phonons()