¹³C MAS NMR of Glycine (CSA) [960 Hz]

The following is a sideband least-squares fitting example of a \(^{13}\text{C}\) MAS NMR spectrum of Glycine spinning at 960 Hz. The following experimental dataset is a part of DMFIT 1 examples. We thank Dr. Dominique Massiot for sharing the dataset.

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 BlochDecaySpectrum
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

host = "https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/"
filename = "13C MAS 960Hz - Glycine.csdf"
experiment = cp.load(host + filename)

# standard deviation of noise from the dataset
sigma = 3.822249

# 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.
plt.figure(figsize=(8, 4))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.set_xlim(280, -10)
plt.grid()
plt.tight_layout()
plt.show()

Create a fitting model

Spin System

C1 = Site(
    isotope="13C",
    isotropic_chemical_shift=176.0,  # in ppm
    shielding_symmetric={"zeta": 70, "eta": 0.6},  # zeta in Hz
)
C2 = Site(
    isotope="13C",
    isotropic_chemical_shift=43.0,  # in ppm
    shielding_symmetric={"zeta": 30, "eta": 0.5},  # zeta in Hz
)

spin_systems = [SpinSystem(sites=[C1], name="C1"), SpinSystem(sites=[C2], name="C2")]

Method

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

MAS = BlochDecaySpectrum(
    channels=["13C"],
    magnetic_flux_density=7.05,  # in T
    rotor_frequency=960,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=experiment,  # experimental dataset
)

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

Guess Model Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[MAS])
sim.config.decompose_spectrum = "spin_system"
sim.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Exponential(FWHM="20 Hz", dv_index=0),  # spin system 0
        sp.apodization.Exponential(FWHM="200 Hz", dv_index=1),  # spin system 1
        sp.FFT(),
        sp.Scale(factor=100),
    ]
)
processed_data = processor.apply_operations(data=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 4))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.plot(processed_data, linewidth=2, alpha=0.6)
ax.set_xlim(280, -10)
plt.grid()
plt.legend()
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, include={"rotor_frequency"})
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Exponential_FWHM            20     -inf      inf     True     None
SP_0_operation_2_Exponential_FWHM           200     -inf      inf     True     None
SP_0_operation_4_Scale_factor               100     -inf      inf     True     None
mth_0_rotor_frequency                       960      860     1060     True     None
sys_0_abundance                              50        0      100     True     None
sys_0_site_0_isotropic_chemical_shift       176     -inf      inf     True     None
sys_0_site_0_shielding_symmetric_eta        0.6        0        1     True     None
sys_0_site_0_shielding_symmetric_zeta        70     -inf      inf     True     None
sys_1_abundance                              50        0      100    False 100-sys_0_abundance
sys_1_site_0_isotropic_chemical_shift        43     -inf      inf     True     None
sys_1_site_0_shielding_symmetric_eta        0.5        0        1     True     None
sys_1_site_0_shielding_symmetric_zeta        30     -inf      inf     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   = 85
    # data points      = 4096
    # variables        = 11
    chi-square         = 4246.76287
    reduced chi-square = 1.03959923
    Akaike info crit   = 170.054533
    Bayesian info crit = 239.549961
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift:  176.131413 +/- 0.00127988 (0.00%) (init = 176)
    sys_0_site_0_shielding_symmetric_zeta:  71.2620768 +/- 0.25475162 (0.36%) (init = 70)
    sys_0_site_0_shielding_symmetric_eta:   0.91217925 +/- 0.00555631 (0.61%) (init = 0.6)
    sys_0_abundance:                        54.5029597 +/- 0.24448335 (0.45%) (init = 50)
    sys_1_site_0_isotropic_chemical_shift:  43.3284518 +/- 0.00887254 (0.02%) (init = 43)
    sys_1_site_0_shielding_symmetric_zeta:  23.0234065 +/- 0.18114229 (0.79%) (init = 30)
    sys_1_site_0_shielding_symmetric_eta:   0.50291614 +/- 0.04141365 (8.23%) (init = 0.5)
    sys_1_abundance:                        45.4970403 +/- 0.24448335 (0.54%) == '100-sys_0_abundance'
    mth_0_rotor_frequency:                  962.151219 +/- 0.04197724 (0.00%) (init = 960)
    SP_0_operation_1_Exponential_FWHM:      33.8471657 +/- 0.27755222 (0.82%) (init = 20)
    SP_0_operation_2_Exponential_FWHM:      179.799993 +/- 1.81583654 (1.01%) (init = 200)
    SP_0_operation_4_Scale_factor:          172.469361 +/- 0.82542295 (0.48%) (init = 100)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_1_site_0_shielding_symmetric_zeta, sys_1_site_0_shielding_symmetric_eta) = -0.601
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_site_0_shielding_symmetric_eta) = -0.448
    C(SP_0_operation_1_Exponential_FWHM, SP_0_operation_4_Scale_factor)            =  0.442
    C(sys_0_abundance, SP_0_operation_2_Exponential_FWHM)                          = -0.441
    C(SP_0_operation_2_Exponential_FWHM, SP_0_operation_4_Scale_factor)            =  0.409
    C(sys_0_abundance, SP_0_operation_1_Exponential_FWHM)                          =  0.398
    C(sys_0_abundance, SP_0_operation_4_Scale_factor)                              = -0.242
    C(sys_0_abundance, sys_1_site_0_shielding_symmetric_zeta)                      = -0.224
    C(sys_1_site_0_shielding_symmetric_zeta, SP_0_operation_4_Scale_factor)        =  0.209
    C(sys_0_site_0_isotropic_chemical_shift, SP_0_operation_1_Exponential_FWHM)    =  0.178
    C(sys_0_site_0_shielding_symmetric_zeta, SP_0_operation_4_Scale_factor)        =  0.149
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_abundance)                      =  0.136
    C(sys_0_abundance, sys_1_site_0_shielding_symmetric_eta)                       = -0.134
    C(sys_1_site_0_shielding_symmetric_eta, SP_0_operation_4_Scale_factor)         =  0.125

The best fit solution

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

plt.figure(figsize=(8, 4))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.plot(residuals, color="gray", linewidth=0.5, label="Residual")
ax.plot(best_fit, linewidth=2, alpha=0.6)
ax.set_xlim(280, -10)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()

1

D.Massiot, F.Fayon, M.Capron, I.King, S.Le Calvé, B.Alonso, J.O.Durand, B.Bujoli, Z.Gan, G.Hoatson, ‘Modelling one and two-dimensional solid-state NMR spectra.’, Magn. Reson. Chem. 40 70-76 (2002) DOI: 10.1002/mrc.984