Method¶
While mrsimulator’s organization of the SpinSystem object and its composite objects, Site and Coupling, are easily understood by anyone familiar with the underlying physical concepts, the organization of the Method object in mrsimulator and its related composite objects require a more detailed explanation of their design. This section assumes that you are already familiar with the topics covered in the Introduction sections Getting Started, Isotopomers Example, and LeastSquares Fitting Example.
Note
Before writing your own custom Method, check if any of our prebuilt methods in the Methods Library can serve your needs.
Overview¶
An experimental NMR method involves a sequence of rf pulses, free evolution periods, and sample motion. The Method object in mrsimulator models the spectrum from an NMR pulse sequence. The Method object is designed to be versatile in its ability to model spectra from various multipulse NMR methods using concepts from the symmetry pathway approach where a pulse sequence is understood in terms of a set of desired (and undesired) transition pathways. Each transition pathway is associated with a single resonance in a multidimensional NMR spectrum. The transition pathway signal encodes information about the spin system interactions in its amplitude and correlated frequencies. Consider the illustration of a 2D pulse sequence shown below, where a desired signal for the method is associated with a particular transition pathway, \({\hat{A} \rightarrow \hat{B} \rightarrow \hat{C} \rightarrow \hat{D} \rightarrow \hat {E} \rightarrow \hat{F}}\).
Here, the first spectral dimension, i.e., the Fourier transform of the transition pathway signal as a function of \(t_1\), derives its average frequency, \(\overline{\Omega}_1\), from a weighted average of the \(\hat{A}\), \(\hat{B}\), and \(\hat{C}\) transition frequencies. The second spectral dimension, i.e., the FT with respect to \(t_2\), derives its average frequency, \(\overline{\Omega}_2\), from a weighted average of the \(\hat{E}\), and \(\hat{F}\) transition frequencies. Much of the experimental design and implementation of an NMR method is in identifying the desired transition pathways and finding ways to acquire their signals while eliminating all undesired transition pathway signals.
While NMR measurements occur in the time domain, mrsimulator simulates the corresponding multidimensional spectra directly in the frequency domain. The Method object in mrsimulator needs only a few details of the NMR pulse sequence to generate the spectrum. It mimics the result of the pulse sequence given the desired transition pathways and their complex amplitudes and average frequencies in each spectroscopic dimension of the dataset. To this end, a Method object is organized according to the UML diagram below.
Note
In UML (Unified Modeling Language) diagrams, each class is represented with a box that contains two compartments. The top compartment has the class’s name, and the bottom compartment contains the class’s attributes. Default attribute values are shown as assignments. A composition is depicted as a binary association decorated with a filled black diamond. Inheritance is shown as a line with a hollow triangle as an arrowhead.
At the heart of a Method object, assigned to its attribute
spectral_dimensions
, is an ordered list of SpectralDimension objects
in the same order as the time evolution dimensions of the experimental NMR
sequence. In each SpectralDimension object is an ordered list of Events
objects assigned to the attribute events
; Event objects are divided
into three types: (1) SpectralEvent()
, (2)
DelayEvent()
, and (3)
MixingEvent()
. This ordered list of Event objects
is used to select the desired transition pathways and determine their average
frequency and complex amplitude in the SpectralDimension.
Warning
DelayEvent objects are not available in version 0.7 of mrsimulator.
SpectralEvent and DelayEvent objects define which transitions are
observed during the event and under which transitiondependent frequency
contributions they evolve. No coherence transfer among transitions or
populations occurs in a spectral or delay event. The transitiondependent
frequency contributions during an Event are selected from a list of
enumeration literals and placed in the freq_contrib
attribute of the event. If freq_contrib
is left unspecified, i.e., the
value of freq_contrib
is set to None
, a default list holding the
enumeration literals for all contributions is generated for the event.
Note
All frequency contributions from direct and indirect spinspin couplings are calculated in the weakcoupling limit in mrsimulator.
Additionally, the user can affect transition frequencies during a spectral or
delay event by changing other measurement attributes: rotor_frequency
,
rotor_angle
, and magnetic_flux_density
. If left unspecified, these
attributes default to the values of the identically named global attributes in
the Method object. SpectralEvent objects use the fraction
attribute to
calculate the weighted average frequency during the spectral dimension for each
selected transition pathway.
Inside SpectralEvent and DelayEvent objects, is a list of
TransitionQuery()
objects (vide infra)
which determine which transitions are observed during the event. Method
objects in mrsimulator are generalpurpose because they are designed for an
arbitrary spin system, i.e., a method does not know the spin system in advance.
When designing a Method object, you cannot identify and select a transition
through its initial and final eigenstate quantum numbers. Transition selection
is done through TransitionQuery and
SymmetryQuery()
objects during individual
spectral or delay events. TransitionQuery objects can hold a
SymmetryQuery object in the attributes ch1
, ch2
, or ch3
, which act on
specific isotopes defined by the channels
attribute in Method. It is
only during a simulation that the Method object uses its TransitionQuery
objects to determine the selected transition pathways for a given SpinSystem
object by the initial and final eigenstate quantum numbers of each transition.
Between adjacent SpectralEvent or DelayEvent objects, mrsimulator defaults
to total mixing, i.e., connecting all selected transitions in the two adjacent
spectral or delay events. This default behavior can be overridden by placing an
explicit MixingEvent object between such events. Inside MixingEvent
objects is a MixingQuery()
object, which
determines the coherence transfer amplitude between transitions. A
MixingQuery object holds
RotationQuery()
objects acting on specific
isotopes in the spin system. As before, the isotope upon which the
RotationQuery objects act is determined by the channels
attribute in the
Method object.
In this guide to designing custom Method objects, we begin with a brief review of the relevant Symmetry Pathway concepts employed in mrsimulator. This review is necessary for understanding (1) how transitions are selected during spectral and delay events and (2) how average signal frequencies and amplitudes in each spectral dimension are determined. We outline the procedures for designing and creating TransitionQuery and MixingQuery objects for single and multispin transitions and how to use them to select the transition pathways with the desired frequency and amplitudes in each spectral dimension of your custom Method object. In multidimensional spectra, we illustrate how the desired frequency correlation can sometimes be achieved by using an appropriate affine transformation. We also examine how changing the frequency contributions in SpectralEvent of DelayEvent objects can be used to obtain the desired frequency and amplitude behavior. The ability to select frequency contributions can often reduce the number of events needed in the design of your custom Method object.
Theoretical Background¶
Before giving details on how to create a custom Method object, we review a few key concepts about spin transitions and transition symmetry functions.
The number of quantized energy eigenstates for \(N\) coupled nuclei is
where \(I_u\) is the total spin angular momentum quantum number of the \(u\text{th}\) nucleus. The transition from quantized energy level \(E_i\) to \(E_j\) is one of
possible transitions between the \(\Upsilon\) levels. Here we count \(i \rightarrow j\) and \(j \rightarrow i\) as different transitions. For example, a single spin with angular momentum quantum number \(I=3/2\) will have \(\Upsilon_{\left\{ 3/2 \right\}} = 2I+1 = 4\) energy levels and \(\mathcal{N}_{\left\{ 3/2 \right\}} = 2I(2I+1) = 12\) possible NMR transitions. A two spin system, with quantum numbers \(I = 1/2\) and \(S = 1/2\), will have
energy levels and
possible NMR transitions. We write a transition (coherence) from state \(i\) to \(j\) using the outer product notation \(\ketbra{j}{i}\). In mrsimulator, all simulations are performed in the highfield limit and further, assume that all spinspin couplings are in the weak limit.
To write a custom Method in mrsimulator, you’ll need to determine the desired transition pathways and select the desired transitions during each SpectralEvent or DelayEvent. Keep in mind, however, that Method objects are designed without any details of the spin systems upon which they will act. For example, in the density matrix of a spin system ensemble, one could easily identify a transition by its row and column indexes. However, those indexes depend on the spin system and how the spins and their eigenstates have been assigned to those indexes. Instead, we need a spinsystem agnostic approach for selecting transitions.
Spin Transition Symmetry Functions¶
One way you can select a subset of singlespin transitions if you don’t know the energy eigenstate quantum numbers is to request all transitions whose singlespin transition symmetry function, \(\text{p}_I\) is \(1\), i.e.,
The \(\text{p}_I\) singlespin transition symmetry function is also known as the singlespin coherence order of the transition.
Note
In the high field limit, only singlespin transitions with \({\text{p}_I = \pm 1}\) are directly observed. Since the evolution frequencies of the \(\ketbra{j}{i}\) and \(\ketbra{i}{j}\) transitions are equal in magnitude but opposite in sign, the convention is to only present the \({\text{p}_I =  1}\) transition resonances in singlequantum spectra.
By selecting only singlespin transitions with \(\text{p}_I = 1\), you get all the “observed” transitions from the set of all possible transitions. Similarly, you can use \(\text{p}_I\) to select any subset of singlespin transitions, such as doublequantum \((\text{p}_I = \pm 2)\) transitions, triplequantum \((\text{p}_I = \pm 3)\) transitions, etc.
While specifying \(\text{p}_I\) alone is not enough to select an individual singlespin transition, any individual singlespin transition can be identified by a combination of \(\text{p}_I\) and the singlespin transition symmetry function \(\text{d}_I\), given by
You can verify this from the values of \(\text {p}_I\) and \(\text{d}_I\) for all singlespin transitions for \(I=1\), \(I=3/2\) and \(I=5/2\) shown below. Note that \(\text{d}_I = 0\) for all transitions in a \(I=1/2\) nucleus.
For a summary on spin transition symmetry functions in NMR, click on the disclosure button below.
Spin Transition Symmetry Functions
Note
In the symmetry pathway approach, the idea of coherence order is extended to form a complete set of spin transition symmetry functions, \(\xi_\ell (i,j)\), given by
where the \(\hat{T}_{\ell,0}\) are irreducible tensor operators. The function symbol \(\xi_\ell(i,j)\) is replaced with the lowercase symbols \(\mathbb{p}(i,j)\), \(\mathbb{d}(i,j)\), \(\mathbb{f} (i,j)\), \(\ldots\), i.e., we follow the spectroscopic subshell letter designations:
To simplify usage in figures and discussions, we scale the transition symmetry functions to integers values according to
The \(\ell=0\) function is dropped as it always evaluates to zero. For a single spin, \(I\), a complete set of functions is defined up to \(\ell = 2I\).
For weakly coupled nuclei, we define the transition symmetry functions
Replacing the symmetry function symbol using subshell letter designations becomes more cumbersome in this case. When the \(\ell\) are zero on all nuclei except one, we identify these functions as
For weakly coupled homonuclear spins it is also convenient to define
When the \(\ell\) are zero on all nuclei except two, then we identify these functions using a concatenation of subshell letter designations, e.g.,
Below is an energy level diagram of two coupled spin \(I=1/2\) nuclei with transition labeled according to their transition symmetry function values. Note that each transition has a unique set of transition symmetry function values.
SingleSpin Queries¶
Based on the review above, we now know for the spin \(I=1\), the transition \(\ketbra{1}{0}\) can be selected with \({(\text{p}_I,\text{d}_I) = (1,1)}\). In mrsimulator, this transition is selected during a SpectralEvent using the SymmetryQuery and TransitionQuery objects, as defined in the code below.
from mrsimulator.method.query import SymmetryQuery, TransitionQuery
from mrsimulator.method import SpectralEvent
symm_query = SymmetryQuery(P=[1], D=[1])
trans_query = TransitionQuery(ch1=symm_query)
spec_event = SpectralEvent(transition_queries=[trans_query])
Note
Python dictionaries can also be used to create and initialize mrsimulator objects. To do this, the dictionary must use the object’s attribute names as the key strings and be passed to a higher level object. Since a SpectralEvent object holds a list of TransitionQuery objects, the above code could have been written as
# Both SymmetryQuery and TransitionQuery object as dict
symm_query_dict = {"P": [1], "D": [1]}
trans_query_dict = {"ch1": symm_query_dict}
# Dictionary of TransitionQuery passed to SpectralEvent
spec_event = SpectralEvent(transition_queries=[trans_query_dict])
In the example above, the SymmetryQuery object is created and assigned to
the TransitionQuery attribute ch1
, i.e., it acts on the isotope in the
“first channel”. Recall that the channels
attribute of the Method object
holds an ordered list of isotope strings. This list’s first, second, and third
isotopes are associated with ch1
, ch2
, and ch3
, respectively.
Currently, mrsimulator only supports up to three channels, although this may
be increased in future versions.
The TransitionQuery object goes into a list in the
transition_queries
attribute of a SpectralEvent object. The
SpectralEvent object, in turn, is added to an ordered list in the events
attribute of a SpectralDimension object. All this is illustrated in the code
below.
from mrsimulator import Site, Coupling, SpinSystem, Simulator
from mrsimulator import Method, SpectralDimension
from mrsimulator import signal_processor as sp
import matplotlib.pyplot as plt
import numpy as np
# Create single Site and Spin System
deuterium = Site(
isotope="2H",
isotropic_chemical_shift=10, # in ppm
shielding_symmetric={"zeta": 80, "eta": 0.25}, # zeta in ppm
quadrupolar={"Cq": 10e3, "eta": 0.0, "alpha": 0, "beta": np.pi / 2, "gamma": 0},
)
deuterium_system = SpinSystem(sites=[deuterium])
# This method selects all observable (p_I = –1) transitions
method_both_transitions = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=40000, # in Hz
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1]}}])],
)
],
)
# This method selects observable (p_I = –1) transitions with d_I = 1
method_transition1 = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=40000, # in Hz
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [1]}}])],
)
],
)
# This method selects observable (p_I = –1) transitions with d_I = 1
method_transition2 = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=40000, # in Hz
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [1]}}])],
)
],
)
# Simulate spectra for all three method with spin system
sim = Simulator(
spin_systems=[deuterium_system],
methods=[method_both_transitions, method_transition1, method_transition2],
)
sim.run()
# Create SignalProcessor for Gaussian Convolution
processor = sp.SignalProcessor(
operations=[sp.IFFT(), sp.apodization.Gaussian(FWHM="100 Hz"), sp.FFT()]
)
# Plot spectra from all three methods
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), subplot_kw={"projection": "csdm"})
ax[0].plot(
processor.apply_operations(dataset=sim.methods[0].simulation).real,
label="$p_I = 1$ transition",
)
ax[0].set_title("SingleQuantum Spectrum All Transitions")
ax[0].legend()
ax[0].grid()
ax[0].invert_xaxis() # reverse xaxis
ax[1].plot(
processor.apply_operations(dataset=sim.methods[1].simulation).real,
label="$(p_I,d_I) = (1,+1)$ transitions",
)
ax[1].plot(
processor.apply_operations(dataset=sim.methods[2].simulation).real,
label="$(p_I,d_I) = (1,1)$ transitions",
)
ax[1].set_title("SingleQuantum Spectrum Single Transitions")
ax[1].legend()
ax[1].grid()
ax[1].invert_xaxis() # reverse xaxis
ax[0].set_ylim(0.02, 0.34) # Set ylimits to be the same
ax[1].set_ylim(0.02, 0.34) # on both plots
plt.tight_layout()
plt.show()
Note
Whenever the D
attribute is omitted, the SymmetryQuery allows
transitions with all values of \(\text{d}_I\). On the other hand, whenever
the P
attribute is omitted, it defaults to P=[0]
, i.e., no selected
transitions on the assigned channel.
Selecting Symmetric SingleSpin Transitions¶
A notable case, particularly useful for halfinteger quadrupolar nuclei, is that \(\text{d}_I = 0\) for all symmetric \((m \rightarrow  m)\) transitions, as these transitions are unaffected by the firstorder quadrupolar coupling frequency contribution. The MQMAS experiment involves a 2D correlation of the two symmetric (\(\text{d}_I = 0\)) transitions, \(\ketbra{\tfrac{1}{2}}{\tfrac{1}{2}}\), the socalled “central transition,” and \(\ketbra{\tfrac{3}{2}}{\tfrac{3}{2}}\), the symmetric triple quantum transition. The code below is an example of a custom 2D method using two SpectralDimension objects, each holding a single SpectralEvent. The TransitionQuery objects select each transition in their respective SpectralDimension objects.
my_mqmas = Method(
channels=["87Rb"],
magnetic_flux_density=9.4,
rotor_frequency=np.inf, # in Hz (here, set to infinity)
spectral_dimensions=[
SpectralDimension(
count=128,
spectral_width=6e3, # in Hz
reference_offset=9e3, # in Hz
label="Symmetric 3Q Frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [3], "D": [0]}}])],
),
SpectralDimension(
count=256,
spectral_width=6e3, # in Hz
reference_offset=5e3, # in Hz
label="Central Transition Frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [0]}}])],
),
],
)
# Create three sites in RbNO3
site1 = Site(
isotope="87Rb",
isotropic_chemical_shift=27.4, # ppm
quadrupolar={"Cq": 1.68e6, "eta": 0.2}, # Cq in Hz
)
site2 = Site(
isotope="87Rb",
isotropic_chemical_shift=28.5, # ppm
quadrupolar={"Cq": 1.94e6, "eta": 1}, # Cq in Hz
)
site3 = Site(
isotope="87Rb",
isotropic_chemical_shift=31.3, # ppm
quadrupolar={"Cq": 1.72e6, "eta": 0.5}, # Cq in Hz
)
# No Couplings, so create a separate SpinSystem for each site.
sites = [site1, site2, site3]
RbNO3_spin_systems = [SpinSystem(sites=[s]) for s in sites]
sim = Simulator(spin_systems=RbNO3_spin_systems, methods=[my_mqmas])
sim.run()
# Apply Gaussian line broadening along both dimensions via convolution
gauss_convolve = sp.SignalProcessor(
operations=[
sp.IFFT(dim_index=(0, 1)),
sp.apodization.Gaussian(FWHM="0.08 kHz", dim_index=0),
sp.apodization.Gaussian(FWHM="0.22 kHz", dim_index=1),
sp.FFT(dim_index=(0, 1)),
]
)
dataset = gauss_convolve.apply_operations(dataset=sim.methods[0].simulation)
In the code below, we use the imshow() from the
matplotlib.pyplot module to
return an image of the dataset on a 2D regular raster. We also use "gist_ncar_r"
from
matpltolib’s included colormaps
to map the dataset amplitude to colors; the colorbar()
function provides the visualization of the dataset mapping to color to the right
of the plot.
plt.figure(figsize=(4, 3))
ax = plt.subplot(projection="csdm")
cb = ax.imshow(dataset.real / dataset.real.max(), aspect="auto", cmap="gist_ncar_r")
plt.colorbar(cb)
ax.invert_xaxis()
ax.invert_yaxis()
plt.tight_layout()
plt.show()
Warning
This custom method, as well as the builtin MultiQuantum VAS methods, assumes uniform excitation and mixing of the multiplequantum transition. In an experimental MQMAS measurement, excitation and mixing efficiencies depend on the ratio of the quadrupolar coupling constant to the rf field strength. Therefore, the relative integrated intensities of this simulation may not agree with the experiment.
Inspecting Transition and Symmetry Pathways¶
You can view the symmetry pathways that will be selected by your custom method
in a given spin system using the function
get_symmetry_pathways()
as shown below.
from pprint import pprint
pprint(my_mqmas.get_symmetry_pathways("P"))
pprint(my_mqmas.get_symmetry_pathways("D"))
Out:
[SymmetryPathway(
ch1(87Rb): [3] ⟶ [1]
total: 3.0 ⟶ 1.0
)]
[SymmetryPathway(
ch1(87Rb): [0] ⟶ [0]
total: 0.0 ⟶ 0.0
)]
Similarly, you can view the transition pathway that will be selected by your custom method in a given
spin system using the function get_transition_pathways()
as shown
below.
from pprint import pprint
pprint(my_mqmas.get_transition_pathways(SpinSystem(sites=[site1])))
Out:
[1.5⟩⟨1.5 ⟶ 0.5⟩⟨0.5, weight=(1+0j)]
MultiSpin Queries¶
When there is more than one site in a spin system, things get a little more complicated with the SymmetryQuery objects. Here we review some important concepts associated with transition symmetry functions in coupled spin systems, and see how SymmetryQuery objects are designed to work in such cases.
SingleSpin SingleQuantum Transitions¶
Consider the case of three weakly coupled proton sites. Here, the selection rule for observable transitions is
These corresponds to the singlespin singlequantum transitions labeled \(\hat{A}_1\), \(\hat{A}_2\), \(\hat{A}_3\), \(\hat{A}_4\), \(\hat{M}_1\), \(\hat{M}_2\), \(\hat{M}_3\), \(\hat{M}_4\), \(\hat{X}_1\), \(\hat{X}_2\), \(\hat{X}_3\), and \(\hat{X}_4\) in the energy level diagram below.
Keep in mind that the Method object does not know, in advance, the number of sites in a spin system.
The TransitionQuery for selecting these 12 singlespin singlequantum transitions is given in the code below.
# Create Site, Coupling and SpinSystem objects
site_A = Site(isotope="1H", isotropic_chemical_shift=0.5)
site_M = Site(isotope="1H", isotropic_chemical_shift=2.5)
site_X = Site(isotope="1H", isotropic_chemical_shift=4.5)
sites = [site_A, site_M, site_X]
coupling_AM = Coupling(site_index=[0, 1], isotropic_j=12)
coupling_AX = Coupling(site_index=[0, 2], isotropic_j=12)
coupling_MX = Coupling(site_index=[1, 2], isotropic_j=12)
couplings = [coupling_AM, coupling_AX, coupling_MX]
proton_system = SpinSystem(sites=sites, couplings=couplings)
# Custom method object emulating a onepulse acquire experiment
method = Method(
channels=["1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=1800, # in Hz
reference_offset=1000, # in Hz
label="$^{1}$H frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1]}}])],
)
],
)
sim = Simulator(spin_systems=[proton_system], methods=[method])
sim.run()
# Add line broadening
processor = sp.SignalProcessor(
operations=[sp.IFFT(), sp.apodization.Exponential(FWHM="1 Hz"), sp.FFT()]
)
plt.figure(figsize=(10, 3)) # set the figure size
ax = plt.subplot(projection="csdm")
ax.plot(processor.apply_operations(dataset=sim.methods[0].simulation).real)
ax.invert_xaxis() # reverse xaxis
plt.tight_layout()
plt.grid()
plt.show()
The assignment of transitions in the spectrum above are, from left to right, are \(\hat{X}_4, (\hat{X}_3, \hat{X}_2)\), and \(\hat{X}_1\) centered at 4.5 ppm, \(\hat{M}_4, (\hat{M}_3, \hat{M}_2)\), and \(\hat{M}_1\) centered at 2.5 ppm, and \(\hat{A}_4, (\hat{A}_3, \hat{A}_2)\), and \(\hat{A}_1\) centered at 0.5 ppm.
It is essential to realize that all sites having the same isotope are
“indistinguishable” to a TransitionQuery object. Recall that ch1
is
associated with the first isotope in the list of isotope strings assigned to the
Method attribute channels
. When the TransitionQuery above is combined with
the SpinSystem object with three \(^1\text{H}\) Sites, it must first expand
its SymmetryQuery into an intermediate set of spinsystemspecific symmetry
queries, illustrated by each row in the table below.
Transitions 
\(\text{p}_A\) 
\(\text{p}_M\) 
\(\text{p}_X\) 

