⁸⁷Rb 2D 3QMAS NMR of RbNO₃

The following is a 3QMAS fitting example for \(\text{RbNO}_3\). The dataset was acquired and shared by Brendan Wilson.

import numpy as np
import csdmpy as cp
import matplotlib.pyplot as plt
from lmfit import Minimizer

from mrsimulator import Simulator
from mrsimulator.method.lib import ThreeQ_VAS
from mrsimulator import signal_processor as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions
from mrsimulator.utils.collection import single_site_system_generator

Import the dataset

filename = "https://ssnmr.org/sites/default/files/mrsimulator/RbNO3_MQMAS.csdf"
experiment = cp.load(filename)

# standard deviation of noise from the dataset
sigma = 175.5476

# For spectral fitting, we only focus on the real part of the complex dataset
experiment = experiment.real

# Convert the coordinates along each dimension from Hz to ppm.
_ = [item.to("ppm", "nmr_frequency_ratio") for item in experiment.dimensions]

# plot of the dataset.
max_amp = experiment.max()
levels = (np.arange(24) + 1) * max_amp / 25  # contours are drawn at these levels.
options = dict(levels=levels, alpha=0.75, linewidths=0.5)  # plot options

plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.contour(experiment, colors="k", **options)
ax.set_xlim(-20, -50)
ax.set_ylim(-45, -65)
plt.grid()
plt.tight_layout()
plt.show()
plot 3 RbNO3 MQMAS

Create a fitting model

Guess model

Create a guess list of spin systems.

shifts = [-26.4, -28.5, -31.3]  # in ppm
Cq = [1.7e6, 2.0e6, 1.7e6]  # in  Hz
eta = [0.2, 1.0, 0.6]
abundance = [33.33, 33.33, 33.33]  # in %

spin_systems = single_site_system_generator(
    isotope="87Rb",
    isotropic_chemical_shift=shifts,
    quadrupolar={"Cq": Cq, "eta": eta},
    abundance=abundance,
)

Method

Create the 3QMAS method.

# Get the spectral dimension parameters from the experiment.
spectral_dims = get_spectral_dimensions(experiment)

MQMAS = ThreeQ_VAS(
    channels=["87Rb"],
    magnetic_flux_density=9.395,  # in T
    spectral_dimensions=spectral_dims,
    experiment=experiment,  # add the measurement to the method.
)

# Optimize the script by pre-setting the transition pathways for each spin system from
# the das method.
for sys in spin_systems:
    sys.transition_pathways = MQMAS.get_transition_pathways(sys)

