⁸⁷Rb 2D QMAT NMR of Rb₂SO₄

The following is an illustration for fitting 2D QMAT/QPASS datasets. The example dataset is a \(^{87}\text{Rb}\) 2D QMAT spectrum of \(\text{Rb}_2\text{SO}_4\) from Walder et al. [1]

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

from mrsimulator import Simulator, SpinSystem, Site
from mrsimulator.method.lib import SSB2D
from mrsimulator import signal_processor as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions
from mrsimulator.spin_system.tensors import SymmetricTensor

Import the dataset

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

# standard deviation of noise from the dataset
sigma = 6.530634

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

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

# plot of the dataset.
max_amp = qmat_dataset.max()
levels = (np.arange(31) + 0.15) * max_amp / 32  # contours are drawn at these levels.
options = dict(levels=levels, alpha=1, linewidths=0.5)  # plot options

plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(qmat_dataset.T, colors="k", **options)
ax.set_xlim(200, -200)
ax.set_ylim(75, -120)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 Rb2SO4 QMAT

Create a fitting model

Guess model

Create a guess list of spin systems.

Rb_1 = Site(
    isotope="87Rb",
    isotropic_chemical_shift=16,  # in ppm
    quadrupolar=SymmetricTensor(Cq=5.3e6, eta=0.1),  # Cq in Hz
)
Rb_2 = Site(
    isotope="87Rb",
    isotropic_chemical_shift=40,  # in ppm
    quadrupolar=SymmetricTensor(Cq=2.2e6, eta=0.95),  # Cq in Hz
)

spin_systems = [SpinSystem(sites=[s]) for s in [Rb_1, Rb_2]]

Method

Create the SSB2D method.

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

PASS = SSB2D(
    channels=["87Rb"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=2604,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=qmat_dataset,  # add the measurement to the method.
)

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

Guess Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[PASS])
sim.run()

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

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(qmat_dataset.T, colors="k", **options)
ax.contour(processed_dataset.T, colors="r", linestyles="--", **options)
ax.set_xlim(200, -200)
ax.set_ylim(75, -120)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 Rb2SO4 QMAT

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)
params["SP_0_operation_1_Gaussian_FWHM"].min = 0
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM              100        0      inf     True     None
SP_0_operation_3_Scale_factor             1e+07     -inf      inf     True     None
sys_0_abundance                              50        0      100     True     None
sys_0_site_0_isotropic_chemical_shift        16     -inf      inf     True     None
sys_0_site_0_quadrupolar_Cq             5.3e+06     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta                0.1        0        1     True     None
sys_1_abundance                              50        0      100    False 100-sys_0_abundance
sys_1_site_0_isotropic_chemical_shift        40     -inf      inf     True     None
sys_1_site_0_quadrupolar_Cq             2.2e+06     -inf      inf     True     None
sys_1_site_0_quadrupolar_eta               0.95        0        1     True     None
None

Solve the minimizer using LMFIT

Fit Statistics

fitting methodleastsq
# function evals132
# data points65536
# variables9
chi-square 217831.436
reduced chi-square 3.32430045
Akaike info crit. 78734.7264
Bayesian info crit. 78816.5396

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift 13.8163037 0.01882009 (0.14%) 16.0 -inf inf True
sys_0_site_0_quadrupolar_Cq 5231087.06 989.456999 (0.02%) 5300000.0 -inf inf True
sys_0_site_0_quadrupolar_eta 0.12875115 3.8569e-04 (0.30%) 0.1 0.00000000 1.00000000 True
sys_0_abundance 55.8942185 0.05573846 (0.10%) 50.0 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift 40.4226839 0.00658482 (0.02%) 40.0 -inf inf True
sys_1_site_0_quadrupolar_Cq 2697888.13 996.535912 (0.04%) 2200000.0 -inf inf True
sys_1_site_0_quadrupolar_eta 0.86991599 0.00143127 (0.16%) 0.95 0.00000000 1.00000000 True
sys_1_abundance 44.1057815 0.05573846 (0.13%) 50.0 0.00000000 100.000000 False 100-sys_0_abundance
SP_0_operation_1_Gaussian_FWHM 142.791704 2.83967727 (1.99%) 100.0 0.00000000 inf True
SP_0_operation_3_Scale_factor 25710912.1 31754.3563 (0.12%) 10000000.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_quadrupolar_Cq0.8275
sys_1_site_0_quadrupolar_Cqsys_1_site_0_quadrupolar_eta-0.8178
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_Cq0.6335
sys_0_abundanceSP_0_operation_3_Scale_factor0.6296
sys_1_site_0_quadrupolar_etaSP_0_operation_1_Gaussian_FWHM0.3560
sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_quadrupolar_eta-0.2628
sys_0_site_0_quadrupolar_etaSP_0_operation_3_Scale_factor0.2396
sys_0_site_0_quadrupolar_etasys_0_abundance0.2347
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_eta-0.1984
SP_0_operation_1_Gaussian_FWHMSP_0_operation_3_Scale_factor0.1967
sys_0_site_0_quadrupolar_Cqsys_0_site_0_quadrupolar_eta-0.1908
sys_0_abundancesys_1_site_0_isotropic_chemical_shift-0.1782
sys_0_site_0_quadrupolar_CqSP_0_operation_3_Scale_factor0.1639
sys_1_site_0_isotropic_chemical_shiftSP_0_operation_3_Scale_factor0.1618
sys_0_abundancesys_1_site_0_quadrupolar_Cq-0.1567
sys_1_site_0_quadrupolar_CqSP_0_operation_1_Gaussian_FWHM-0.1560
sys_0_site_0_quadrupolar_Cqsys_0_abundance0.1376
sys_1_site_0_isotropic_chemical_shiftSP_0_operation_1_Gaussian_FWHM0.1242
sys_1_site_0_quadrupolar_CqSP_0_operation_3_Scale_factor0.1041
sys_0_site_0_isotropic_chemical_shiftSP_0_operation_3_Scale_factor0.1041


The best fit solution

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

# Plot of the best fit solution
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(qmat_dataset.T, colors="k", **options)
ax.contour(best_fit.T, colors="r", linestyles="--", **options)
ax.set_xlim(200, -200)
ax.set_ylim(75, -120)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 Rb2SO4 QMAT

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

Gallery generated by Sphinx-Gallery