\(\hat{A}_1, \hat{A}_2, \hat{A}_3, \hat{A}_4\) 
–1 
0 
0 
\(\hat{M}_1, \hat{M}_2, \hat{M}_3, \hat{M}_4\) 
0 
–1 
0 
\(\hat{X}_1, \hat{X}_2, \hat{X}_3, \hat{X}_4\) 
0 
0 
–1 
The intermediate spinsystemspecific symmetry query in each row selects a subset of transitions from the complete set of transitions. The final set of selected transitions is obtained from the union of transition subsets from each spinsystemspecific symmetry query.
The get_transition_pathways()
function will allow
you to inspect the transitions selected by the Method
in terms of the initial and final Zeeman eigenstate quantum numbers.
from pprint import pprint
pprint(method.get_transition_pathways(proton_system))
Out:
[0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j)]
To further illustrate how the TransitionQuery and SymmetryQuery objects work in a multisite spin system, let’s examine a few more examples in the case of three weakly coupled proton sites.
TwoSpin DoubleQuantum Transitions¶
In this spin system there are six twospin doublequantum transitions where \(\text{p}_{AMX} = \text{p}_{A} + \text{p}_{M} + \text{p}_{X} = 2\) and another six twospin doublequantum transitions where \(\text{p}_{AMX} = \text{p}_{A} + \text{p}_{M} + \text{p}_{X} = +2\). The \(\text{p}_{AMX} = 2\) transitions are illustrated in the energylevel diagram below.
The code below will select the six twospin doublequantum transitions where \(\text{p}_{AMX} = 2\).
method = Method(
channels=["1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=2000, # in Hz
reference_offset=2000, # in Hz
label="$^{1}$H frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1, 1]}}])],
)
],
)
sim = Simulator(spin_systems=[proton_system], methods=[method])
sim.run()
plt.figure(figsize=(10, 3)) # set the figure size
ax = plt.subplot(projection="csdm")
ax.plot(processor.apply_operations(dataset=sim.methods[0].simulation).real)
ax.invert_xaxis() # reverse xaxis
plt.tight_layout()
plt.grid()
plt.show()
The assignment of transitions in the spectrum above are, from left to right, are \(\hat{D}_{2,MX}\), \(\hat{D}_{1,MX}\), \(\hat{D}_{2,AX}\), \(\hat{D}_{1,AX}\), \(\hat{D}_{2,AM}\), and \(\hat{D}_{1,AM}\),
As before, when this generic TransitionQuery is combined with the threesite SpinSystem object, the SymmetryQuery is expanded into an intermediate set of spinsystemspecific symmetry queries illustrated in the table below.
Transitions 
\(\text{p}_A\) 
\(\text{p}_M\) 
\(\text{p}_X\) 

