Note
Click here to download the full example code
Amorphous material, ²⁹Si (I=1/2)¶
²⁹Si (I=1/2) simulation of amorphous material.
One of the advantages of the mrsimulator package is that it is a fast NMR spectrum simulation library. We can exploit this feature to simulate bulk spectra and eventually model amorphous materials. In this section, we illustrate how the mrsimulator library may be used in simulating the NMR spectrum of amorphous materials.
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mrsimulator import Simulator
from mrsimulator.method.lib import BlochDecaySpectrum
from mrsimulator.utils.collection import single_site_system_generator
from mrsimulator.method import SpectralDimension
Generating tensor parameter distribution¶
We model the amorphous material by assuming a distribution of interaction tensors. For example, a tri-variate normal distribution of the shielding tensor parameters, i.e., the isotropic chemical shift, the anisotropy parameter, \(\zeta\), and the asymmetry parameter, \(\eta\). In the following, we use pure NumPy and SciPy methods to generate the three-dimensional distribution, as follows,
mean = [-100, 50, 0.15] # given as [isotropic chemical shift in ppm, zeta in ppm, eta].
covariance = [[3.25, 0, 0], [0, 26.2, 0], [0, 0, 0.002]] # same order as the mean.
# range of coordinates along the three dimensions
iso_range = np.arange(100) * 0.3055 - 115 # in ppm
zeta_range = np.arange(30) * 2.5 + 10 # in ppm
eta_range = np.arange(21) / 20
# The coordinates grid
iso, zeta, eta = np.meshgrid(iso_range, zeta_range, eta_range, indexing="ij")
pos = np.asarray([iso, zeta, eta]).T
# Three-dimensional probability distribution function.
pdf = multivariate_normal(mean=mean, cov=covariance).pdf(pos).T
Here, iso
, zeta
, and eta
are the isotropic chemical shift, nuclear
shielding anisotropy, and nuclear shielding asymmetry coordinates of the 3D-grid
system over which the multivariate normal probability distribution is evaluated. The
mean of the distribution is given by the variable mean
and holds a value of -100
ppm, 50 ppm, and 0.15 for the isotropic chemical shift, nuclear shielding anisotropy,
and nuclear shielding asymmetry parameter, respectively. Similarly, the variable
covariance
holds the covariance matrix of the multivariate normal distribution.
The two-dimensional projections from this three-dimensional distribution are shown
below.
_, ax = plt.subplots(1, 3, figsize=(9, 3))
# isotropic shift v.s. shielding anisotropy
ax[0].contourf(zeta_range, iso_range, pdf.sum(axis=2))
ax[0].set_xlabel(r"shielding anisotropy, $\zeta$ / ppm")
ax[0].set_ylabel("isotropic chemical shift / ppm")
# isotropic shift v.s. shielding asymmetry
ax[1].contourf(eta_range, iso_range, pdf.sum(axis=1))
ax[1].set_xlabel(r"shielding asymmetry, $\eta$")
ax[1].set_ylabel("isotropic chemical shift / ppm")
# shielding anisotropy v.s. shielding asymmetry
ax[2].contourf(eta_range, zeta_range, pdf.sum(axis=0))
ax[2].set_xlabel(r"shielding asymmetry, $\eta$")
ax[2].set_ylabel(r"shielding anisotropy, $\zeta$ / ppm")
plt.tight_layout()
plt.show()
Create the Simulator object¶
Spin system:
Let’s create the sites and single-site spin system objects from these parameters.
Use the single_site_system_generator()
utility
function to generate single-site spin systems. # Here, iso
, zeta
, and eta
are the array of tensor parameter coordinates, and pdf
is the array of the
corresponding amplitudes.
spin_systems = single_site_system_generator(
isotope="29Si",
isotropic_chemical_shift=iso,
shielding_symmetric={"zeta": zeta, "eta": eta},
abundance=pdf,
)
Method:
Let’s also create a Bloch decay spectrum method.
method = BlochDecaySpectrum(
channels=["29Si"],
spectral_dimensions=[
SpectralDimension(spectral_width=25000, reference_offset=-7000) # values in Hz
],
)
The above method simulates a static \(^{29}\text{Si}\) spectrum at 9.4 T field (default value).
Simulator:
Now that we have the spin systems and the method, create the simulator object and add the respective objects.
Static spectrum¶
Observe the static \(^{29}\text{Si}\) NMR spectrum simulation.
sim.run()
The plot of the simulation.
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(sim.methods[0].simulation.real, color="black", linewidth=1)
ax.invert_xaxis()
plt.tight_layout()
plt.show()
Note
The broad spectrum seen in the above figure is a result of spectral averaging of spectra arising from a distribution of shielding tensors. There is no line-broadening filter applied to the spectrum.
Spinning sideband simulation at \(90^\circ\)¶
Here is an example of a sideband simulation, spinning at a 90-degree angle.
sim.methods[0] = BlochDecaySpectrum(
channels=["29Si"],
rotor_frequency=5000, # in Hz
rotor_angle=1.57079, # in rads, equivalent to 90 deg.
spectral_dimensions=[
SpectralDimension(spectral_width=25000, reference_offset=-7000) # values in Hz
],
)
sim.config.number_of_sidebands = 8 # eight sidebands are sufficient for this example
sim.run()
The plot of the simulation.
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(sim.methods[0].simulation.real, color="black", linewidth=1)
ax.invert_xaxis()
plt.tight_layout()
plt.show()
Spinning sideband simulation at the magic angle¶
Here is another example of a sideband simulation at the magic angle.
sim.methods[0] = BlochDecaySpectrum(
channels=["29Si"],
rotor_frequency=1000, # in Hz
rotor_angle=54.735 * np.pi / 180.0, # in rads
spectral_dimensions=[
SpectralDimension(spectral_width=25000, reference_offset=-7000) # values in Hz
],
)
sim.config.number_of_sidebands = 16 # sixteen sidebands are sufficient for this example
sim.run()
The plot of the simulation.
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(sim.methods[0].simulation.real, color="black", linewidth=1)
ax.invert_xaxis()
plt.tight_layout()
plt.show()
Total running time of the script: ( 0 minutes 12.467 seconds)