# Copyright © 2024-2025 HQS Quantum Simulations GmbH. All Rights Reserved.
# LICENSE PLACEHOLDER
"""Functions that can identify for a particular spin the most strongly coupled spins.
Used to construct spin-dependent clusters.
"""
from __future__ import annotations
from typing import TYPE_CHECKING
import numpy as np
if TYPE_CHECKING:
from hqs_nmr.datatypes import NMRSolverSettings
from lattice_functions.nmr_tools import (
fiedler_partitioning,
)
[docs]
def get_weight_matrix(hJp: np.ndarray, Jz: np.ndarray, weight_offset: float) -> np.ndarray:
r"""Create weight matrix according to perturbation theory.
.. math::
w_{x,y} = \frac{|J_{x,y}|^2}{|J^z_x - J^z_y| + d}
Args:
hJp: J-coupling matrix.
Jz: Jz vector.
weight_offset: Small factor to avoid division by zero.
Returns:
Matrix with weights from perturbation theory.
"""
weights = np.zeros_like(hJp)
for x in range(0, Jz.size):
for y in range(0, Jz.size):
denominator = np.abs(Jz[x] - Jz[y]) + weight_offset
if denominator < 1e-10:
raise ValueError(
"Choose a weight offset unequal to zero to avoid division by zero."
)
weights[x, y] = hJp[x, y] ** 2 / denominator
return weights
[docs]
def spin_dependent_cluster_finder(
spin_index: int,
weights: np.ndarray,
J_coupling: np.ndarray,
max_cluster_size: int,
minimize_cluster: bool = False,
) -> list[int]:
r"""For one spin return sorted list with indices of spins its strongest coupled with.
Args:
spin_index: Index of the spin for which to find the coupled spins.
weights: Weight matrix of the effective spin system.
J_coupling: Matrix with J-coupling values.
max_cluster_size: Maximum number of spins the spin is allowed to couple to.
minimize_cluster: If True, it is checked, if independent partitions can be
identified in the cluster. Defaults to False.
Returns:
Integer array with sorted indices that make up the cluster for the spin.
"""
idx = np.argsort(weights[spin_index, :])
idx = np.delete(idx, np.where(idx == spin_index))
indices = np.zeros(max_cluster_size, dtype=int)
indices[0] = spin_index
indices[1:] = np.flip(idx[-(max_cluster_size - 1) :])
if minimize_cluster:
# check if there are independent partitions in the molecule
fiedler_partitions = fiedler_partitioning(J_coupling[indices, :][:, indices], verbose=0)
# only keep partition with spin_index
for partition in fiedler_partitions:
if np.where(indices == spin_index) in partition:
indices = indices[partition]
break
sorted_indices = [int(index) for index in sorted(indices)]
return sorted_indices
[docs]
def spin_dependent_cluster_list(
weights: np.ndarray,
J_coupling: np.ndarray,
max_cluster_size: int,
number_spins: int,
verbose: int = 0,
minimize_cluster: bool = False,
) -> tuple[list, list]:
r"""Identify clusters for each spin.
Args:
weights: Weight matrix of the effective spin system.
J_coupling: Matrix with J-coupling values.
max_cluster_size: Maximum number of spins the spin is allowed to couple to.
number_spins: Number of spin half in the full system.
verbose: Verbosity level of output.
minimize_cluster: If True, it is checked, if partitions can be identified in
the cluster. Defaults to False.
Returns:
List of lists of spins that are associated with identical clusters.
List of integer arrays with sorted indices that make up the clusters.
"""
effective_max_cluster_size = min(number_spins, max_cluster_size)
cluster_list: list[list[int]] = []
spins_with_the_same_cluster: list[list[int]] = []
for spin_index in range(number_spins):
new_cluster = spin_dependent_cluster_finder(
spin_index,
weights,
J_coupling,
effective_max_cluster_size,
minimize_cluster,
)
if new_cluster in cluster_list:
index_new_cluster = cluster_list.index(new_cluster)
spins_with_the_same_cluster[index_new_cluster].append(spin_index)
else:
cluster_list.append(new_cluster)
spins_with_the_same_cluster.append([spin_index])
if verbose > 0:
for cluster_index, cluster in enumerate(cluster_list):
print(
"cluster for spin(s)",
spins_with_the_same_cluster[cluster_index],
": ",
cluster,
)
return spins_with_the_same_cluster, cluster_list
[docs]
def identify_spin_dependent_clusters(
Jz: np.ndarray, hJp: np.ndarray, solver_settings: NMRSolverSettings
) -> tuple[list, list]:
"""Identify for each spin the best cluster and group spins with the same cluster together.
Args:
Jz: Array with S^z terms. dim: (number_spins).
hJp: Array with the S^+ S^- coupling. dim: (number_spins x number_spins).
solver_settings: NMRSolverSettings object containing information on the cluster and solver
method.
Returns:
List of lists of spins that are associated with identical clusters.
List of integer arrays with sorted indices that make up the clusters.
"""
# cluster Hamiltonian
number_spins = Jz.size
full_cluster = list(range(number_spins))
if (number_spins <= solver_settings.max_cluster_size) or solver_settings.solve_exactly:
return [full_cluster], [full_cluster]
weights = get_weight_matrix(hJp, Jz, solver_settings.weight_offset)
return spin_dependent_cluster_list(
weights,
hJp,
solver_settings.max_cluster_size,
number_spins,
verbose=solver_settings.verbose,
minimize_cluster=True,
)