\(\hat{D}_{1,AM}, \hat{D}_{2,AM}\) 
–1 
–1 
0 
\(\hat{D}_{1,MX}, \hat{D}_{2,MX}\) 
0 
–1 
–1 
\(\hat{D}_{1,AX}, \hat{D}_{2,AX}\) 
–1 
0 
–1 
Again, each row’s intermediate spinsystemspecific symmetry query is used to select a subset of transitions from the complete set of transitions. The final set of selected transitions is obtained from the union of transition subsets from each spinsystemspecific symmetry query.
from pprint import pprint
pprint(method.get_transition_pathways(proton_system))
Out:
[0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j)]
ThreeSpin SingleQuantum Transitions¶
Another interesting example in this spin system with three weakly coupled proton sites are the three threespin singlequantum transitions having \(\text{p}_{AMX} = \text{p}_{A} + \text{p}_{M} + \text{p}_{X} = 1\) and the three threespin singlequantum transitions having \(\text{p}_{AMX} = \text{p}_{A} + \text{p}_{M} + \text{p}_{X} = +1\)
The three threespin singlequantum transitions having \(\text{p}_{AMX} = 1\) are illustrated in the energy level diagram below.
The code below will select these threespin singlequantum transitions.
method = Method(
channels=["1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=4000, # in Hz
reference_offset=1000, # in Hz
label="$^{1}$H frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1, 1, +1]}}])]
)
],
)
sim = Simulator(spin_systems=[proton_system], methods=[method])
sim.run()
plt.figure(figsize=(10, 3)) # set the figure size
ax = plt.subplot(projection="csdm")
ax.plot(processor.apply_operations(dataset=sim.methods[0].simulation).real)
ax.invert_xaxis() # reverse xaxis
plt.tight_layout()
plt.grid()
plt.show()
The assignment of transitions in the spectrum above are, from left to right, are \(\hat{S}_{3,AMX}\), \(\hat{S}_{2,AMX}\), and \(\hat{S}_{1,AMX}\)
Again, combined with the threesite SpinSystem object, the SymmetryQuery is expanded into the set of spinsystemspecific symmetry queries illustrated in the table below.
Transitions 
\(\text{p}_A\) 
\(\text{p}_M\) 
\(\text{p}_X\) 

