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.

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(cluster_size: int, linewidth: float = 4) → Tuple[numpy.array, numpy.array]¶

Simulate the NMR spectrum of the spin system represented by self

We specify a cluster_size that specfies the maximum size of clusters to solve exactly. 4 seems to be a good number where the simulation is not too onerous and the accuracy is good. 1 will turn off all couplings of spins

Parameters

cluster_size – the maximum size of a given 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(cluster_size: int) → scipy.sparse._csr.csr_matrix¶

Create the clusters to be solved exactly by breaking the system into all (potentially overlapping) sets of clusters of size at most cluster_size.

Parameters: cluster_size – max size of clusters to fragment system into

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._csc.csc_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._csc.csc_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._coo.coo_matrix, scipy.sparse._coo.coo_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

schrodinger.application.jaguar.nmr_spectrum_sim.expv_mat(t: float, A: scipy.sparse._csr.csr_matrix, v: scipy.sparse._csr.csr_matrix) → scipy.sparse._csr.csr_matrix¶

Driver function for expv, allows for multi-vector inputs.

Maybe there is a way that this and expv could be combined to do more work simultaneously and be faster.

Parameters

t – scaling factor in exp(t*A)*v
A – matrix to exponentiate in exp(t*A)*v
v – matrix to contract with matrix exponential in exp(t*A)*v

schrodinger.application.jaguar.nmr_spectrum_sim.expv(t: float, A: scipy.sparse._csr.csr_matrix, v: scipy.sparse._csr.csr_matrix, tol: Optional[float] = None, m: Optional[int] = None) → scipy.sparse._csr.csr_matrix¶

A Krylov space method of computing the expm_multiply, i.e. exp(t*A)*v

Only works on vectors (as opposed to matrices). Was much slower than the scipy.sparse.linalg method.

Leaving this in case it is of use to someone else.

Parameters

t – scaling factor in exp(t*A)*v
A – matrix to exponentiate in exp(t*A)*v
v – vector to contract with matrix exponential in exp(t*A)*v, in the form of n x 1 sparse matrix
tol – tolerance to converge to
m – size of Krylov space

schrodinger.application.jaguar.nmr_spectrum_sim.better_kron(A: scipy.sparse._coo.coo_matrix, B: scipy.sparse._coo.coo_matrix) → scipy.sparse._coo.coo_matrix¶

Kronecker product of sparse matrices A and B assuming that input and output are coo_matrices

COO matrices are good for creating sparse matrices efficiently. Since we are doing repeated, nested Kronecker products, it is much more efficient to not convert between matrix types before and after each operation. This optimization works much better for our case.

Parameters

A – first matrix of the product
B – second matrix of the product

Returns

kronecker product in a sparse coo_matrix format