Guess Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[MQMAS])
sim.config.number_of_sidebands = 1
sim.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        # Gaussian convolution along both dimensions.
        sp.IFFT(dim_index=(0, 1)),
        sp.apodization.Gaussian(FWHM="0.08 kHz", dim_index=0),
        sp.apodization.Gaussian(FWHM="0.1 kHz", dim_index=1),
        sp.FFT(dim_index=(0, 1)),
        sp.Scale(factor=1e7),
    ]
)
processed_dataset = processor.apply_operations(dataset=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.contour(experiment, colors="k", **options)
ax.contour(processed_dataset, colors="r", linestyles="--", **options)
ax.set_xlim(-20, -50)
ax.set_ylim(-45, -65)
plt.grid()
plt.tight_layout()
plt.show()
plot 3 RbNO3 MQMAS

Least-squares minimization with LMFIT

Use the make_LMFIT_params() for a quick setup of the fitting parameters.

params = sf.make_LMFIT_params(sim, processor)
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM             0.08     -inf      inf     True     None
SP_0_operation_2_Gaussian_FWHM              0.1     -inf      inf     True     None
SP_0_operation_4_Scale_factor             1e+07     -inf      inf     True     None
sys_0_abundance                           33.33        0      100     True     None
sys_0_site_0_isotropic_chemical_shift     -26.4     -inf      inf     True     None
sys_0_site_0_quadrupolar_Cq             1.7e+06     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta                0.2        0        1     True     None
sys_1_abundance                           33.33        0      100     True     None
sys_1_site_0_isotropic_chemical_shift     -28.5     -inf      inf     True     None
sys_1_site_0_quadrupolar_Cq               2e+06     -inf      inf     True     None
sys_1_site_0_quadrupolar_eta                  1        0        1     True     None
sys_2_abundance                           33.33        0      100    False 100-sys_0_abundance-sys_1_abundance
sys_2_site_0_isotropic_chemical_shift     -31.3     -inf      inf     True     None
sys_2_site_0_quadrupolar_Cq             1.7e+06     -inf      inf     True     None
sys_2_site_0_quadrupolar_eta                0.6        0        1     True     None
None

Solve the minimizer using LMFIT

Fit Statistics

fitting methodleastsq
# function evals287
# data points262144
# variables14
chi-square 514166.208
reduced chi-square 1.96149318
Akaike info crit. 176621.988
Bayesian info crit. 176768.661

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift -26.7425947 3.3482e-04 (0.00%) -26.4 -inf inf True
sys_0_site_0_quadrupolar_Cq 1694125.71 114.895763 (0.01%) 1700000.0 -inf inf True
sys_0_site_0_quadrupolar_eta 0.21788122 3.1374e-04 (0.14%) 0.2 0.00000000 1.00000000 True
sys_0_abundance 41.9246619 0.02246211 (0.05%) 33.333333333333336 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift -28.1022004 0.00103119 (0.00%) -28.5 -inf inf True
sys_1_site_0_quadrupolar_Cq 2026012.18 416.618303 (0.02%) 2000000.0 -inf inf True
sys_1_site_0_quadrupolar_eta 0.90758492 8.5980e-04 (0.09%) 1.0 0.00000000 1.00000000 True
sys_1_abundance 24.2906054 0.02570929 (0.11%) 33.333333333333336 0.00000000 100.000000 True
sys_2_site_0_isotropic_chemical_shift -31.0115064 4.1666e-04 (0.00%) -31.3 -inf inf True
sys_2_site_0_quadrupolar_Cq 1721914.37 151.906451 (0.01%) 1700000.0 -inf inf True
sys_2_site_0_quadrupolar_eta 0.56031029 4.5970e-04 (0.08%) 0.6 0.00000000 1.00000000 True
sys_2_abundance 33.7847327 0.02136581 (0.06%) 33.33333333333332 0.00000000 100.000000 False 100-sys_0_abundance-sys_1_abundance
SP_0_operation_1_Gaussian_FWHM 0.06205343 2.8075e-04 (0.45%) 0.08 -inf inf True
SP_0_operation_2_Gaussian_FWHM 0.15097947 9.3595e-05 (0.06%) 0.1 -inf inf True
SP_0_operation_4_Scale_factor 13383848.1 6961.20143 (0.05%) 10000000.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_1_site_0_quadrupolar_Cqsys_1_site_0_quadrupolar_eta-0.8497
sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_quadrupolar_Cq-0.7934
sys_2_site_0_quadrupolar_Cqsys_2_site_0_quadrupolar_eta-0.6463
sys_0_abundancesys_1_abundance-0.6139
SP_0_operation_2_Gaussian_FWHMSP_0_operation_4_Scale_factor0.4577
sys_2_site_0_isotropic_chemical_shiftsys_2_site_0_quadrupolar_Cq-0.4173
sys_0_site_0_quadrupolar_Cqsys_0_site_0_quadrupolar_eta-0.2966
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_eta-0.2849
sys_0_abundanceSP_0_operation_4_Scale_factor-0.2756
sys_2_site_0_isotropic_chemical_shiftsys_2_site_0_quadrupolar_eta-0.2687
sys_1_abundanceSP_0_operation_4_Scale_factor0.2478
sys_0_site_0_quadrupolar_etaSP_0_operation_1_Gaussian_FWHM-0.2232
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_Cq-0.2162
sys_1_abundanceSP_0_operation_2_Gaussian_FWHM-0.1318
sys_2_site_0_quadrupolar_CqSP_0_operation_1_Gaussian_FWHM-0.1303
sys_1_site_0_quadrupolar_Cqsys_1_abundance0.1204
sys_0_site_0_quadrupolar_etasys_0_abundance0.1127
sys_2_site_0_isotropic_chemical_shiftSP_0_operation_1_Gaussian_FWHM0.1066
SP_0_operation_1_Gaussian_FWHMSP_0_operation_4_Scale_factor0.1064


The best fit solution

best_fit = sf.bestfit(sim, processor)[0].real

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.contour(experiment, colors="k", **options)
ax.contour(best_fit, colors="r", linestyles="--", **options)
ax.set_xlim(-20, -50)
ax.set_ylim(-45, -65)
plt.grid()
plt.tight_layout()
plt.show()
plot 3 RbNO3 MQMAS

Image plots with residuals

residuals = sf.residuals(sim, processor)[0].real

fig, ax = plt.subplots(
    1, 3, sharey=True, figsize=(10, 3.0), subplot_kw={"projection": "csdm"}
)
vmax, vmin = experiment.max(), experiment.min()
for i, dat in enumerate([experiment, best_fit, residuals]):
    ax[i].imshow(dat, aspect="auto", cmap="gist_ncar_r", vmax=vmax, vmin=vmin)
    ax[i].set_xlim(-20, -50)
ax[0].set_ylim(-45, -65)
plt.tight_layout()
plt.show()
plot 3 RbNO3 MQMAS

Total running time of the script: ( 0 minutes 19.617 seconds)

Gallery generated by Sphinx-Gallery