\(\hat{S}_{1,AMX}\) 
–1 
+1 
–1 
\(\hat{S}_{2,AMX}\) 
–1 
–1 
+1 
\(\hat{S}_{3,AMX}\) 
+1 
–1 
–1 
from pprint import pprint
pprint(method.get_transition_pathways(proton_system))
Out:
[0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j),
0.5, 0.5, 0.5⟩⟨0.5, 0.5, 0.5, weight=(1+0j)]
As you can surmise from the examples, the attributes of SymmetryQuery, P
and
D
, hold a list of singlespin transition symmetry function values, and the
length of the list is the desired number of spins that are involved in the
transition.
Heteronuclear multiplespin transitions¶
How does D
fit into the multisite SymmetryQuery story? Consider the
case of two coupled hydrogens, except we replace one of the \(^1H\) with
\(^2H\). Let’s focus on the singlespin singlequantum transitions, shown
below as \(\hat{A}_{1\pm}\) and \(\hat{A}_{2\pm}\) on the left, and the
twospin triplequantum transition, shown below as \(\hat{T}_{AX}\) on
the right.
from mrsimulator.spin_system.tensors import SymmetricTensor
import numpy as np
site_A = Site(
isotope="2H",
isotropic_chemical_shift=0.5,
quadrupolar=SymmetricTensor(
Cq=100000, # in Hz
eta=0.2,
alpha=5 * np.pi / 180,
beta=np.pi / 2,
gamma=70 * np.pi / 180,
),
)
site_X = Site(isotope="1H", isotropic_chemical_shift=4.5)
sites = [site_A, site_X]
coupling_AX = Coupling(site_index=[0, 1], dipolar={"D": 20000})
couplings = [coupling_AX]
system_AX = SpinSystem(sites=sites, couplings=couplings)
methodAll1Q = Method(
channels=["2H", "1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=200000, # in Hz
reference_offset=0, # in Hz
label="$^{2}$H frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1]}}])],
)
],
)
methodHalf1Q = Method(
channels=["2H", "1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=200000, # in Hz
reference_offset=0, # in Hz
label="$^{2}$H frequency",
events=[SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [1]}}])],
)
],
)
method3Q = Method(
channels=["2H", "1H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=16000,
spectral_width=10000, # in Hz
reference_offset=5000, # in Hz
label="$^{2}$H frequency",
events=[
SpectralEvent(
transition_queries=[{"ch1": {"P": [2]}, "ch2": {"P": [1]}}]
)
],
)
],
)
processor = sp.SignalProcessor(
operations=[sp.IFFT(), sp.apodization.Gaussian(FWHM="100 Hz"), sp.FFT()]
)
sim = Simulator(spin_systems=[system_AX], methods=[methodAll1Q, methodHalf1Q, method3Q])
sim.config.integration_volume = "hemisphere"
sim.run()
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), subplot_kw={"projection": "csdm"})
ax[0].plot(processor.apply_operations(dataset=sim.methods[0].simulation).real)
ax[0].set_title("Full SingleQuantum Spectrum")
ax[0].grid()
ax[0].invert_xaxis() # reverse xaxis
ax[1].plot(processor.apply_operations(dataset=sim.methods[1].simulation).real)
ax[1].set_title("Half SingleQuantum Spectrum")
ax[1].grid()
ax[1].invert_xaxis() # reverse xaxis
plt.tight_layout()
plt.show()
The deuterium spectrum of a staticpolycrystalline sample is shown on the left for all singlespin singlequantum transitions on deuterium, \(\hat{A}_{1\pm}\) and \(\hat{A}_{2\pm}\). The spectrum on the right is for half of the singlespin singlequantum transitions on deuterium: \(\hat{A}_{1}\) and \(\hat{A}_{2}\).
plt.figure(figsize=(10, 3)) # set the figure size
ax = plt.subplot(projection="csdm")
ax.set_title("Heteronuclear TwoSpin ($^2$H$^1$H) TripleQuantum Spectrum")
ax.plot(processor.apply_operations(dataset=sim.methods[2].simulation).real)
plt.tight_layout()
plt.grid()
plt.show()
Transitions 
\(\text{p}_A\) 
\(\text{d}_A\) 
\(\text{p}_X\) 

