schrodinger.application.jaguar.nmr_spectrum_sim module

Simulate NMR spectrum from chemical shifts and coupling constants.

This module implements the method of Kuprov, Wagner-Rundell, and Hore: I. Kuprov, Polynomially scaling spin dynamics simulation algorithm based on adaptive state-space restriction, J.Magn. Reson. 189 (2007) 241–250 DOI: 10.1016/j.jmr.2007.09.014

Basically, the system is broken up into overlapping clusters of spins. We then simulate the actual NMR experiment, preparing a state with all spins aligned on the Z-axis, knocking them all into the xy-plane with a pi/2 pulse, and allowing them to relax back to the Z-axis, sampling along the way to obtain the time-domain version of the NMR spectrum, which we then Fourier transform to the usual frequency-domain spectrum. We do it in this way as it allows a natural way to account for the effects of the overlapping cluster, without having to diagonalize a 2^N x 2^N matrix where N is the total number of spins in the molecule. We also still account for second order effects in the spectrum.

Copyright Schrodinger, LLC. All rights reserved.

schrodinger.application.jaguar.nmr_spectrum_sim.PauliMatrices

alias of schrodinger.application.jaguar.nmr_spectrum_sim.PauliMats

schrodinger.application.jaguar.nmr_spectrum_sim.build_pauli(mult: int) schrodinger.application.jaguar.nmr_spectrum_sim.PauliMats

Build Pauli matrices for the given multiplicity.

We use these so often it’s worth it to cache these for the two values we might concievably allow (mult 2 and 3).

Parameters

mult – multiplicity to create the Pauli matrices of

Returns

namedtuple with the x, y, z, plus, and minus Pauli matrices

class schrodinger.application.jaguar.nmr_spectrum_sim.SpinSystem(chemical_shifts: List[float], J_couplings: List[Tuple[int, int, float]], field: float, addl_couplings: Iterable[Tuple[int, int, float]] = (), min_ppm: float = - 0.5, max_ppm: float = 9.0)

Bases: object

SPEC_BUFFER = 0.9
__init__(chemical_shifts: List[float], J_couplings: List[Tuple[int, int, float]], field: float, addl_couplings: Iterable[Tuple[int, int, float]] = (), min_ppm: float = - 0.5, max_ppm: float = 9.0)

An object storing a spin system which can be used to simulate an NMR spectrum

Parameters
  • chemical_shifts – list of chemical shifts for spins in system (in ppm)

  • J_couplings – list of J-coupling (in Hz) tuples, the first two values are indices in the chemical_shifts list of which are coupled with a value of the third value (i.e. (2, 3, 4.5) indicates that the third and fourth spins in the chemical_shifts list are coupled with a J-coupling value of 4.5 Hz)

  • field – magnetic field strength (used to convert between ppm and Hz)

  • addl_couplings – list of J-coupling (in Hz) to spins that should not be displayed, the first value is the real, displayed spin in chemical_shifts which is coupled to a phantom spin (the second value) by the third value (i.e. (1, 0, 10.1) indicates that the second spin in chemical_shifts is coupled to the first undisplayed spin with a J-coupling value of 10.1 Hz). This scenario may arise from 31P or 19F splitting of a 1H or 13C spectrum.

  • min_ppm – This is the min ppm to consider for a real displayed spin. If a spin is specified in chemical_shifts outside of this window, the chemical_shifts list will be deferred to. To adjust the plot window, use the plot_x keyword

  • max_ppm – This is the max ppm to consider for a real displayed spin. If a spin is specified in chemical_shifts outside of this window, the chemical_shifts list will be deferred to. To adjust the plot window, use the plot_x keyword

  • plot_x – Specifies the bounds of the plot window in ppm (larger ppm first)

  • plot_y – Specifies the bounds of the plot window in arbitrary intensity units

simulate_spectrum(nth_neighbors: int, linewidth: float = 4) Tuple[numpy.array, numpy.array]

Simulate the NMR spectrum of the spin system represented by self

We specify a nth_neighbors that describes, regarding the spin system as a graph where spins are nodes and edges are between spins that have non-zero J-coupling, how many nearest neighbors to include in clusters that are treated exactly. 4 seems to be a good number where the simulation is not too onerous and the accuracy is good. 0 will turn off all couplings of spins

