Source code for hqs_nmr.hamiltonian_tools

# Copyright © 2024-2025 HQS Quantum Simulations GmbH. All Rights Reserved.

# LICENSE PLACEHOLDER

"""Tools to obtain and manipulate a NMR Hamiltonian."""

from __future__ import annotations
from typing import cast, Any, Dict

from scipy.sparse import coo_matrix
from numpy import typing as npt
import numpy as np

from struqture_py.spins import SpinHamiltonianSystem, PauliProduct

from hqs_nmr_parameters import GYROMAGNETIC_RATIOS, nmr_hamiltonian, Isotope

from hqs_nmr import datatypes, utils
from hqs_nmr.solver.implementations.system_tools import identify_symmetry_groups
from hqs_nmr.spin_dependent_clustering import (
    identify_spin_dependent_clusters,
)

DEFAULT_REFERENCE_ISOTOPE = Isotope(1, "H")
GyromagneticRatios = cast(dict[Isotope, float], GYROMAGNETIC_RATIOS)


[docs] def extract_coupling( H: SpinHamiltonianSystem, ) -> tuple[np.ndarray, np.ndarray, np.ndarray]: """Extract the coupling matrices from a Struqture spin Hamiltonian. Args: H: Struqture Spin Hamiltonian. Returns: Array with S^z term with prefactor 2, array with the S^z-S^z coupling, array with the S^+-S^- coupling. """ nspins_total = H.number_spins() Jz = np.zeros(nspins_total) hJz = np.zeros((nspins_total, nspins_total)) hJp = np.zeros((nspins_total, nspins_total)) for r in range(nspins_total): Jz[r] = 2 * H.get(PauliProduct().z(r)) # 2 Jz S_z = Jz sigma_z for s in range(nspins_total): if r != s: hJz[r, s] = 2 * H.get(PauliProduct().z(r).z(s)) hJp[r, s] = 2 * H.get(PauliProduct().x(r).x(s)) return Jz, hJz, hJp
[docs] def extract_cluster_nmr_parameters( molecule_parameters: datatypes.NMRParameters, cluster_indices: list[int] ) -> datatypes.NMRParameters: """Extract the NMR parameters of a cluster from the NMR parameters of the molecule. Args: molecule_parameters: The molecular isotopes, shifts and J-coupling values. cluster_indices: Indices of the cluster. Returns: An object storing the NMR data of the cluster Hamiltonian. """ cluster_isotopes = [] cluster_shifts = [] cluster_j_couplings = [] for i in cluster_indices: cluster_isotopes.append(molecule_parameters.isotopes[i]) cluster_shifts.append(molecule_parameters.shifts[i]) for j_coupling in molecule_parameters.j_couplings: if (j_coupling[0][0] in cluster_indices) and (j_coupling[0][1] in cluster_indices): cluster_j_coupling = ( ( cluster_indices.index(j_coupling[0][0]), cluster_indices.index(j_coupling[0][1]), ), j_coupling[1], ) cluster_j_couplings.append(cluster_j_coupling) return datatypes.NMRParameters( isotopes=cluster_isotopes, shifts=cluster_shifts, j_couplings=cluster_j_couplings, )
[docs] def sh_to_matrix(H: SpinHamiltonianSystem) -> npt.NDArray[np.complexfloating]: """Convert a Struqture spin Hamiltonian to its matrix representation. Args: H: Struqture Spin Hamiltonian Returns: The Spin Hamiltonian as a dense matrix (numpy array). """ H_coo = H.sparse_matrix_coo() return coo_matrix(H_coo).toarray()
[docs] def cluster_molecule( molecule_parameters: datatypes.NMRParameters, gyromagnetic_ratios: dict[Isotope, float], field: float, solver_settings: datatypes.NMRSolverSettings, reference_isotope: Isotope, ) -> list[Dict[str, Any]]: """Cluster the molecule into spin-dependent cluster. Args: molecule_parameters: The molecular isotopes, shifts and J-coupling values. gyromagnetic_ratios: Dictionary of gyromagnetic ratios in rad / (T s). field: Magnetic field in T. solver_settings: NMRSolverSettings object containing information on the clustering and solver methods. reference_isotope: Isotope specified as Isotope(mass, symbol) to define the frequency (w=gamma*field) of the rotating frame. Defaults to Isotope(1, 'H'). Returns: List with a dictionary defining each spin-dependent cluster. """ # [hamiltonian]: rad / s hamiltonian = nmr_hamiltonian( molecule_parameters, field, reference_isotope=reference_isotope, gyromagnetic_ratios=gyromagnetic_ratios, ) # Extract Hamiltonian matrices from spin Hamiltonian. Jz, _, hJp = extract_coupling(hamiltonian) # Identify the spin-dependent clusters. spins_with_same_cluster, cluster_list = identify_spin_dependent_clusters( Jz, hJp, solver_settings ) # Create cluster dictionaries. cluster_dict_list = [] for spin_contribution_indices, cluster_indices in zip(spins_with_same_cluster, cluster_list): # Extract cluster nmr parameters from the molecule nmr parameters. cluster_parameters = extract_cluster_nmr_parameters(molecule_parameters, cluster_indices) # Reduce spin contributions indices to only include relevant contributions. relevant_spin_contribution_indices = utils.reduce_contributions_to_relevant_spins( spin_contribution_indices, molecule_parameters, reference_isotope, solver_settings ) # Evaluate which indices the spin contributions have with respect to the cluster indices. spin_contribution_indices_in_cluster = [ cluster_indices.index(i) for i in relevant_spin_contribution_indices ] # Check if local SU2 symmetry should be exploited for this cluster cluster_size = len(cluster_indices) use_local_su2 = utils.check_local_su2(cluster_size, solver_settings) # Save the information on the cluster in a dictionary. cluster_dict = { "spin_contribution_indices_in_molecule": relevant_spin_contribution_indices, "spin_contribution_indices_in_cluster": spin_contribution_indices_in_cluster, "cluster_indices": cluster_indices, "nmr_parameters": cluster_parameters, "use_local_su2": use_local_su2, } cluster_dict_list.append(cluster_dict) return cluster_dict_list
[docs] def magnetically_equivalent_spins( molecule_parameters: datatypes.NMRParameters, calculation_parameters: datatypes.NMRCalculationParameters, tolerance_couplings: float = 1.0, tolerance_shifts: float = 1.0, verbose: int = 0, ) -> list[np.ndarray]: """Identify magnetically equivalent spins and return them as groups. Args: molecule_parameters: The molecular isotopes, shifts and J-coupling values. calculation_parameters: Object storing all parameters specifying how a calculation will run. tolerance_couplings: Tolerance for the J-couplings in the group identifier in percent. Defaults to 1. tolerance_shifts: Tolerance for the shifts in the group identifier in percent. Defaults to 1. verbose: Verbosity level of output. Defaults to 0. Returns: A list of integer arrays with the indices of spins in the individual groups. """ hamiltonian = nmr_hamiltonian( molecule_parameters, calculation_parameters.field_T, reference_isotope=calculation_parameters.reference_isotope, gyromagnetic_ratios=calculation_parameters.gyromagnetic_ratios, ) Jz, hJz, _ = extract_coupling(hamiltonian) groups_in_partition = identify_symmetry_groups( hJz, Jz, tolerance_couplings_percent=tolerance_couplings, tolerance_shifts_percent=tolerance_shifts, verbose=verbose, ) return groups_in_partition
[docs] def normalized_spin_dependent_gyromagnetic_ratios( molecule_parameters: datatypes.NMRParameters, gyromagnetic_ratios: dict[Isotope, float] = GyromagneticRatios, ) -> np.ndarray: """Gyromagnetic ratios of each spin, normalized w.r.t. the gyromagnetic ratio of hydrogen. Args: molecule_parameters: The molecular isotopes, shifts and J-coupling values. gyromagnetic_ratios: Dictionary of gyromagnetic ratios in rad / (T s). Returns: Scaled gyromagnetic ratios of each spin. """ gammas: np.ndarray = np.zeros(len(molecule_parameters.shifts)) for i, isotope in enumerate(molecule_parameters.isotopes): gammas[i] = gyromagnetic_ratios[isotope] / gyromagnetic_ratios[Isotope(1, "H")] return gammas