\(\hat{T}_{AX}\) 
–2 
0 
–1 
The single transition in the heteronuclear twospin (\(^2\text{H}\)\(^1\text{H}\)) triplequantum spectrum is unaffected by the dipolar and quadrupolar frequency anisotropies.
Frequency Contributions¶
The NMR frequency, \(\Omega(\Theta,i,j)\), of an \(i \rightarrow j\) transition between the eigenstates of the stationarystate semiclassical Hamiltonian in a sample with a lattice spatial orientation, \(\Theta\), can be written as a sum of components,
with each component, \(\Omega_k(\Theta,i,j)\), separated into three parts:
where \({\xi}^{(k)}_\ell(i,j)\) are the spin transition symmetry functions described earlier, \({\Xi}^{(k)}_L(\Theta)\) are the spatial symmetry functions, and \(\omega_k\) gives the size of the kth frequency component. The experimentalist indirectly influences a frequency component \(\Omega_k\) by direct manipulation of the quantum transition, \(i \rightarrow j\), and the spatial orientation, \(\Theta\), of the sample.
The function symbol \(\Xi_\ell(\Theta)\) is replaced with the uppercase symbols \(\mathbb{S}\), \(\mathbb{P}(\Theta)\), \(\mathbb{D}(\Theta)\), \(\mathbb{F}(\Theta)\), \(\mathbb{G}(\Theta)\), \(\ldots\), i.e., following the spectroscopic subshell letter designations for \(L\). Consult the Symmetry Pathways paper for more details on the form of the spatial symmetry functions. In short, the \(\mathbb{S}\) function is independent of sample orientation, i.e., it will appear in all isotropic frequency contributions. The \(\mathbb{D}(\Theta)\) function has a secondrank dependence on sample orientation, and can be averaged away with fast magicangle spinning, i.e., spinning about an angle, \(\theta_R\), that is the root of the secondrank Legendre polynomial \(P_2(\cos \theta_R)\). The other spatial symmetry functions are removed by spinning the sample about the corresponding root of the \(L^\text{th}\)rank Legendre polynomial \(P_L(\cos \theta_R)\)
Note
For 2ndorder quadrupolar coupling contributions, it is convenient to define “hybrid” spin transition functions as linear combinations of the spin transition functions
These transition symmetry functions play an essential role in evaluating the
individual frequency contributions to the overall transition frequency, given in
the table below and in
FrequencyEnum()
. They also aid in
pulse sequence design by identifying how different frequency contributions
refocus through the transition pathways.
For a summary on echo symmetry classification in NMR, click on the disclosure button below.
Echo Symmetry Classification
Note
The wellknown Hahnecho can occur whenever the \(p_I\) values of transitions in a transition pathway change sign. This is because the changing sign of \(p_I\) leads to a sign change for every \(p_I\)dependent transition frequency contribution. Thus, a Hahn echo forms whenever
assuming a frequency contribution’s spatial symmetry function, \({\Xi}\), remains constant during this period. As seen in the table in the Frequency Contributions table, sign changes in other symmetry functions can also lead to corresponding sign changes for dependent frequency contributions. Thus, a problem with showing only the \(p_I\) symmetry pathway for an NMR method is that it does not explain the formation of other classes of echoes that result when other symmetry functions change sign in a transition pathway. To fully understand when and which frequency contributions refocus into echoes, we must follow all relevant spatial, transition, or spatialtransition product symmetries through an NMR experiment. Thus, we generally classify echoes that refocus during a time interval as a transition symmetry echo (at constant \({\Xi}_k\)) when
and as a spatial symmetry echo (at constant \({\xi}_k\)) when
and as a spatialtransition symmetry product echo when
Within the class of transition echoes we find subclasses such as \(\text{p}\) echoes, which include the Hahn echo and the stimulated echo; \(\text{d}\) echoes, which include the solid echo and Solomon echoes, \(\text{c}_4\) echoes, used in MQMAS and SatelliteTransition MagicAngle Spinning (STMAS); \(\text{c}_2\) echoes, used in Correlation Of Anisotropies Separated Through Echo Refocusing (COASTER); and \(\text{c}_0\) echoes, used in MultipleQuantum DOuble Rotation (MQDOR).
Within the class of spatial echoes we find subclasses such as \(\mathbb{D}\) rotary echoes, which occur during sample rotation, and \(\mathbb{D}_0\) and \(\mathbb{G}_0\) echoes, which are designed to occur simultaneously during the DynamicAngle Spinning (DAS) experiment.
Interactions 
perturbation 
anisotropy 

Expression 

order 
rank 

shielding 
1st 
0th 

\(\omega_0 \sigma_\text{iso} \cdot \mathbb{p}_I\) 
shielding 
1st 
2nd 

\(\omega_0 \zeta_\sigma \cdot \mathbb{D}^{\{\sigma\}} \cdot \mathbb{p}_I\) 
weak J 
1st 
0th 

\(2 \pi J_\text{iso} \, (\mathbb{pp})_{IS}\) 
weak J 
1st 
2nd 

\(2 \pi \zeta_J \cdot \mathbb{D}^{\{d_{IS}\}} \cdot (\mathbb{pp})_{IS}\) 
weak dipolar 
1st 
2nd 

\(\omega_d \cdot \mathbb{D}^{\{d_{IS}\}} \cdot (\mathbb{pp})_{IS}\) 
quadrupolar 
1st 
2nd 

\(\omega_q \cdot \mathbb{D}^{\{q\}} \cdot \mathbb{d}_I\) 
quadrupolar 
2nd 
0th 

\(\displaystyle \frac{\omega_q^2}{\omega_0} \cdot \mathbb{S}^{\{qq\}} \cdot \mathbb{c}_0\) 
quadrupolar 
2nd 
2nd 

\(\displaystyle\frac{\omega_q^2}{\omega_0} \cdot \mathbb{D}^{\{qq\}} \cdot \mathbb{c}_2\) 
quadrupolar 
2nd 
4th 