Parameters
  • nth_neighbors – the nth_neighbors-th nearest neighbors are in one cluster, must be a non-negative integer

  • linewidth – approximate linewidth (in Hz) for spectra, 4 Hz is a good value for 1H spectra and 40Hz is good for 13C

make_clusters(nth_neighbors: int) scipy.sparse._csr.csr_matrix

Create the clusters to be solved exactly by breaking the system into potentially overlapping sets of n nearest neighbors (regarding the spins as a graph where there is an edge present if there is a J-coupling between those spins)

Parameters

nth_neighbors – number of nearest neighbors to include in the clusters, 0 = ignore all couplings, 1 = clusters consist only of nearest neighbors, 2 = first and second nearest neighbors, and so on

make_state_lists() scipy.sparse._csr.csr_matrix

Make the state lists. These are just the combination of all possible states (numbered 0 through multipl^2) scatterized to the right place in a big sparse matrix.

That is, if you had 3 spins with multiplicity 2 (i.e. spin 1/2) and the first, and third were in one cluster and the second was in another, we are generating:

000 001 002 003 100 101 102 103 200 201 … 302 303 010 020 030

Note, the uniquification step. It’s possible for one spin to be in more than one cluster, so, for example, the state where that spin is in state 1 and all other spins are state 0 appears more than once in the concatenation of lists generated for each cluster.

Returns

sparse matrix where each row is a state that are treated exactly in this method (using clustering) and each column is the state of any given spin in the system in the system-level states.

make_cross_membership_matrix() Tuple[numpy.array]

Build the cross-membership matrix which connects clusters and the states that involve them.

Specifically, the cross-membership matrix is 1 for (i, j) if state i involves cluster j.

For our purpose, we don’t need the whole matrix; all we need is essentially a dictionary keyed by the cluster number whose values are the states that cluster is involved in.

Returns

map from cluster index to array of state indices that involve that cluster

compute_density_matrices() List[Dict[int, scipy.sparse._csr.csr_matrix]]

Compute cluster density matrices for each state in the cluster.

We store these as a list of dicts because the first index is the cluster which will have values from 0 to n_clusters, while the second index is the state index which is generally only a subset of the ints from 0 to n_states.

Returns

2D map from cluster, state to density matrix (e.g. return[cluster_idx][state_idx] is the density matrix corresponding to cluster_idx and state_idx)

convert_state_to_density_mat(state_number: int, cluster_number: int) scipy.sparse._csr.csr_matrix

Convert a given state to a density matrix.

Parameters
  • state_number – state to convert to density matrix

  • cluster_number – cluster in which to write density matrix

Returns

the density matrix for the given state

compute_cluster_momentum_operators() List[scipy.sparse._csr.csr_matrix]

Compute cluser representation of the overall transverse angular momentum on y-axis.

For 1D NMR, we only need one these and so arbitrarily we do Ly. If 2D NMR is someday desired, cluster Lx can be created by swapping L.x for L.y in the expressions below.

Returns

map from cluster index to its transverse angular momentum operator

get_mats_for_cluster_placement(pos: int, pos_vector: List[int]) Tuple[scipy.sparse._dia.dia_matrix, scipy.sparse._dia.dia_matrix]

Create the matrices that can be used to relate a given spin’s Pauli-sized matrix to the locations in the entire cluster which that spin affects.

Parameters
  • pos – position within the pos_vector to expand into the cluster basis

  • pos_vector – vector of spins which are involved in the given cluster

Returns

a pair of sparse identity matrices which indicates the locations in the full cluster basis that are affected by the pos-th spin in the cluster

build_cluster_hamilitonians() List[scipy.sparse._csr.csr_matrix]

Build the cluster Hamiltonians.

Returns

list of Hamiltonians for each cluster

build_hamil_for_cluster(cluster_number: int) scipy.sparse._csr.csr_matrix

Build the Hamiltonian for the given cluster.

Parameters

cluster_number – index of cluster to build Hamiltonian of

Returns

Hamiltonian matrix for this cluster

