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), False)]
                           for iq in range(self.n_q)}

save_status(fname)

Save checkpoint state (H_q_dict + metadata) to an NPZ file.

Only the master MPI rank writes the file.

Parameters:

Name	Type	Description	Default
fname	str	Path to the output file. '.npz' extension is added if missing.	required

Source code in tdscha/QSpaceHessian.py

def save_status(self, fname):
    """Save checkpoint state (H_q_dict + metadata) to an NPZ file.

    Only the master MPI rank writes the file.

    Parameters
    ----------
    fname : str
        Path to the output file. '.npz' extension is added if missing.
    """
    Parallel.barrier()
    if Parallel.am_i_the_master():
        if not fname.lower().endswith(".npz"):
            fname += ".npz"

        uci = self.qlanc.uci_structure
        save_dict = dict(
            T=np.float64(self.qlanc.T),
            n_q=np.int64(self.n_q),
            n_bands=np.int64(self.n_bands),
            ignore_v3=np.bool_(self.ignore_v3),
            ignore_v4=np.bool_(self.ignore_v4),
            unit_cell=uci.unit_cell,
            coords=uci.coords,
            atom_types=uci.get_atomic_types(),
        )

        completed = sorted(self.H_q_dict.keys())
        save_dict["completed_irr_qpoints"] = np.array(completed,
                                                       dtype=np.int64)
        if len(completed) > 0:
            save_dict["H_q_keys"] = np.array(completed, dtype=np.int64)
            save_dict["H_q_values"] = np.stack(
                [self.H_q_dict[iq] for iq in completed])
        else:
            save_dict["H_q_keys"] = np.array([], dtype=np.int64)
            save_dict["H_q_values"] = np.zeros((0, self.n_bands,
                                                 self.n_bands),
                                                dtype=np.complex128)

        np.savez_compressed(fname, **save_dict)

load_status(fname)

Load checkpoint state from an NPZ file and validate metadata.

Master rank reads the file and broadcasts to all MPI ranks. Populates H_q_dict with the saved Hessian matrices.

Parameters:

Name	Type	Description	Default
fname	str	Path to the NPZ file. '.npz' extension is added if missing.	required

Returns:

Name	Type	Description
completed_iqs	set of int	Set of irreducible q-point indices that have been completed.

Raises:

Type	Description
ValueError	If any metadata field in the checkpoint does not match the current object.

Source code in tdscha/QSpaceHessian.py