\(\displaystyle\frac{\omega_q^2}{\omega_0} \cdot \mathbb{G}^{\{qq\}} \cdot \mathbb{c}_4\) 
Affine Transformations¶
The ability to refocus different spatial and transition symmetries into echoes with different paths in timeresolved NMR experiments creates opportunities for generating multidimensional spectra that correlate different interactions. These spectra can be made easier to interpret through similarity transformations. Most similarity transformations in NMR are affine transformations, as they preserve the colinearity of points and ratios of distances. Essential in any similarity transformation is whether to implement the transformation actively or passively. Active transformations change the appearance of the signal while leaving the coordinate system unchanged, whereas passive transformations leave the appearance of the signal unchanged while changing the coordinate system. Both active and passive transformations are used extensively in NMR.
The general form of the affine transformation of a ndimensional spectrum is
In the twodimensional case, this is given by
Note
For the multiplequantum MAS experiment, a shear and scale transformation is often applied to the spectrum to create a 2D spectrum correlating the MQMAS isotropic frequency to the anisotropic central transition frequency. This correlation can be achieved by adding an affine matrix to the method.
For 3QMAS on a spin \(I=3/2\) nucleus, where the shear factor is \(\kappa^{(\omega_2)} = 21/27\), the affine matrix giving the appropriate shear and scale transformation is given by
After the affine transformation, the position of the resonance in the isotropic projection is a weighted average of the multiple quantum and central transition isotropic frequencies given by
If the spectrum is to be referenced to a frequency other than the rf carrier frequency (i.e. zero is not defined in the middle of the spectrum), then the reference offset used in the singlequantum dimension must be multiplied by a factor of \({\left({\text{p}_I^{[1]}}/{\text{p}_I^{[2]}} + \kappa^{(\omega_1)} \right)/(1+ \kappa^{(\omega_1)} )}\) when used in the isotropic dimension.
See the “Symmetry Pathways in SolidState NMR” paper for a more detailed discussion on affine transformations in NMR.
In the code below, the 3QMAS method described earlier is modified to include an affine matrix to perform this shear transformation.
my_sheared_mqmas = Method(
channels=["87Rb"],
magnetic_flux_density=9.4,
rotor_frequency=np.inf, # in Hz (here, set to infinity)
spectral_dimensions=[
SpectralDimension(
count=128,
spectral_width=6e3, # in Hz
reference_offset=9e3, # in Hz
label="3QMAS isotropic dimension",
events=[
SpectralEvent(transition_queries=[{"ch1": {"P": [3], "D": [0]}}])
],
),
SpectralDimension(
count=256,
spectral_width=6e3, # in Hz
reference_offset=5e3, # in Hz
label="Central Transition Frequency",
events=[
SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [0]}}])
],
),
],
affine_matrix=[[9 / 16, 7 / 16], [0, 1]],
)
sim = Simulator(spin_systems=RbNO3_spin_systems, methods=[my_sheared_mqmas])
sim.run()
dataset = gauss_convolve.apply_operations(dataset=sim.methods[0].simulation)
plt.figure(figsize=(4, 3))
ax = plt.subplot(projection="csdm")
cb = ax.imshow(dataset.real / dataset.real.max(), aspect="auto", cmap="gist_ncar_r")
plt.colorbar(cb)
ax.invert_xaxis()
ax.invert_yaxis()
plt.tight_layout()
plt.show()
Note
For MQMAS, a second shear and scale can be applied to remove isotropic chemical shift component along the \(\Omega^{[2]''}\) axis. For a spin \(I=3/2\) nucleus, with a second shear factor of \(\kappa^{(\omega_1)} =  8/17\), the affine matrix is given by
and the product of the two affine transformations is
Below is the code for simulating a 3QMAS spectrum with a double shear transformation.
my_twice_sheared_mqmas = Method(
channels=["87Rb"],
magnetic_flux_density=9.4,
rotor_frequency=np.inf, # in Hz (here, set to infinity)
spectral_dimensions=[
SpectralDimension(
count=128,
spectral_width=6e3, # in Hz
reference_offset=9e3, # in Hz
label="3QMAS isotropic dimension",
events=[
SpectralEvent(transition_queries=[{"ch1": {"P": [3], "D": [0]}}])
],
),
SpectralDimension(
count=256,
spectral_width=6e3, # in Hz
reference_offset=0, # in Hz
label="CT QuadOnly Frequency",
events=[
SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [0]}}])
],
),
],
affine_matrix=[[9 / 16, 7 / 16], [9 / 50, 27 / 50]],
)
sim = Simulator(spin_systems=RbNO3_spin_systems, methods=[my_twice_sheared_mqmas])
sim.run()
dataset = gauss_convolve.apply_operations(dataset=sim.methods[0].simulation)
plt.figure(figsize=(4, 3))
ax = plt.subplot(projection="csdm")
cb = ax.imshow(dataset.real / dataset.real.max(), aspect="auto", cmap="gist_ncar_r")
plt.colorbar(cb)
ax.invert_xaxis()
ax.invert_yaxis()
plt.tight_layout()
plt.show()
Average Frequency & Multiple Events¶
To illustrate the versatility of the Method object, we can also design an MQMAS method that correlates the isotropic MQMAS frequency to the central transition without the need for an affine transformation. Recall that the 3QMAS isotropic frequency on spin \(I=3/2\) is given by
As we saw at the beginning of this section, the first spectral dimension derives its average frequency, \(\overline{\Omega}_1\), from a weighted average of multiple transition frequencies. Thus, this weighted average frequency can be obtained through the use of multiple SpectralEvent objects in the SpectralDimension associated with the isotropic dimension, as shown in the code below.
my_three_event_mqmas = Method(
channels=["87Rb"],
magnetic_flux_density=9.4,
rotor_frequency=np.inf, # in Hz (here, set to infinity)
spectral_dimensions=[
SpectralDimension(
count=128,
spectral_width=6e3, # in Hz
reference_offset=9e3, # in Hz
label="3QMAS isotropic dimension",
events=[
SpectralEvent(
fraction=9 / 16, transition_queries=[{"ch1": {"P": [3], "D": [0]}}]
),
SpectralEvent(
fraction=7 / 16, transition_queries=[{"ch1": {"P": [1], "D": [0]}}]
),
],
),
SpectralDimension(
count=256,
spectral_width=6e3, # in Hz
reference_offset=5e3, # in Hz
label="Central Transition Frequency",
events=[
SpectralEvent(transition_queries=[{"ch1": {"P": [1], "D": [0]}}])
],
),
],
)
sim = Simulator(spin_systems=RbNO3_spin_systems, methods=[my_three_event_mqmas])
sim.run()
dataset = gauss_convolve.apply_operations(dataset=sim.methods[0].simulation)
plt.figure(figsize=(4, 3))
ax = plt.subplot(projection="csdm")
cb = ax.imshow(dataset.real / dataset.real.max(), aspect="auto", cmap="gist_ncar_r")
plt.colorbar(cb)
ax.invert_xaxis()
ax.invert_yaxis()
plt.tight_layout()
plt.show()
We could apply an affine transformation to remove the isotropic chemical shift
from the central transition (horizontal) dimension. If you go back to the
previous discussion, you will find that the required value for the
affine_matrix
in the Method object to do this shear is given by
affine_matrix=[[1,0],[8/25, 17/25]]
Mixing Queries¶
The amplitude of a transition pathway signal derives from the product of mixing amplitudes associated with each transfer between transitions in a transition pathway, e.g.,
Here, \(a_{0A}\) is the amplitude of the initial \(\hat{t}_A\) transition, \(a_{AB}\) is the mixing amplitude for the transfer from \(\hat{t}_A\) to \(\hat{t}_B\), \(a_{BC}\) is the mixing amplitude for the transfer from \(\hat{t}_B\) to \(\hat{t}_C\), and so on. The growing product \((a_{0A}a_{AB}a_{BC} \cdots)\) is the transition pathway amplitude. Eliminating a transition with a TransitionQuery in a spectral or delay event sets the eliminated transition’s pathway amplitude to zero, i.e., it prunes that transition pathway branch.
Default Total Mixing between Adjacent Spectral or Delay Events¶
In previous discussions, we did not mention the efficiency of transfer between selected transitions in adjacent SpectralEvent objects. This is because, as default behavior, mrsimulator does a total mixing, i.e., connects all selected transitions in the two adjacent spectral or delay events. In other words, if the first of two adjacent SpectralEvent objects has three selected transitions, and the second has two selected transitions, then mrsimulator will make \(3 \times 2 = 6\) connections, i.e., six transition pathways passing from the first to second SpectralEvent objects.
Additionally, this total mixing assumes that every connection has a mixing amplitude of 1. This is unrealistic, but if used correctly gives a significant speedup in the simulation by avoiding the need to calculate mixing amplitudes.
Warning
The use of total mixing, i.e., the default mixing, can complicate the comparison of integrated intensities between different methods, depending on the selected transition pathways.
If this default mixing behavior had been explicitly shown in the previous example, the events list in the first SpectralDimension would have looked like the code below.
from mrsimulator.method import MixingEvent
events = [
SpectralEvent(fraction=9 / 16, transition_queries=[{"ch1": {"P": [3], "D": [0]}}]),
MixingEvent(query="TotalMixing"),
SpectralEvent(fraction=7 / 16, transition_queries=[{"ch1": {"P": [1], "D": [0]}}]),
MixingEvent(query="TotalMixing"),
]
Since only one transition was selected in each SpectralEvent, the expected (and default) behavior is that there is a mixing (transfer) of coherence between the symmetric triplequantum and central transitions, forming the desired transition pathway.
However, when multiple transition pathways are present in a method, you may need
more accurate mixing amplitudes when connecting selected transitions of adjacent
events. You may also need to prevent the undesired mixing of specific
transitions between two adjacent events. As described below, you can avoid a
"TotalMixing"
event by inserting MixingEvent object with a certain rotation
query.
Rotation Query¶
A rotation of \(\theta\) about an axis defined by \(\phi\) in the \(x\)\(y\) plane on a selected transition, \(\ketbra{I, m_f}{I, m_i}\), in a spectral or delay event transfers it to all selected transitions, \(\ketbra{I,m_f'}{I,m_i'}\) in the next spectral or delay event, according to
where \(\Delta p_I = p_I'  p_I\). From this expression, we obtain the complex mixing amplitude from \(\ketbra{I, m_f}{I, m_i}\) to \(\ketbra{I, m_f'}{I, m_i'}\) due to a rotation to be
From this expression, we note a few interesting and useful cases. One is the coherence transfer under a \(\pi\) rotation, given by
That is, a \(\pi\) rotation can make only one connection between transitions in adjacent events. It is also a special connection because the \(\text{p}_I\) transition symmetry value for the two transitions are equal but opposite in sign. Additionally, the \(\text{d}_I\) transition symmetry remains unchanged (\(\Delta \text{d}_I = 0\)) for the two transitions. (In fact, this behavior under a \(\pi\) rotation is generally true for odd (\(\text{p}_I, \text{f}_I, \ldots)\) and even (\(\text{d}_I, \text{g}_I, \ldots)\) rank spin transition symmetry functions.)
Another interesting result is that, while a rotation can transfer a transition into many other transitions, the \(\text{d}_I\) transition symmetry value cannot remain unchanged (\(\Delta \text{d}_I \neq 0\)) between two connected transitions under a \(\pi/2\) rotation.
Finally, another useful result is
While it’s not surprising that a rotation through an angle of zero does nothing
to the transition, this turns out to help act as the opposite of a total mixing
event, i.e., a "NoMixing"
event. As a convenience, this is defined as a
"NoMixing"
query and can be implemented with the code below.
MixingEvent(query="NoMixing")
The MixingEvent object holds the rotation details in a MixingQuery object as
a RotationQuery object associated with a channels
attribute. This is
illustrated in the sample code below.
import numpy as np
from mrsimulator.method.query import RotationQuery
rot_query_90 = RotationQuery(angle=np.pi/2, phase=0)
rot_query_180 = RotationQuery(angle=np.pi, phase=0)
rot_mixing = MixingEvent(query={
"ch1": rot_query_90,
"ch2": rot_query_180
}
)
p and d Echoes on Deuterium¶
Here, we examine two examples in a deuterium spin system that illustrate the importance of echo classification in understanding how transitionfrequency contributions can be eliminated or separated based on their dependence on different transition symmetry functions.
First, we implement two Method objects that follow the design of two experimental pulse sequences. In this effort, we use RotationQuery objects to select the desired transition pathways and obtain spectra with the desired average frequencies. Then, we implement two simpler Method objects that produce identical spectra and illustrate how frequency contributions can be used to reduce the number of events needed in a custom method.
Consider the Hahn and Solid Echo pulse sequences on the left and right, respectively.
The Hahn Echo sequence, with \(\pi/2\tau\pit\rightarrow\), leads to the formation of a \(\text{p}_I\) echo at \(t = \tau\). The two transition pathways created by this experiment on a deuterium nucleus are illustrated beneath the sequence. Remember that a \(\pi\) rotation is a special because it connects transitions with equal but opposite signs of \(\text{p}_I\) while \(\text{d}_I\) remains invariant.
The Solid Echo sequence, with \(\pi/2\tau\pi/2t\rightarrow\), leads to the formation of a \(\text{d}_I\) echo at \(t = \tau\). The two transition pathways created by this experiment on a deuterium nucleus are illustrated beneath the sequence. Here, also recall that the \(\text{d}_I\) transition symmetry value cannot remain unchanged (\(\Delta \text{d}_I \neq 0\)) between two connected transitions under a \(\pi/2\) rotation.
Below are two custom Method objects for simulating the Hahn and Solid Echo experiments. There is only one SpectralDimension object in each method, and the average frequency during each spectral dimension is derived from equal fractions of two SpectralEvent objects. Between these two SpectralEvent objects is a MixingEvent with a RotationQuery object. The RotationQuery object is created with a \(\pi\) rotation in the Hahn Echo method, and a \(\pi/2\) rotation in the Solid Echo method.
Note
The transition_queries
attribute of SpectralEvent holds a list of
TransitionQuery objects. Each TransitionQuery in the list applies to
the full set of transitions in the spin system. The union of these transition
subsets becomes the final set of selected transitions during the
SpectralEvent.
We use the deuterium Site defined earlier in this document.
from mrsimulator.method import MixingEvent
deuterium = Site(
isotope="2H",
isotropic_chemical_shift=10, # in ppm
shielding_symmetric={"zeta": 80, "eta": 0.25}, # zeta in ppm
quadrupolar={"Cq": 10e3, "eta": 0.0, "alpha": 0, "beta": np.pi / 2, "gamma": 0},
)
deuterium_system = SpinSystem(sites=[deuterium])
hahn_echo = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=2e4, # in Hz
events=[
SpectralEvent(
fraction=0.5,
transition_queries=[
{"ch1": {"P": [1], "D": [1]}},
{"ch1": {"P": [1], "D": [1]}},
],
),
MixingEvent(query={"ch1": {"angle": 3.141592, "phase": 0}}),
SpectralEvent(
fraction=0.5,
transition_queries=[
{"ch1": {"P": [1], "D": [1]}},
{"ch1": {"P": [1], "D": [1]}},
],
),
],
)
],
)
solid_echo = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=2e4, # in Hz
events=[
SpectralEvent(
fraction=0.5,
transition_queries=[
{"ch1": {"P": [1], "D": [1]}},
{"ch1": {"P": [1], "D": [1]}},
],
),
MixingEvent(query={"ch1": {"angle": 3.141592 / 2, "phase": 0}}),
SpectralEvent(
fraction=0.5,
transition_queries=[
{"ch1": {"P": [1], "D": [1]}},
{"ch1": {"P": [1], "D": [1]}},
],
),
],
)
],
)
We can check the resulting transition pathways using these TransitionQuery objects with the
code below for the hahn_echo
method,
from pprint import pprint
pprint(hahn_echo.get_transition_pathways(deuterium_system))
Out:
[1.0⟩⟨0.0 ⟶ 1.0⟩⟨0.0, weight=(1+0j)
0.0⟩⟨1.0 ⟶ 0.0⟩⟨1.0, weight=(1+0j)]
and for the solid_echo
method with the code below.
pprint(solid_echo.get_transition_pathways(deuterium_system))
Out:
[1.0⟩⟨0.0 ⟶ 0.0⟩⟨1.0, weight=(0.5+0j)
0.0⟩⟨1.0 ⟶ 1.0⟩⟨0.0, weight=(0.5+0j)]
Notice that the weights of the transition pathways in the solidecho method are half of those in the Hahnecho method. This is because the \(\pi\) pulse in the Hahnecho method gives perfect transfer between the two transitions in the adjacent spectral events. In contrast, while the \(\pi/2\) pulse in the solidecho method prevents the undesired transition pathways with \(\Delta \text{d}_I = 0\), it also connects the selected transitions during the first spectral event to undesired transitions in the second spectral event, which are eliminated by its symmetry query.
Next, we simulate both methods, and perform a Gaussian line shape convolution on each output spectrum, and plot the datasets.
sim = Simulator()
sim.spin_systems = [deuterium_system]
sim.methods = [hahn_echo, solid_echo]
sim.run()
processor = sp.SignalProcessor(
operations=[
sp.IFFT(),
sp.apodization.Gaussian(FWHM="100 Hz"),
sp.FFT(),
]
)
hahn_dataset = processor.apply_operations(dataset=sim.methods[0].simulation)
solid_dataset = processor.apply_operations(dataset=sim.methods[1].simulation)
fig, ax = plt.subplots(1, 2, subplot_kw={"projection": "csdm"}, figsize=[8.5, 3])
ax[0].set_title("HahnEcho Spectrum")
ax[0].plot(hahn_dataset.real)
ax[0].invert_xaxis()
ax[0].grid()
ax[1].set_title("SolidEcho Spectrum")
ax[1].plot(solid_dataset.real)
ax[1].invert_xaxis()
ax[1].grid()
plt.tight_layout()
plt.show()
In the Hahnecho spectrum, the \(\text{p}_I\)dependent frequency contributions (i.e., the shielding contributions) were averaged to zero, leaving only the \(\text{d}_I\)dependent frequency contributions (i.e., the firstorder quadrupolar contribution). Conversely, in the solidecho spectrum, the \(\text{d}_I\)dependent frequency contributions (i.e., the firstorder quadrupolar contribution) were averaged to zero, leaving only the \(\text{p}_I\)dependent frequency contributions (i.e., the shielding contributions).
While these two examples nicely illustrate numerous important concepts for
building custom methods, it should also be noted that identical spectra could
have been obtained with a simpler custom method that used the freq_contrib
to remove the undesired frequency contributions. The code for these two methods
is illustrated below.
quad_only = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=2e4, # in Hz
events=[
SpectralEvent(
transition_queries=[{"ch1": {"P": [1]}}],
freq_contrib=["Quad1_2"]
)
],
)
],
)
shielding_only = Method(
channels=["2H"],
magnetic_flux_density=9.4, # in T
spectral_dimensions=[
SpectralDimension(
count=512,
spectral_width=2e4, # in Hz
events=[
SpectralEvent(
transition_queries=[{"ch1": {"P": [1]}}],
freq_contrib=["Shielding1_0", "Shielding1_2"],
)
],
)
],
)
sim = Simulator()
sim.spin_systems = [SpinSystem(sites=[deuterium])]
sim.methods = [quad_only, shielding_only]
sim.run()
processor = sp.SignalProcessor(
operations=[
sp.IFFT(),
sp.apodization.Gaussian(FWHM="100 Hz"),
sp.FFT(),
]
)
quad_only_dataset = processor.apply_operations(dataset=sim.methods[0].simulation)
shielding_only_dataset = processor.apply_operations(dataset=sim.methods[1].simulation)
fig, ax = plt.subplots(1, 2, subplot_kw={"projection": "csdm"}, figsize=[8.5, 3])
ax[0].set_title("Quad. Only Spectrum")
ax[0].plot(quad_only_dataset.real)
ax[0].invert_xaxis()
ax[0].grid()
ax[1].set_title("Shielding Only Spectrum")
ax[1].plot(shielding_only_dataset.real)
ax[1].invert_xaxis()
ax[1].grid()
plt.tight_layout()
plt.show()
Origin and Reference Offset¶
SpectralDimension()
has additional
attributes that have already been discussed in earlier sections of the documentation.
Notably, origin_offset
and reference_offset
are important for converting
the frequency coordinate into a dimensionless frequency ratio coordinate. For
spectra where all the spectral dimensions are associated with singlequantum
transitions on a single isotope, the convention for defining origin_offset
and reference_offset
is well established;
the origin_offset
, \(o_k\), is interpreted as the NMR spectrometer
frequency and the reference_offset
, \(b_k\), as the reference
frequency. Given the frequency coordinate, \({X}\), the corresponding
dimensionlessfrequency ratio follows,
In the case of multiple quantum dimensions, however, there appear
to be no formal conventions for defining origin_offset
and reference_offset
.
Attribute Summaries¶
Attribute Name 
Type 
Description 

