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

from mrsimulator import Simulator, SpinSystem, Site
from mrsimulator.methods import SSB2D
from mrsimulator import signal_processing as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions

Import the dataset

filename = "https://sandbox.zenodo.org/record/834704/files/Rb2SO4_QMAT.csdf"
qmat_data = 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_data = qmat_data.real

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

# plot of the dataset.
max_amp = qmat_data.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_data.T, colors="k", **options)
ax.set_xlim(200, -200)
ax.set_ylim(75, -120)
plt.grid()
plt.tight_layout()
plt.show()

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={"Cq": 5.5e6, "eta": 0.1},  # Cq in Hz
)
Rb_2 = Site(
    isotope="87Rb",
    isotropic_chemical_shift=40,  # in ppm
    quadrupolar={"Cq": 2.1e6, "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_data)

PASS = SSB2D(
    channels=["87Rb"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=2604,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=qmat_data,  # 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(axis=0),
        sp.apodization.Gaussian(FWHM="50 Hz"),
        sp.IFFT(axis=0),
        sp.Scale(factor=1e4),
    ]
)
processed_data = processor.apply_operations(data=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(qmat_data.T, colors="k", **options)
ax.contour(processed_data.T, colors="r", linestyles="--", **options)
ax.set_xlim(200, -200)
ax.set_ylim(75, -120)
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               50     -inf      inf     True     None
SP_0_operation_3_Scale_factor             1e+04     -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.5e+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.1e+06     -inf      inf     True     None
sys_1_site_0_quadrupolar_eta               0.95        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   = 133
    # data points      = 65536
    # variables        = 9
    chi-square         = 217831.432
    reduced chi-square = 3.32430040
    Akaike info crit   = 78734.7253
    Bayesian info crit = 78816.5385
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift:  13.8247093 +/- 0.01883732 (0.14%) (init = 16)
    sys_0_site_0_quadrupolar_Cq:            5229369.47 +/- 989.315964 (0.02%) (init = 5500000)
    sys_0_site_0_quadrupolar_eta:           0.12875964 +/- 3.8584e-04 (0.30%) (init = 0.1)
    sys_0_abundance:                        55.8944405 +/- 0.05573973 (0.10%) (init = 50)
    sys_1_site_0_isotropic_chemical_shift:  40.4488902 +/- 0.00658995 (0.02%) (init = 40)
    sys_1_site_0_quadrupolar_Cq:            2697016.78 +/- 996.758989 (0.04%) (init = 2100000)
    sys_1_site_0_quadrupolar_eta:           0.86990373 +/- 0.00143241 (0.16%) (init = 0.95)
    sys_1_abundance:                        44.1055595 +/- 0.05573973 (0.13%) == '100-sys_0_abundance'
    SP_0_operation_1_Gaussian_FWHM:        -142.754756 +/- 2.83847927 (1.99%) (init = 50)
    SP_0_operation_3_Scale_factor:          6277.07909 +/- 7.75261506 (0.12%) (init = 10000)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_Cq)    =  0.827
    C(sys_1_site_0_quadrupolar_Cq, sys_1_site_0_quadrupolar_eta)             = -0.818
    C(sys_1_site_0_isotropic_chemical_shift, sys_1_site_0_quadrupolar_Cq)    =  0.633
    C(sys_0_abundance, SP_0_operation_3_Scale_factor)                        =  0.630
    C(sys_1_site_0_quadrupolar_eta, SP_0_operation_1_Gaussian_FWHM)          = -0.356
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_eta)   = -0.263
    C(sys_0_site_0_quadrupolar_eta, SP_0_operation_3_Scale_factor)           =  0.240
    C(sys_0_site_0_quadrupolar_eta, sys_0_abundance)                         =  0.235
    C(sys_1_site_0_isotropic_chemical_shift, sys_1_site_0_quadrupolar_eta)   = -0.198
    C(SP_0_operation_1_Gaussian_FWHM, SP_0_operation_3_Scale_factor)         = -0.197
    C(sys_0_site_0_quadrupolar_Cq, sys_0_site_0_quadrupolar_eta)             = -0.191
    C(sys_0_abundance, sys_1_site_0_isotropic_chemical_shift)                = -0.178
    C(sys_0_site_0_quadrupolar_Cq, SP_0_operation_3_Scale_factor)            =  0.164
    C(sys_1_site_0_isotropic_chemical_shift, SP_0_operation_3_Scale_factor)  =  0.162
    C(sys_0_abundance, sys_1_site_0_quadrupolar_Cq)                          = -0.157
    C(sys_1_site_0_quadrupolar_Cq, SP_0_operation_1_Gaussian_FWHM)           =  0.156
    C(sys_0_site_0_quadrupolar_Cq, sys_0_abundance)                          =  0.138
    C(sys_1_site_0_isotropic_chemical_shift, SP_0_operation_1_Gaussian_FWHM) = -0.124
    C(sys_1_site_0_quadrupolar_Cq, SP_0_operation_3_Scale_factor)            =  0.104
    C(sys_0_site_0_isotropic_chemical_shift, SP_0_operation_3_Scale_factor)  =  0.104

The best fit solution

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

# Plot of the best fit solution
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(qmat_data.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()

1

B. J. Walder, K. K. Dey, D. C. Kaseman, J. H. Baltisberger, and P. J. Grandinetti, Sideband separation experiments in NMR with phase incremented echo train acquisition, J. Phys. Chem. 2013, 138, 174203-1-12. DOI: 10.1063/1.4803142