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 thechemical_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 inchemical_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 theplot_x
keywordmax_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 theplot_x
keywordplot_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 basispos_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