Source code for hqs_nmr.solver.implementations.local_SU2_conserved_routines

# 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