Overview of the nmrsim Top-Level API

This notebook gives a tour of the top level classes the nmrsim API provides. These are conveniences that abstract away lower-level API functions. Users wanting more control can consult the full API documentation.

[1]:
import os
import sys
import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.dpi']= 300
%matplotlib inline
[2]:
%config InlineBackend.figure_format = 'svg'  # makes inline plot look less blurry
[3]:
home_path = os.path.abspath(os.path.join('..', '..', '..'))
if home_path not in sys.path:
    sys.path.append(home_path)

tests_path = os.path.abspath(os.path.join('..', '..', '..', 'tests'))
if tests_path not in sys.path:
    sys.path.append(tests_path)

Definitions

In naming classes, functions, methods, data types etc. certain phrases, taken from NMR nomenclature, have the following interpretations:

  • multiplet (e.g. the nmrsim.Multiplet class): a first-order simulation for one signal (i.e. one or more chemical shift-equivalent nuclei). Examples: doublet, triplet, doublet of triplets, but not an AB quartet (which is a second-order pattern for two nuclei).
  • spin system (e.g. the SpinSystem class): a simulation of a set of coupled nuclei.
  • spectrum (e.g. the Spectrum class): a complete collection of first- and/or second-order components for simulating a total NMR spectrum. ‘Spectrum’ can also refer in general to the simulation results for the system, e.g a peaklist or lineshape (see below).
  • peak: a pair of frequency (Hz), intensity values corresponding to a resonance in an NMR spectrum. For example, a 1H triplet centered at 100 Hz with J = 10 Hz would have the following peaks: (110, 0.25), (100, 0.5), (90, 0.25).
  • peaklist: a list of peaks (e.g. [(110, 0.25), (100, 0.5), (90, 0.25)] for the above triplet).
  • lineshape: a pair of [x_coordinates…], [y_coordinates] arrays for plotting the lineshape of a spectrum.