channels 

A required list of isotopes given as strings over which the given method applies.
For example, 
magnetic_flux_density 

An optional float describing the macroscopic magnetic flux density of the applied
external magnetic field in tesla. For example, 
rotor_frequency 

An optional float describing the sample rotation frequency in Hz. For example, 
rotor_angle 

An optional float describing the angle between the sample rotation axis and the external
magnetic field in radians. The default value is the magic angle,

spectral_dimensions 

A list of SpectralDimension objects describing the spectral dimensions for the method. 
affine_matrix 

An optional ( 
simulation 
CSDM object 
A CSDM object representing the spectrum simulated by the method. By default, the value is

experiment 
CSDM object 
An optional CSDM object holding an experimental measurement of the method. The default
value is 
Attribute Name 
Type 
Description 

count 

An optional integer representing the number of points, \(N\), along the spectroscopic
dimension. For example, 
spectral_width 

An optional float representing the width, \(\Delta x\), of the spectroscopic dimension
in Hz. For example, 
reference_offset 

An optional float representing the reference offset, \(x_0\), of the spectroscopic
dimension in Hz. For example, 
origin_offset 

An optional float representing the origin offset, or Larmor frequency, along the
spectroscopic dimension in units of Hz. The default value is 
events 

An optional list of Events used to emulate an experiment. The default value is an empty list. 
Attribute Name 
Type 
Description 

magnetic_flux_density 

An optional float describing the macroscopic magnetic flux density of the applied
external magnetic field in tesla. For example, 
rotor_angle 

An optional float describing the angle between the sample rotation axis and the external
magnetic field in radians. The default is 
rotor_frequency 

An optional float describing the sample rotation frequency in Hz. For example, 
freq_contrib 

An optional list of FrequencyEnum (list of allowed strings) selecting which
contributions to include when calculating a transition frequency. For example,

transition_queries 

An optional 
Attribute Name 
Type 
Description 

query 

A 