def load_status(self, fname):
    """Load checkpoint state from an NPZ file and validate metadata.

    Master rank reads the file and broadcasts to all MPI ranks.
    Populates ``H_q_dict`` with the saved Hessian matrices.

    Parameters
    ----------
    fname : str
        Path to the NPZ file. '.npz' extension is added if missing.

    Returns
    -------
    completed_iqs : set of int
        Set of irreducible q-point indices that have been completed.

    Raises
    ------
    ValueError
        If any metadata field in the checkpoint does not match
        the current object.
    """
    Parallel.barrier()
    if not fname.lower().endswith(".npz"):
        fname += ".npz"

    if Parallel.am_i_the_master():
        if not os.path.exists(fname):
            raise IOError("Checkpoint file not found: {}".format(fname))
        data = dict(np.load(fname, allow_pickle=True))
    else:
        data = None
    data = Parallel.broadcast(data)

    # Validate metadata
    uci = self.qlanc.uci_structure
    mismatches = []

    if not np.isclose(float(data["T"]), self.qlanc.T, atol=1e-10):
        mismatches.append("T: checkpoint={}, current={}".format(
            float(data["T"]), self.qlanc.T))
    if int(data["n_q"]) != self.n_q:
        mismatches.append("n_q: checkpoint={}, current={}".format(
            int(data["n_q"]), self.n_q))
    if int(data["n_bands"]) != self.n_bands:
        mismatches.append("n_bands: checkpoint={}, current={}".format(
            int(data["n_bands"]), self.n_bands))
    if bool(data["ignore_v3"]) != self.ignore_v3:
        mismatches.append("ignore_v3: checkpoint={}, current={}".format(
            bool(data["ignore_v3"]), self.ignore_v3))
    if bool(data["ignore_v4"]) != self.ignore_v4:
        mismatches.append("ignore_v4: checkpoint={}, current={}".format(
            bool(data["ignore_v4"]), self.ignore_v4))
    if not np.allclose(data["unit_cell"], uci.unit_cell, atol=1e-10):
        mismatches.append("unit_cell differs")
    if not np.allclose(data["coords"], uci.coords, atol=1e-10):
        mismatches.append("atomic coordinates differ")
    current_types = uci.get_atomic_types()
    if not np.array_equal(data["atom_types"], current_types):
        mismatches.append("atom_types differ")

    if mismatches:
        raise ValueError(
            "Checkpoint metadata mismatch:\n  " +
            "\n  ".join(mismatches))

    # Restore H_q_dict
    keys = data["H_q_keys"].astype(int)
    values = data["H_q_values"]
    for i, iq in enumerate(keys):
        self.H_q_dict[int(iq)] = values[i]

    completed = set(data["completed_irr_qpoints"].astype(int).tolist())
    return completed

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]
            # Zero out entries for modes in degenerate blocks
            # (different irreps → zero by Schur's lemma).
            # This prevents GMRES noise from breaking degeneracy
            # after Hermitian symmetrization.
            for _, other_block in solve_schedule:
                if len(other_block) >= 2:
                    for m in other_block:
                        G_q[m, band_i] = 0.0
        else:
            d = len(block)
            # Extract Schur-consistent scalars only (not the full
            # noisy GMRES column). This ensures all columns within
            # the degenerate block are filled identically, preserving
            # perfect block structure and preventing degeneracy
            # breaking after symmetrization.
            c_diag = x[band_i]  # within-block diagonal constant

            # Fill ALL columns in this block (including rep) uniformly
            for j in range(d):
                col = block[j]
                # Within-block diagonal
                G_q[col, col] = c_diag
                # 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], col] = 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, checkpoint=None)

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
checkpoint	str or None	If a string path, auto-save after each completed irreducible q-point and resume from the checkpoint file on restart. The '.npz' extension is added automatically if missing.	None

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,
                         checkpoint=None):
    """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.
    checkpoint : str or None
        If a string path, auto-save after each completed irreducible
        q-point and resume from the checkpoint file on restart.
        The '.npz' extension is added automatically if missing.

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

    # Load checkpoint if it exists
    completed_iqs = set()
    if checkpoint is not None:
        ckpt_file = checkpoint
        if not ckpt_file.lower().endswith(".npz"):
            ckpt_file += ".npz"
        if os.path.exists(ckpt_file):
            completed_iqs = self.load_status(checkpoint)
            if self.verbose:
                print("Restarting: {} q-points already done".format(
                    len(completed_iqs)))

    for iq_irr in self.irr_qpoints:
        if iq_irr in completed_iqs:
            if self.verbose:
                print("Irreducible q-point {} / {} (loaded from "
                      "checkpoint)".format(
                          self.irr_qpoints.index(iq_irr) + 1,
                          len(self.irr_qpoints)))
            continue
        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)
        if checkpoint is not None:
            self.save_status(checkpoint)

    # 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, is_tr 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,
                                     time_reversal=is_tr)
            if is_tr:
                # Time-reversal: H(-Rq) = D_TR @ conj(H(q)) @ D_TR^dag
                self.H_q_dict[iq_rot] = D @ np.conj(H_irr) @ D.conj().T
            else:
                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()

from_qspace_lanczos(qlanc, verbose=True, use_symmetries=True) classmethod

Create a QSpaceHessian from an existing QSpaceLanczos object.

Parameters:

Name	Type	Description	Default
qlanc	QSpaceLanczos	An initialized QSpaceLanczos object (e.g., from load_distributed_tdscha).	required
verbose	bool	If True, print progress information.	True
use_symmetries	bool	If True, use symmetries to find irreducible q-points.	True

Returns:

Type	Description
QSpaceHessian	A QSpaceHessian object with the given qlanc as its Lanczos engine.

Source code in tdscha/QSpaceHessian.py

@classmethod
def from_qspace_lanczos(cls, qlanc, verbose=True, use_symmetries=True):
    """Create a QSpaceHessian from an existing QSpaceLanczos object.

    Parameters
    ----------
    qlanc : QSpaceLanczos
        An initialized QSpaceLanczos object (e.g., from load_distributed_tdscha).
    verbose : bool
        If True, print progress information.
    use_symmetries : bool
        If True, use symmetries to find irreducible q-points.

    Returns
    -------
    QSpaceHessian
        A QSpaceHessian object with the given qlanc as its Lanczos engine.
    """
    # Create instance using __new__ (bypass __init__)
    hess = cls.__new__(cls)

    # Set the qlanc reference (shares data, no copy)
    hess.qlanc = qlanc

    # Copy shortcuts from qlanc
    hess.ensemble = qlanc.ensemble
    hess.n_q = qlanc.n_q
    hess.n_bands = qlanc.n_bands
    hess.q_points = qlanc.q_points
    hess.w_q = qlanc.w_q
    hess.pols_q = qlanc.pols_q

    # Initialize caches to None
    hess.H_q_dict = {}
    hess.irr_qpoints = None
    hess.q_star_map = None
    hess._cached_w_qp_sq = None
    hess._cached_inv_w_qp_sq = None
    hess._cached_inv_lambda = None
    hess._cached_lambda = None
    hess._cached_yw_outer = None
    hess._has_symmetry_data = False

    # Set verbose flag
    hess.verbose = verbose

    # Copy ignore flags from qlanc
    hess.ignore_v3 = qlanc.ignore_v3
    hess.ignore_v4 = qlanc.ignore_v4

    # Initialize symmetries if requested
    if use_symmetries and __SPGLIB__:
        hess._find_irreducible_qpoints()
    else:
        hess.irr_qpoints = list(range(hess.n_q))
        hess.q_star_map = {iq: [(iq, np.eye(3), np.zeros(3), False)]
                           for iq in range(hess.n_q)}

    return hess

load_distributed_hessian(data_dir, population_id, dyn, T, lo_to_split=None, use_symmetries=True, n_configs=None, final_dyn=None, final_T=None, verbose=True, ignore_v3=False, ignore_v4=False, **kwargs)

Load QSpaceHessian with distributed configurations across MPI ranks.

Loads the ensemble on master rank only, then distributes configuration data across all ranks. Combines load_distributed_tdscha and QSpaceHessian.from_qspace_lanczos.

Parameters:

Name	Type	Description	Default
data_dir	str	Directory containing the ensemble data files.	required
population_id	int	Population ID of the ensemble to load.	required
dyn	Phonons	Dynamical matrix object (used to create the ensemble).	required
T	float	Temperature in Kelvin.	required
lo_to_split	str, ndarray, or None	LO-TO splitting mode.	None
use_symmetries	bool	If True, use q-space symmetries.	True
n_configs	int or None	Number of configs to load.	None
final_dyn	Phonons	Final dynamical matrix from the SSCHA calculation. If provided, the ensemble weights are updated using update_weights(final_dyn, final_T).	None
final_T	float	Temperature for weight updates. Defaults to T if not specified.	None
verbose	bool	If True, print progress information.	True
ignore_v3	bool	If True, exclude cubic (D3) anharmonic contributions.	False
ignore_v4	bool	If True, exclude quartic (D4) anharmonic contributions.	False
**kwargs		Additional arguments passed to QSpaceLanczos.	{}

Returns:

Type	Description
QSpaceHessian	QSpaceHessian with distributed ensemble. Already initialized.

Usage

mpirun -np 8 python your_script.py

Example with weight update (recommended for production): final_dyn = CC.Phonons.Phonons("final_dyn_", nqirr=3) hess = load_distributed_hessian( "data/", 1, initial_dyn, 300, final_dyn=final_dyn, final_T=300 ) hessian_dyn = hess.compute_full_hessian()

Example without weight update (for testing/debugging): hess = load_distributed_hessian("data/", 1, dyn, 300) hessian_dyn = hess.compute_full_hessian()

Source code in tdscha/QSpaceHessian.py

def load_distributed_hessian(data_dir, population_id, dyn, T, lo_to_split=None,
                            use_symmetries=True, n_configs=None,
                            final_dyn=None, final_T=None,
                            verbose=True, ignore_v3=False, ignore_v4=False,
                            **kwargs):
    """Load QSpaceHessian with distributed configurations across MPI ranks.

    Loads the ensemble on master rank only, then distributes configuration data
    across all ranks. Combines load_distributed_tdscha and QSpaceHessian.from_qspace_lanczos.

    Parameters
    ----------
    data_dir : str
        Directory containing the ensemble data files.
    population_id : int
        Population ID of the ensemble to load.
    dyn : CC.Phonons.Phonons
        Dynamical matrix object (used to create the ensemble).
    T : float
        Temperature in Kelvin.
    lo_to_split : str, ndarray, or None
        LO-TO splitting mode.
    use_symmetries : bool
        If True, use q-space symmetries.
    n_configs : int or None
        Number of configs to load.
    final_dyn : CC.Phonons.Phonons, optional
        Final dynamical matrix from the SSCHA calculation. If provided,
        the ensemble weights are updated using update_weights(final_dyn, final_T).
    final_T : float, optional
        Temperature for weight updates. Defaults to T if not specified.
    verbose : bool
        If True, print progress information.
    ignore_v3 : bool
        If True, exclude cubic (D3) anharmonic contributions.
    ignore_v4 : bool
        If True, exclude quartic (D4) anharmonic contributions.
    **kwargs
        Additional arguments passed to QSpaceLanczos.

    Returns
    -------
    QSpaceHessian
        QSpaceHessian with distributed ensemble. Already initialized.

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

    Example with weight update (recommended for production):
        final_dyn = CC.Phonons.Phonons("final_dyn_", nqirr=3)
        hess = load_distributed_hessian(
            "data/", 1, initial_dyn, 300,
            final_dyn=final_dyn, final_T=300
        )
        hessian_dyn = hess.compute_full_hessian()

    Example without weight update (for testing/debugging):
        hess = load_distributed_hessian("data/", 1, dyn, 300)
        hessian_dyn = hess.compute_full_hessian()
    """
    from tdscha.QSpaceLanczos import load_distributed_tdscha

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

    hess = QSpaceHessian.from_qspace_lanczos(
        qlanc, verbose=verbose, use_symmetries=use_symmetries
    )
    hess.ignore_v3 = ignore_v3
    hess.ignore_v4 = ignore_v4
    hess.init(use_symmetries=use_symmetries)

    return hess