# Copyright © 2023-2025 HQS Quantum Simulations GmbH. All Rights Reserved.
# LICENSE PLACEHOLDER
"""Methods to calculate a spectrum taking into account local SU2 symmetries and Sz conservation."""
from __future__ import annotations
from typing import cast, TYPE_CHECKING
import numpy as np
from lattice_functions.fci import (
StateMap,
get_spin_triplets,
)
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import eigsh, ArpackNoConvergence
from scipy.linalg import eigh
import opt_einsum as oe
from itertools import product
from hqs_nmr.solver.implementations.system_tools import (
EffectiveSpinSystem,
HigherSpinState,
S2_term_higher_spin,
)
from hqs_nmr.solver.implementations.Sz_conserved_routines import (
calc_correlator_function_parallel_Sz_conserved,
set_I_minus_operator,
)
if TYPE_CHECKING:
from hqs_nmr.datatypes import NMRSolverSettings
[docs]
def determine_energy_shift_using_symmetry(
effective_spin_system: EffectiveSpinSystem, beta: float = 0.0
) -> float:
r"""Calculate the energy shift if beta > 0.
Args:
effective_spin_system: Object containing all the information on a given effective spin
system.
beta: Inverse temperature. Defaults to 0.0.
Returns:
The energy shift if beta > 0, otherwise zero.
"""
if beta > 1e-12:
return calc_groundstate_energy_using_symmetry(effective_spin_system)
else:
return 0.0
[docs]
def calc_groundstate_energy_using_symmetry(
effective_spin_system: EffectiveSpinSystem,
) -> float:
r"""Calculate the ground state energy by iteratively going through all S^2 and Sz sectors.
Args:
effective_spin_system: Object containing all the information on an effective spin system.
Returns:
The ground state energy.
"""
num_groups = effective_spin_system.number_of_groups
# j sectors of the group
all_spin_states: list[list[HigherSpinState]] = []
for n in range(num_groups):
spin_states = list(effective_spin_system.groups[n].coupled_spin_states)
all_spin_states.append(spin_states)
# There are a lot of combinations, which are permutations of one another. There might be a
# way to exploit this, however I did not think about this further yet.
all_spin_state_combinations: list[tuple[HigherSpinState, ...]] = list(
product(*all_spin_states)
)
nspins_total = effective_spin_system.reduced_size
statemap = StateMap()
current_lowest_energy = None
for higher_spin_state_combination in all_spin_state_combinations:
j_groups = []
for higher_spin_state in higher_spin_state_combination:
j_groups.append(int(2 * higher_spin_state.J))
representations = np.ones(nspins_total, dtype=int)
if len(j_groups) > 0:
representations[-len(j_groups) :] = np.array(j_groups, dtype=int)
E_S2_term = 0.0
for higher_spin_state in higher_spin_state_combination:
E_S2_term += S2_term_higher_spin(higher_spin_state)
representations = np.ones(nspins_total, dtype=int)
if len(j_groups) > 0:
representations[-len(j_groups) :] = np.array(j_groups, dtype=int)
Sz_max = np.sum(representations)
for Sz in range(-Sz_max, Sz_max + 1, 2):
# state map of the Hamiltonian in the current iteration
statemap.init_spins_S(nspins_total, representations, Sz, 0)
possible_lowest_energy = lowest_energy_Sz_and_S2_conserved_hamiltonian(
statemap,
effective_spin_system.reduced_J_coupling,
effective_spin_system.reduced_Jz,
)
possible_lowest_energy += E_S2_term
if (current_lowest_energy is None) or (possible_lowest_energy < current_lowest_energy):
current_lowest_energy = possible_lowest_energy
return cast("float", current_lowest_energy)
[docs]
def lowest_energy_Sz_and_S2_conserved_hamiltonian(
statemap: StateMap, J_coupling: np.ndarray, Jz: np.ndarray
) -> float:
r"""Find lowest energy of Hamiltonian in one Sz / S^2 sector.
Args:
statemap: State map of the Sz / S^2 sector.
J_coupling: J coupling matrix (possibly reduced due to local SU(2) symmetry).
Jz: Matrix with magnetic fields.
Returns:
Lowest eigenvalue in this sector.
"""
triplets = get_spin_triplets(
Jz=Jz,
Jp=np.zeros_like(Jz),
Jm=np.zeros_like(Jz),
hJz=J_coupling,
hJp=0.5 * J_coupling,
hJd=J_coupling,
sm=statemap,
sm_to=statemap,
)
hamiltonian_in_Sz_S2_sector = csr_matrix(
(triplets.values(), (triplets.cols(), triplets.rows())),
shape=[statemap.size(), statemap.size()],
)
if statemap.size() == 1:
lowest_energy = hamiltonian_in_Sz_S2_sector.todense()
else:
try:
lowest_energy = eigsh(
hamiltonian_in_Sz_S2_sector, which="SA", k=1, return_eigenvectors=False
)
except ArpackNoConvergence:
energies = eigh(
hamiltonian_in_Sz_S2_sector.todense(),
eigvals_only=True,
subset_by_index=[0, 0],
)
lowest_energy = energies[0]
return lowest_energy
[docs]
def diagonalize_Sz_and_S2_conserved_hamiltonian(
statemap: StateMap,
J_coupling: np.ndarray,
Jz: np.ndarray,
higher_spin_state_list: tuple[HigherSpinState, ...],
beta: float = 0.0,
) -> tuple[np.ndarray, np.ndarray]:
r"""Diagonalize hamiltonian in one Sz / S^2 sector.
Args:
statemap: State map of the Sz / S^2 sector.
J_coupling: J coupling matrix (possibly reduced due to local SU(2) symmetry).
Jz: Matrix with magnetic fields.
higher_spin_state_list: List of higher spins.
beta: Inverse temperature. Defaults to 0.0.
Returns:
Tuple[np.ndarray, np.ndarray]: eigenvalues and eigenvectors in this sector
"""
triplets = get_spin_triplets(
Jz=Jz,
Jp=np.zeros_like(Jz),
Jm=np.zeros_like(Jz),
hJz=J_coupling,
hJp=0.5 * J_coupling,
hJd=np.zeros_like(J_coupling),
sm=statemap,
sm_to=statemap,
)
OpH = csr_matrix(
(triplets.values(), (triplets.cols(), triplets.rows())),
shape=[statemap.size(), statemap.size()],
)
energies, vectors = eigh(OpH.todense())
# only energy differences in one j sector matter if we are at infinite temperature,
# hence we don't need this term, as it is the same within a j sector
if beta > 1e-12:
E_S2_term = 0.0
for higher_spin_state in higher_spin_state_list:
E_S2_term += S2_term_higher_spin(higher_spin_state)
energies += E_S2_term
return energies, vectors
[docs]
def calc_correlator_function_using_symmetry_variable_spin_number(
calc_greens_function: bool,
spin_indices: list[int],
effective_spin_system: EffectiveSpinSystem,
gyromagnetic_ratios: np.ndarray,
omegas: np.ndarray,
eta: np.ndarray,
magnetic_shift: float,
threshold_matrix_elements: float,
) -> np.ndarray:
"""Calculate an NMR spectrum for an effective spin system using Sz and local SU2 conservation.
Args:
calc_greens_function: If true the Green's function is calculated, otherwise the spectral
function.
spin_indices: List of spin indices for which to calculate the spectral function.
effective_spin_system: Object containing all the information on an effective spin system.
gyromagnetic_ratios: Array of the gyromagnetic factors per site.
omegas: Desired frequencies.
eta: Explicit spin-dependent broadening of the peaks.
magnetic_shift: Magnetic shift of the Hamiltonian.
threshold_matrix_elements: Threshold until which the Matrix elements are evaluated.
Returns:
An array of the correlator function for each frequency.
"""
num_groups = effective_spin_system.number_of_groups
nspins_total = effective_spin_system.reduced_size
number_steps = omegas.size
# j sectors of the group
all_spin_states: list[list[HigherSpinState]] = []
for n in range(num_groups):
spin_states = list(effective_spin_system.groups[n].coupled_spin_states)
all_spin_states.append(spin_states)
# There are a lot of combinations, which are permutations of one another. There might be a
# way to exploit this, however I did not think about this further yet.
all_spin_state_combinations: list[tuple[HigherSpinState, ...]] = list(
product(*all_spin_states)
)
statemap = StateMap()
statemap_minus = StateMap()
Z = 2**effective_spin_system.size
if calc_greens_function:
correlator_function_fortran = np.zeros(
(number_steps, len(spin_indices)), dtype=complex, order="F"
)
else:
correlator_function_fortran = np.zeros(
(number_steps, len(spin_indices)), dtype=float, order="F"
)
for higher_spin_state_combination in all_spin_state_combinations:
j_groups = []
for higher_spin_state in higher_spin_state_combination:
j_groups.append(int(2 * higher_spin_state.J))
representations = np.ones(nspins_total, dtype=int)
if len(j_groups) > 0:
representations[-len(j_groups) :] = np.array(j_groups, dtype=int)
Sz_max = np.sum(representations)
statemap.init_spins_S(nspins_total, representations, -Sz_max, 0)
energies, vectors = diagonalize_Sz_and_S2_conserved_hamiltonian(
statemap,
effective_spin_system.reduced_J_coupling,
effective_spin_system.reduced_Jz,
higher_spin_state_combination,
)
dim = statemap.size()
# Obtain the eigenvalues and eigenvectors of the spin Hamiltonian.
for Sz in range(-Sz_max, Sz_max, 2):
# state map of the Hamiltonain in the current iteration
statemap.init_spins_S(nspins_total, representations, Sz + 2, 0)
# state map corresponding to the sector where one ends up by applying I^-
statemap_minus.init_spins_S(nspins_total, representations, Sz, 0)
# Calculate energies in the statemap sector and keep the result from the statemap_minus
# sector.
energies_minus = energies
vectors_minus = vectors
dim_minus = dim
dim = statemap.size()
energies, vectors = diagonalize_Sz_and_S2_conserved_hamiltonian(
statemap,
effective_spin_system.reduced_J_coupling,
effective_spin_system.reduced_Jz,
higher_spin_state_combination,
)
# E_m + 0.5 * (Sz - nspins_total) * magnetic_shift # noqa: ERA001
# - E_n - 0.5 * (Sz + 2 - nspins_total) * magnetic_shift # noqa: ERA001
# = Em - En - magnetic_shift
Enm = (
np.outer(energies_minus, np.ones(dim))
- np.outer(np.ones(dim_minus), energies)
- magnetic_shift
)
# site-resolved response
list_expectation_values_left = []
for spin_index in spin_indices:
J_minus = np.zeros(nspins_total)
J_minus[spin_index] = gyromagnetic_ratios[spin_index]
operator_matrix_left = set_I_minus_operator(J_minus, statemap, statemap_minus)
expectation_values_left = oe.contract(
"ij,ik,kl -> jl",
vectors_minus,
operator_matrix_left.todense(),
vectors,
optimize=True,
)
list_expectation_values_left.append(expectation_values_left)
# We excite all spins, Im = I^m
operators_matrix_right = set_I_minus_operator(
gyromagnetic_ratios, statemap, statemap_minus
)
expectation_values_right = oe.contract(
"ij,ik,kl -> jl",
vectors_minus,
operators_matrix_right.todense(),
vectors,
optimize=True,
)
calc_correlator_function_parallel_Sz_conserved(
Enm,
list_expectation_values_left,
expectation_values_right,
threshold_matrix_elements,
eta,
omegas,
correlator_function_fortran,
)
correlator_function = np.ascontiguousarray(correlator_function_fortran.T) / Z
return correlator_function
[docs]
def calc_correlator_function_from_effective_spin_system(
spin_contribution_indices: list[int],
effective_spin_system: EffectiveSpinSystem,
gyromagnetic_ratios: np.ndarray,
omegas: np.ndarray,
eta: np.ndarray,
magnetic_shift: float,
solver_settings: NMRSolverSettings,
) -> np.ndarray:
"""Evaluate the correlator function exploiting Sz conservation and local SU2 symmetry.
Args:
spin_contribution_indices: List of indices of the spin contributions to calculate.
effective_spin_system: Effective spin system for which to calculate the spectral function.
gyromagnetic_ratios: Array of the gyromagnetic factors per site.
omegas: Desired frequencies.
eta: Explicit spin-dependent broadening of the peaks.
magnetic_shift: Magnetic shift of the hamiltonian.
solver_settings: NMRSolverSettings object containing information for the actual solver
routine.
Returns:
Array with the site-resolved spectral function for the specified spin indices.
"""
# Identify the list of effective spin contributions corresponding to the original ones.
effective_spin_contribution_indices, eta_effective_spin_system = (
effective_spin_system.identify_effective_spin_contribution_indices(
spin_contribution_indices, eta
)
)
# Evaluate correlator function in the effective system, where individual indices can code for
# a group.
effective_correlator_function = calc_correlator_function_using_symmetry_variable_spin_number(
solver_settings.calc_greens_function,
effective_spin_contribution_indices,
effective_spin_system,
gyromagnetic_ratios[effective_spin_system.indices],
omegas,
eta_effective_spin_system,
magnetic_shift,
solver_settings.threshold_matrix_elements,
)
# The correlator function needs to be transformed to the standard non-effective form.
# Create a matrix to store the corresponding result.
if solver_settings.calc_greens_function:
correlator_function = np.zeros((len(spin_contribution_indices), omegas.size), complex)
else:
correlator_function = np.zeros((len(spin_contribution_indices), omegas.size))
# Add effective spin contributions to final correlator function.
for x, spin_index in enumerate(spin_contribution_indices):
effective_spin_index = effective_spin_system.identify_effective_spin_index(spin_index)
correlator_function[x, :] = effective_correlator_function[
effective_spin_contribution_indices.index(effective_spin_index), :
] / effective_spin_system.num_spin_half(spin_index)
return correlator_function