In this notebook the term list is interchangeable with other iterables such as numpy arrays or tuples. As much as possible, nmrsim relies on <”duck typing”>(https://en.wikipedia.org/wiki/Duck_typing) to accept a variety of iterables as inputs, converting them to specific types such as numpy arrays as needed. The term matrix refers to a 2D array-like object in general, e.g. a list of lists or a 2D numpy array. It does not refer specifically to the (marked-for-deprecation) numpy.matrix class.

The following idioms are used for arguments: * v for a frequency or list of frequencies (similar to \(\nu\) ). * I for a signal intensity * J for coupling constant data (exact format depends on the implementation).

Scenario: user wants to plot a spectrum for an ABX 3-spin system.

A spin system can be described using a list of frequencies v and J (coupling constant) data . For this example, a function from nmrsim’s test suite will provide some example data:

[4]:
# This dataset is for the vinyl group of vinyl acetate, as used in:
# http://www.users.csbsju.edu/~frioux/nmr/ABC-NMR-Tensor.pdf
def rioux():
    v = np.array([430.0, 265.0, 300.0])
    J = np.zeros((3, 3))
    J[0, 1] = 7.0
    J[0, 2] = 15.0
    J[1, 2] = 1.50
    J = J + J.T
    return v, J
[5]:
v, J = rioux()
print('v: ', v)  # frequencies in Hz
print('J: \n', J)  # matrix of coupling constants
v:  [430. 265. 300.]
J:
 [[ 0.   7.  15. ]
 [ 7.   0.   1.5]
 [15.   1.5  0. ]]

The J matrix is constructed so that J[a, b] is the coupling constant between v[a] and v[b]. The diagonal elements should be 0.

The SpinSystem class can be used to model a set of coupled nuclei.

[6]:
from nmrsim import SpinSystem
[7]:
abx_system = SpinSystem(v, J)
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

The SpinSystem.peaklist() method returns the peaklist for the simulation:

[8]:
abx_system.peaklist()
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
[8]:
[(260.6615285748296, 0.23011249131787795),
 (291.3191136690316, 0.22882003310401866),
 (419.5193577561387, 0.29107244545559474),
 (292.8468885409388, 0.2138123115725174),
 (426.4877446901904, 0.26629867696733883),
 (262.18930344673686, 0.24876061726413431),
 (434.5231959501799, 0.2300458619680737),
 (267.62991550888137, 0.24855578963215708),
 (306.32295186307283, 0.29251680284079634),
 (441.49158288423155, 0.21257278181929304),
 (307.85072673497996, 0.2648680204713065),
 (269.1576903807885, 0.27256416758689195)]

You can plot this data with the visualization library of your choice. However, the nmrsim.plt library has functions for convenient plotting of common nmrsim data types. The plt.mplplot function will take a peaklist and use matplotlib to plot the corresponding lineshape. The optional keyword argument y_max can be used to set the maximum for the y-axis (and y_min for the minimum).

[9]:
from nmrsim.plt import mplplot
[10]:
mplplot(abx_system.peaklist(), y_max=0.2);
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
[<matplotlib.lines.Line2D object at 0x7efdf6ce00d0>]
../_images/jupyter_1-Introduction-to-API_18_1.svg

To plot the spectra as a “stick” style plot (single lines for each peak, rather than a simulated lineshape), you can use the mplplot_stick function instead of mplplot:

[11]:
from nmrsim.plt import mplplot_stick
[12]:
# The range of the x axis can be specified using the 'limits' keyword argument:
mplplot_stick(abx_system.peaklist(), y_max=0.3, limits=(250, 320));
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
../_images/jupyter_1-Introduction-to-API_21_1.svg

SpinSystem defaults to second-order simulation of a spin system. If the SpinSystem object is instantiated with the second_order=False keyword argument, or if the SpinSystem.second_order attribute is set to False, first-order simulation will be performed instead.

[13]:
abx_system.second_order = False
mplplot(abx_system.peaklist(), y_max=0.2);
[<matplotlib.lines.Line2D object at 0x7efdf714ba90>]
../_images/jupyter_1-Introduction-to-API_23_1.svg

Depending on the resolution of the plot and how the data points for the lineshape are interpolated, the peak heights may not look identical. The correct relative intensities can be seen in the stick plot, however:

[14]:
mplplot_stick(abx_system.peaklist(), y_max=0.3);
../_images/jupyter_1-Introduction-to-API_25_0.svg

Scenario: User wants to simulate individual first-order multiplets

The Multiplet class can be used to represent an individual first-order multiplet.

[15]:
from nmrsim import Multiplet

Required arguments for Multiplet are the central frequency v, the intensity I (“integration”) in the absence of coupling, and a list of coupling data J. These arguments become attributes of Multiplet. Each list entry is a tuple of (J value in Hz, number of nuclei causing the coupling). For example, the following Multiplet represents: 1200 Hz, 2H, td, J = 7.1, 1.1 Hz.

[16]:
# 1200 Hz, 2H, td, J= 7.1, 1.1 Hz
td = Multiplet(1200.0, 2, [(7.1, 2), (1.1, 1)])
print(td.v)
print(td.I)
print(td.J)
1200.0
2
[(7.1, 2), (1.1, 1)]

The Multiplet.peaklist() method returns the peaklist for the multiplet:

[17]:
mplplot_stick(td.peaklist());
../_images/jupyter_1-Introduction-to-API_32_0.svg
[18]:
mplplot(td.peaklist());
[<matplotlib.lines.Line2D object at 0x7efdf70e51d0>]
../_images/jupyter_1-Introduction-to-API_33_1.svg

Multiplet attributes can be modified.

[19]:
td2 = Multiplet(1200.0, 2, [(7.1, 2), (1.1, 1)])
td2.v = 1100
mplplot(td2.peaklist());
[<matplotlib.lines.Line2D object at 0x7efdf7014bd0>]
../_images/jupyter_1-Introduction-to-API_35_1.svg

If a Multiplet is multiplied by a scalar, a new Multiplet is returned that has all intensities multiplied by the scalar. In-place multiplication (*=) modifies the original Multiplet object.

[20]:
td3 = td2 * 2
td2 *= 2
assert td2 is not td3
mplplot(td2.peaklist());
[<matplotlib.lines.Line2D object at 0x7efdf7002350>]
../_images/jupyter_1-Introduction-to-API_37_1.svg

Multiplets are equal to each other if their peaklists are equal.

[21]:
assert td2 == td3

Division and division in place is also possible:

[22]:
td4 = td2 / 2
td2 /= 2
assert td4 == td2

If two multiplets are added together, the result is a Spectrum object. See the next Scenario for the usage of Spectrum.

Scenario: User wants to simulate a spectrum built from individual components

Any object that has a .peaklist() method can be used to create a Spectrum object.

A Spectrum object can be specifically created by providing a list of components as the first argument:

[23]:
from nmrsim import Spectrum
[24]:
two_td = Spectrum([td, td3])
[25]:
mplplot(two_td.peaklist());
[<matplotlib.lines.Line2D object at 0x7efdf72ee2d0>]
../_images/jupyter_1-Introduction-to-API_47_1.svg

A Spectrum object is also returned from certain binary operations, such as addition:

[26]:
td3.v = 1000
td4.v = 900

all_tds = td + td2 + td3 + td4
mplplot(all_tds.peaklist());
[<matplotlib.lines.Line2D object at 0x7efdf6ca2890>]
../_images/jupyter_1-Introduction-to-API_49_1.svg

A Spectrum can be composed from both first- and second-order components:

[27]:
combo_spectrum = abx_system + td3 + td4

# mplplot has an optional y_max keyword argument to set the max range of the y-axis
mplplot(combo_spectrum.peaklist(), y_max=0.4);
[<matplotlib.lines.Line2D object at 0x7efdf74afd50>]
../_images/jupyter_1-Introduction-to-API_51_1.svg

Scenario: User wants to model a specific spin system using an explicit (non-qm) solution

The nmrsim.partial module contains “canned” mathematical solutions for second-order systems.

Example: simulate the AB part of an ABX3 system

[28]:
from nmrsim.discrete import ABX3
[29]:
help(ABX3)
Help on function ABX3 in module nmrsim.discrete:

ABX3(Jab, Jax, Jbx, Vab, Vcentr)
    Simulation of the AB part of an ABX3 spin system.

    Parameters
    ---------
    Jab : float
        the Ha-Hb coupling constant (Hz).
    Jax : float
        the Ha-Hb coupling constant (Hz).
    Jbx : float
        the Ha-Hb coupling constant (Hz).
    Vab : float
        the difference in the frequencies (Hz) of Ha and Hb in the absence of
        coupling. Positive when vb > va.
    Vcentr : float
        the frequency (Hz) for the center of the AB signal.

    Returns
    -------
    [(float, float)...]
        a list of (frequency, intensity) tuples.

[30]:
abx3_peaklist = ABX3(-12, 7, 7, 14, 150)
mplplot(abx3_peaklist, y_max=0.25);
[<matplotlib.lines.Line2D object at 0x7efdf7396490>]
../_images/jupyter_1-Introduction-to-API_56_1.svg

Here is an alternate, non-qm simulation for the ABX system from the SpinSystem demonstration:

[31]:
from nmrsim.discrete import ABX
[32]:
help(ABX)
Help on function ABX in module nmrsim.discrete:

ABX(Jab, Jax, Jbx, Vab, Vcentr, vx, normalize=True)
    Non-QM approximation for an ABX spin system. The approximation assumes
    that Hx is very far away in chemical shift from Ha/Hb.

    Parameters
    ---------
    Jab : float
        The Ha-Hb coupling constant (Hz).
    Jax : float
        The Ha-Hx coupling constant (Hz).
    Jbx : float
        The Hb-Hx coupling constant (Hz).
    Vab : float
        The difference in the frequencies (in the absence of
        coupling) of Ha and Hb (Hz).
    Vcentr : float
        The frequency (Hz) for the center of the AB signal.
    vx : float
        The frequency (Hz) for Hx in the absence of coupling.

    normalize: bool (optional)
        whether the signal intensity should be normalized. If false, the total
        signal intensity happens to be ~12.

    Returns
    -------
    [(float, float)...]
        a list of (frequency, intensity) tuples.

[33]:
abx_peaklist = ABX(1.5, 7, 15, 35, 282.5, 430)
mplplot(abx_peaklist, y_max=0.4);
[<matplotlib.lines.Line2D object at 0x7efdf6ee2610>]
../_images/jupyter_1-Introduction-to-API_60_1.svg

Scenario: User wants to model DNMR two-spin exchange, without and with coupling

The nmrsim.dnmr library provides functions for calculating DNMR lineshapes, and classes to describe these systems. Currently, models for two uncoupled nuclei and two coupled nuclei are provided.

[34]:
from nmrsim.dnmr import DnmrTwoSinglets, DnmrAB

For: va = 165 Hz, vb = 135 Hz, k = 65.9 s-1, line widths (at the slow exchange limit) wa and wb = 0.5 Hz, and population of state a = 0.5 (i.e. 50%):

[35]:
two_singlet_system = DnmrTwoSinglets(165.00, 135.00, 65.9, 0.50, 0.50, 0.50)
[36]:
from nmrsim.plt import mplplot_lineshape
[37]:
mplplot_lineshape(*two_singlet_system.lineshape());
../_images/jupyter_1-Introduction-to-API_67_0.svg

Class attributes can be changed. In the previous case, k = 65.9 -1 corresponds to the point of coalescence. When the rate of exchange is lower, two separate peaks are observed.

[38]:
two_singlet_system.k = 5
[39]:
mplplot_lineshape(*two_singlet_system.lineshape());
../_images/jupyter_1-Introduction-to-API_70_0.svg

What if the relative populations of states a and b are 75% and 25%, respectively?

[40]:
two_singlet_system.pa = 0.75
mplplot_lineshape(*two_singlet_system.lineshape());
../_images/jupyter_1-Introduction-to-API_72_0.svg

To model an AB-like system of two coupled nuclei undergoing exchange, use the DnmrAB class. In the following example, the frequencies are the same as for the previous system. J = 5 Hz, k = 12 -1, and the line width (at the slow exchange limit) is 0.5 Hz.

[41]:
from nmrsim.dnmr import DnmrAB
[42]:
AB = DnmrAB(165, 135, 5, 10, 0.5)
[43]:
mplplot_lineshape(*AB.lineshape());
../_images/jupyter_1-Introduction-to-API_76_0.svg