# Copyright © 2022-2025 HQS Quantum Simulations GmbH. All Rights Reserved.
"""Routines to calculate a spectrum, if only the Sz conservation is taken into account."""
from __future__ import annotations
from typing import Union, TYPE_CHECKING
import numpy as np
from scipy.linalg import eigh
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import ArpackNoConvergence, eigsh
from lattice_functions.fci import (
StateMap,
csr_matrix_eigen,
get_spin_triplets,
get_spin_triplets_from_expression,
)
from lattice_functions.nmr_tools import (
calculate_spectral_function_for_nmr_spin_dependent,
calculate_greens_function_for_nmr_spin_dependent,
)
if TYPE_CHECKING:
from lattice_functions import Fermions as fermions
[docs]
def determine_energy_shift(
Jz: np.ndarray, hJz: np.ndarray, hJp: np.ndarray, beta: float = 0.0
) -> float:
"""Determine the energy shift corresponding to the overall groundstate energy.
Note that the ground state energy is only calculated if beta is larger than zero.
Args:
Jz: Array with S^z terms. dim: (number_spins).
hJz: Array with the S^z S^z coupling. dim: (number_spins x number_spins).
hJp: Array with the S^+ S^- coupling. dim: (number_spins x number_spins).
beta: Inverse temperature. Defaults to 0.0.
Returns:
The ground state energy if beta > 0, otherwise zero.
"""
if beta > 1e-12:
return calc_groundstate_energy(Jz, hJz, hJp)
else:
return 0.0
[docs]
def calc_groundstate_energy(Jz: np.ndarray, hJz: np.ndarray, hJp: np.ndarray) -> float:
r"""Calculate ground state by iterating over all Sz sectors.
Args:
Jz: Array with S^z terms. dim: (number_spins).
hJz: Array with the S^z S^z coupling. dim: (number_spins x number_spins).
hJp: Array with the S^+ S^- coupling. dim: (number_spins x number_spins).
Returns:
The ground state energy.
"""
nspins_total = Jz.shape[0]
statemap = StateMap()
statemap.init_spins(nspins_total, -nspins_total, 0)
current_lowest_energy = lowest_energy_in_Sz_sector(statemap, Jz, hJz, hJp)
# Obtain the eigenvalues and eigenvectors of the spin Hamiltonian.
for Sz in range(-nspins_total, nspins_total, 2):
# state map of the Hamiltonian in the current iteration
statemap.init_spins(nspins_total, Sz + 2, 0)
possible_lowest_energy = lowest_energy_in_Sz_sector(statemap, Jz, hJz, hJp)
if possible_lowest_energy < current_lowest_energy:
current_lowest_energy = possible_lowest_energy
return current_lowest_energy
[docs]
def lowest_energy_in_Sz_sector(
statemap: StateMap, Jz: np.ndarray, hJz: np.ndarray, hJp: np.ndarray
) -> float:
r"""Builds the Hamiltonian in one Sz sector and finds the lowest eigenvalue.
Args:
statemap: State map of the Sz sector.
Jz: Array with S^z terms. dim: (number_spins).
hJz: Array with the S^z S^z coupling. dim: (number_spins x number_spins).
hJp: Array with the S^+ S^- coupling. dim: (number_spins x number_spins).
Returns:
Returns the smallest algebraic eigenvalue.
"""
triplets = get_spin_triplets(
Jz=Jz,
Jp=np.zeros_like(Jz),
Jm=np.zeros_like(Jz),
hJz=hJz,
hJp=0.5 * hJp,
hJd=np.zeros_like(hJz),
sm=statemap,
sm_to=statemap,
)
hamiltonian_in_Sz_sector = csr_matrix(
(triplets.values(), (triplets.cols(), triplets.rows())),
shape=[statemap.size(), statemap.size()],
)
if statemap.size() == 1:
lowest_energy = hamiltonian_in_Sz_sector.todense()
else:
try:
lowest_energy = eigsh(
hamiltonian_in_Sz_sector, which="SA", k=1, return_eigenvectors=False
)
except ArpackNoConvergence:
energies = eigh(
hamiltonian_in_Sz_sector.todense(),
eigvals_only=True,
subset_by_index=[0, 0],
)
lowest_energy = energies[0]
return lowest_energy
[docs]
def set_I_minus_operator(
J_minus: np.ndarray, statemap: StateMap, statemap_minus: StateMap
) -> csr_matrix_eigen:
r"""Sets the operator I^-.
Args:
J_minus: Array with I^- terms. dim: (number_spins).
statemap: State map of the current Sz sector.
statemap_minus: State map of the Sz sector after applying I^-.
Returns:
Returns the I^- operator applied to the sites specified in J_minus.
"""
nspins_total = J_minus.size
op_I_minus = csr_matrix_eigen(statemap_minus.size(), statemap.size())
op_I_minus.set(
Jz=np.zeros(nspins_total),
Jp=np.zeros(nspins_total),
Jm=J_minus,
hJz=np.zeros((nspins_total, nspins_total)),
hJp=np.zeros((nspins_total, nspins_total)),
hJd=np.zeros((nspins_total, nspins_total)),
sm=statemap,
sm_to=statemap_minus,
)
return op_I_minus
[docs]
def diagonalize_Sz_conserved_hamiltonian(
statemap: StateMap, Jz: np.ndarray, hJz: np.ndarray, hJp: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
r"""Builds and diagonalizes the Hamiltonian in one Sz sector.
Args:
statemap: State map of the Sz sector.
Jz: Array with S^z terms. dim: (number_spins).
hJz: Array with the S^z S^z coupling. dim: (number_spins x number_spins).
hJp: Array with the S^+ S^- coupling. dim: (number_spins x number_spins).
Returns:
Returns the Eigenvalues, and Eigenvectors in one Sz sector.
"""
nspins_total = Jz.shape[0]
Sz_conserved_hamiltonian = csr_matrix_eigen(statemap.size(), statemap.size())
Sz_conserved_hamiltonian.set(
Jz=Jz,
Jp=np.zeros(nspins_total),
Jm=np.zeros(nspins_total),
hJz=hJz,
hJp=0.5 * hJp,
hJd=np.zeros((nspins_total, nspins_total)),
sm=statemap,
sm_to=statemap,
)
energies, vectors = np.linalg.eigh(Sz_conserved_hamiltonian.todense())
return energies, vectors
[docs]
def calc_correlator_function_parallel_Sz_conserved(
Enm: np.ndarray,
list_expectation_values_left: list[np.ndarray],
expectation_values_right: np.ndarray,
threshold_matrix_elements: float,
eta: np.ndarray,
omegas: np.ndarray,
correlator_function_fortran: np.ndarray,
) -> None:
"""Calculate spectral function for one individual spin.
Args:
Enm: Energy differences between the different eigenenergies.
list_expectation_values_left: Expectation values for the ladder operators on the left of
the correlator.
expectation_values_right: Expectation values of total spin ladder operator
threshold_matrix_elements: Threshold until which the Matrix elements are evaluated.
eta: Explicit spin-dependent broadening of the peaks.
omegas: Frequency range.
correlator_function_fortran: Empty array to be overwritten with the spin resolved spectral
function. Has to fortran style.
"""
for x in range(len(list_expectation_values_left)):
# precalculate the terms <n|I_x^+|m><m|I^+|n> to be used in the loop over the frequencies
expectation_values = list_expectation_values_left[x] * expectation_values_right * eta[x]
idx_n, idx_m = np.where(np.abs(expectation_values) > threshold_matrix_elements)
Enm_x = Enm[idx_n, idx_m]
expectation_values = expectation_values[idx_n, idx_m]
if correlator_function_fortran.dtype == complex:
calculate_greens_function_for_nmr_spin_dependent(
x, eta[x], omegas, Enm_x, expectation_values, correlator_function_fortran
)
else:
calculate_spectral_function_for_nmr_spin_dependent(
x, eta[x], omegas, Enm_x, expectation_values, correlator_function_fortran
)
[docs]
def evaluate_braket(
operator: fermions.ExpressionSpinful,
bra: np.ndarray,
ket: np.ndarray,
statemap_bra: StateMap,
statemap_ket: StateMap,
) -> np.ndarray:
"""Evaluate expressions of the form <n|O|m>, where O is some operator.
Args:
operator: Fermion expression for which to evaluate the braket.
bra: Bra vectors, stored in matrix.
ket: Ket vectors, stored in matrix.
statemap_bra: State map for the bra vector.
statemap_ket: State map for the ket vector.
"""
if operator.size() == 0:
return np.zeros((statemap_bra.size(), statemap_ket.size()))
operator_triplets = get_spin_triplets_from_expression(operator, statemap_ket, statemap_bra)
operators_matrix = csr_matrix_eigen(statemap_bra.size(), statemap_ket.size())
operators_matrix.set(operator_triplets)
tmp = np.zeros((operators_matrix.shape[0], ket.shape[1]))
operators_matrix.dot(ket, tmp)
return bra.T @ tmp
[docs]
def calc_correlator_function_variable_spin_number(
calc_greens_function: bool,
spin_indices: list[int],
site_resolved_operators: Union[fermions.ExpressionSpinful, list[fermions.ExpressionSpinful]],
operator_sum: fermions.ExpressionSpinful,
Jz: np.ndarray,
hJz: np.ndarray,
hJp: np.ndarray,
omegas: np.ndarray,
eta: np.ndarray,
magnetic_shift: float,
threshold_matrix_elements: float,
verbose: int,
) -> np.ndarray:
r"""Calculate specified spin contributions of an NMR spectrum.
The spin contributions are calculated exploiting Sz conservation through a direct resolvent
approach.
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.
site_resolved_operators: List of site-resolved operators or operators on one site.
operator_sum: Sum of operators acting on each site.
Jz: Array with Sz terms.
hJz: Array with the Sz Sz coupling.
hJp: Array with the S+ S- coupling.
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.
verbose: Level of verbosity of output.
Returns:
An array of the spectral function for each frequency for this cluster.
"""
if isinstance(site_resolved_operators, list):
site_resolved_operator_list = site_resolved_operators
else:
site_resolved_operator_list = [site_resolved_operators]
nspins_total = Jz.shape[0]
Z = 2**nspins_total
number_steps = omegas.size
statemap = StateMap()
statemap_minus = StateMap()
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"
)
statemap.init_spins(nspins_total, -nspins_total, 0)
if verbose > 1:
print(f"{nspins_total} Start Sz: {-nspins_total}: {-nspins_total / nspins_total}")
energies, vectors = diagonalize_Sz_conserved_hamiltonian(statemap, Jz, hJz, hJp)
dim = statemap.size()
# Obtain the eigenvalues and eigenvectors of the spin Hamiltonian.
for Sz in range(-nspins_total, nspins_total, 2):
if verbose > 1:
print(f"{nspins_total} Start Sz: {Sz + 2}: {Sz + 2 / nspins_total}")
# state map of the Hamiltonian in the current iteration
statemap.init_spins(nspins_total, Sz + 2, 0)
# state map corresponding to the sector where one ends up by applying I^-
statemap_minus.init_spins(nspins_total, Sz, 0)
if verbose > 1:
print(f"{nspins_total} init Sz: {Sz}: {Sz / nspins_total}")
# calculate energies in the statemap sector,
# keep the result from the statemap_minus sector
energies_minus = energies
vectors_minus = vectors
dim_minus = dim
dim = statemap.size()
energies, vectors = diagonalize_Sz_conserved_hamiltonian(statemap, Jz, hJz, hJp)
# 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
)
list_expectation_values_left = [
evaluate_braket(
site_resolved_operator_list[spin_index],
vectors_minus,
vectors,
statemap_minus,
statemap,
)
for spin_index in spin_indices
]
expectation_values_right = evaluate_braket(
operator_sum, vectors_minus, vectors, statemap_minus, statemap
)
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