⁸⁷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, report_fit

from mrsimulator import Simulator
from mrsimulator.methods import ThreeQ_VAS
from mrsimulator import signal_processing 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://sandbox.zenodo.org/record/814455/files/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()

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(
    isotopes="87Rb",
    isotropic_chemical_shifts=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_data = processor.apply_operations(data=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_data, colors="r", linestyles="--", **options)
ax.set_xlim(-20, -50)
ax.set_ylim(-45, -65)
plt.grid()
plt.tight_layout()
plt.show()

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

minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma))
result = minner.minimize()
report_fit(result)

Out:

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 409
    # data points      = 262144
    # variables        = 14
    chi-square         = 514378.952
    reduced chi-square = 1.96230478
    Akaike info crit   = 176730.432
    Bayesian info crit = 176877.105
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift: -26.7600498 +/- 3.3639e-04 (0.00%) (init = -26.4)
    sys_0_site_0_quadrupolar_Cq:            1693623.74 +/- 115.012581 (0.01%) (init = 1700000)
    sys_0_site_0_quadrupolar_eta:           0.21817910 +/- 3.1396e-04 (0.14%) (init = 0.2)
    sys_0_abundance:                        41.9395333 +/- 0.02257962 (0.05%) (init = 33.33333)
    sys_1_site_0_isotropic_chemical_shift: -28.1157726 +/- 9.7428e-04 (0.00%) (init = -28.5)
    sys_1_site_0_quadrupolar_Cq:            2027908.41 +/- 348.654106 (0.02%) (init = 2000000)
    sys_1_site_0_quadrupolar_eta:           0.90003177 +/- 7.1639e-04 (0.08%) (init = 1)
    sys_1_abundance:                        24.2721929 +/- 0.02577385 (0.11%) (init = 33.33333)
    sys_2_site_0_isotropic_chemical_shift: -31.0314607 +/- 4.1538e-04 (0.00%) (init = -31.3)
    sys_2_site_0_quadrupolar_Cq:            1721336.41 +/- 151.036797 (0.01%) (init = 1700000)
    sys_2_site_0_quadrupolar_eta:           0.56032125 +/- 4.5876e-04 (0.08%) (init = 0.6)
    sys_2_abundance:                        33.7882738 +/- 0.02119368 (0.06%) == '100-sys_0_abundance-sys_1_abundance'
    SP_0_operation_1_Gaussian_FWHM:         0.06201469 +/- 2.8251e-04 (0.46%) (init = 0.08)
    SP_0_operation_2_Gaussian_FWHM:         0.15110281 +/- 9.3281e-05 (0.06%) (init = 0.1)
    SP_0_operation_4_Scale_factor:          13385739.8 +/- 6957.71769 (0.05%) (init = 1e+07)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_1_site_0_quadrupolar_Cq, sys_1_site_0_quadrupolar_eta)             = -0.805
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_Cq)    = -0.795
    C(sys_2_site_0_quadrupolar_Cq, sys_2_site_0_quadrupolar_eta)             = -0.646
    C(sys_0_abundance, sys_1_abundance)                               = -0.623
    C(SP_0_operation_2_Gaussian_FWHM, SP_0_operation_4_Scale_factor)         =  0.452
    C(sys_2_site_0_isotropic_chemical_shift, sys_2_site_0_quadrupolar_Cq)    = -0.415
    C(sys_1_site_0_isotropic_chemical_shift, sys_1_site_0_quadrupolar_eta)   = -0.293
    C(sys_0_site_0_quadrupolar_Cq, sys_0_site_0_quadrupolar_eta)             = -0.284
    C(sys_1_site_0_isotropic_chemical_shift, sys_1_site_0_quadrupolar_Cq)    = -0.276
    C(sys_0_abundance, SP_0_operation_4_Scale_factor)                        = -0.274
    C(sys_2_site_0_isotropic_chemical_shift, sys_2_site_0_quadrupolar_eta)   = -0.272
    C(sys_1_abundance, SP_0_operation_4_Scale_factor)                        =  0.245
    C(sys_0_site_0_quadrupolar_eta, SP_0_operation_1_Gaussian_FWHM)          = -0.238
    C(sys_1_site_0_quadrupolar_eta, sys_1_abundance)                         = -0.228
    C(sys_1_site_0_quadrupolar_Cq, sys_1_abundance)                          =  0.228
    C(sys_0_abundance, sys_1_site_0_quadrupolar_eta)                         =  0.176
    C(sys_0_abundance, sys_1_site_0_quadrupolar_Cq)                          = -0.160
    C(sys_1_abundance, SP_0_operation_2_Gaussian_FWHM)                       = -0.145
    C(sys_2_site_0_quadrupolar_Cq, SP_0_operation_1_Gaussian_FWHM)           = -0.130
    C(sys_0_site_0_quadrupolar_eta, sys_0_abundance)                         =  0.117
    C(sys_1_site_0_quadrupolar_eta, SP_0_operation_2_Gaussian_FWHM)          =  0.114
    C(SP_0_operation_1_Gaussian_FWHM, SP_0_operation_4_Scale_factor)         =  0.107
    C(sys_2_site_0_isotropic_chemical_shift, SP_0_operation_1_Gaussian_FWHM) =  0.106
    C(sys_1_site_0_quadrupolar_Cq, SP_0_operation_4_Scale_factor)            =  0.102

The best fit solution

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

# 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()

Image plots with residuals

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

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()