make_liouvillians() Tuple[scipy.sparse._csr.csr_matrix, scipy.sparse._csr.csr_matrix]

Build the Liouvillians for this spin system.

For 1D NMR, we only need, Ly. If 2D is need in the future, Lx can also be generated analogously (although without the -j prefactor).

Returns

the Liouvillians of the system under the evolution of the Hamiltonian and under a transverse angular momentum perturbation

get_Splus_Sz() Tuple[scipy.sparse._csr.csr_matrix, scipy.sparse._csr.csr_matrix]

Determine the states that look like S+ states and Sz states.

Since we numbered Sz as 1 and S+ as 2 (see convert_state_to_dm), this amounts to finding the indices of the the rows in the state list with a single entry of 1 or 2, respectively.

We store the collection of these states as a sparse n_states x 1 matrix (i.e. a sparse vector).

Returns

the Sp_pure_states vector is the detection state vector, the Sz_pure_states vector is the initial state vector (i.e. aligned along the magnetic field, which by convention is the z-axis)

build_propagators(timestep: float)

Compute the system propagator and the pi/2 pulse matrix

Parameters

timestep – time step to use for propagation of the state

propagate_system(linewidth)

Generate spectrum trace by simulating system evolution.

We generate an FID (free induction decay) of the system initially aligned along the z-axis after a pi/2 pulse knocking the system into the xy-plane and allowing it to evolve back to the z-axis. An FFT of this FID is the NMR spectrum.

Parameters

linewidth – approximate linewidth (in Hz) for spectra, 4 Hz is a good value for 1H spectra and 40Hz is good for 13C

Returns

x-coords and y-coords of NMR spectrum (x-coords in ppm, y-coords in arbitrary intensity units)

schrodinger.application.jaguar.nmr_spectrum_sim.uniquify_rows(mat: scipy.sparse._csr.csr_matrix, remove_zero_row: bool = False) scipy.sparse._csr.csr_matrix

Uniquify the rows of a sparse matrix, retaining the order where each row is first encountered

Parameters
  • mat – sparse matrix with rows that are possibly duplicates of each other

  • remove_zero_row – If True, remove a row if it is all zeros

Returns

sparse matrix with duplicate rows removed

schrodinger.application.jaguar.nmr_spectrum_sim.remove_subset_clusters(mat: scipy.sparse._csr.csr_matrix) scipy.sparse._csr.csr_matrix

Remove clusters that are subsets of other clusters.

If a cluster is purely a subset of another cluster, it does not add any value to the simulation.

Parameters

mat – matrix with rows that may be subsets

Returns

matrix with subset rows removed

schrodinger.application.jaguar.nmr_spectrum_sim.scatterize_sparse_mat(target_mat: scipy.sparse._csr.csr_matrix, target_indices: List[int], src_sp_mat: scipy.sparse._csr.csr_matrix)

Scatter the data from one sparse matrix to another using a uniform map from the source matrix to the target spase matrix.

Parameters
  • target_mat – Sparse matrix to be populated with data from the source matrix src_sp_mat

  • target_indices – map from source matrix indices to target matrix indices. In our case, we are uniformly tranforming from both axes by the same map.

  • src_sp_mat – Sparse matrix whose data will be copied to target_mat

schrodinger.application.jaguar.nmr_spectrum_sim.build_cluster_liouvillian(mat: scipy.sparse._csr.csr_matrix, state_vectors: scipy.sparse._csr.csr_matrix) scipy.sparse._csr.csr_matrix

Convert a matrix to the Liouvillian representation.

In our case, this is the Liouvillian of each cluster, rather than the whole system.

Parameters
  • mat – matrix to be converted to a Liouvillian representation.

  • state_vector – matrix of state vectors which form the basis of the Liouvillian representation

Returns

the matrix in the Liouvllian representation.

schrodinger.application.jaguar.nmr_spectrum_sim.remove_vals_below_thresh(mat: scipy.sparse._csr.csr_matrix, thresh: float)

Floor small values to 0 in sparse matrices so that they are not stored explicitly.

Parameters
  • mat – matrix possibly with some small values

  • thresh – threshold for small values to be floored to zeros