¹³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

from mrsimulator import Simulator, SpinSystem, Site
from mrsimulator.method.lib import BlochDecaySpectrum
from mrsimulator import signal_processor 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.982936

# 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()
plot 2 13C glycine 960Hz

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_dataset = processor.apply_operations(dataset=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_dataset, linewidth=2, alpha=0.6)
ax.set_xlim(280, -10)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plot 2 13C glycine 960Hz

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

Fit Statistics

fitting methodleastsq
# function evals85
# data points4096
# variables11
chi-square 3911.01342
reduced chi-square 0.95740843
Akaike info crit.-167.294017
Bayesian info crit.-97.7985895

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift 176.131431 0.00127988 (0.00%) 176.0 -inf inf True
sys_0_site_0_shielding_symmetric_zeta 71.2620847 0.25475165 (0.36%) 70.0 -inf inf True
sys_0_site_0_shielding_symmetric_eta 0.91217923 0.00555631 (0.61%) 0.6 0.00000000 1.00000000 True
sys_0_abundance 54.5029599 0.24448158 (0.45%) 50.0 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift 43.3284564 0.00887254 (0.02%) 43.0 -inf inf True
sys_1_site_0_shielding_symmetric_zeta 23.0234103 0.18113487 (0.79%) 30.0 -inf inf True
sys_1_site_0_shielding_symmetric_eta 0.50291564 0.04141684 (8.24%) 0.5 0.00000000 1.00000000 True
sys_1_abundance 45.4970401 0.24448158 (0.54%) 50.0 0.00000000 100.000000 False 100-sys_0_abundance
mth_0_rotor_frequency 962.151218 0.04197685 (0.00%) 960.0 860.000000 1060.00000 True
SP_0_operation_1_Exponential_FWHM 33.8471657 0.27755222 (0.82%) 20.0 -inf inf True
SP_0_operation_2_Exponential_FWHM 179.799996 1.81583677 (1.01%) 200.0 -inf inf True
SP_0_operation_4_Scale_factor 172.469360 0.82541789 (0.48%) 100.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_1_site_0_shielding_symmetric_zetasys_1_site_0_shielding_symmetric_eta-0.6011
sys_0_site_0_shielding_symmetric_zetasys_0_site_0_shielding_symmetric_eta-0.4482
SP_0_operation_1_Exponential_FWHMSP_0_operation_4_Scale_factor0.4423
sys_0_abundanceSP_0_operation_2_Exponential_FWHM-0.4409
SP_0_operation_2_Exponential_FWHMSP_0_operation_4_Scale_factor0.4092
sys_0_abundanceSP_0_operation_1_Exponential_FWHM0.3976
sys_0_abundanceSP_0_operation_4_Scale_factor-0.2420
sys_0_abundancesys_1_site_0_shielding_symmetric_zeta-0.2238
sys_1_site_0_shielding_symmetric_zetaSP_0_operation_4_Scale_factor0.2090
sys_0_site_0_isotropic_chemical_shiftSP_0_operation_1_Exponential_FWHM0.1780
sys_0_site_0_shielding_symmetric_zetaSP_0_operation_4_Scale_factor0.1486
sys_0_site_0_shielding_symmetric_zetasys_0_abundance0.1361
sys_0_abundancesys_1_site_0_shielding_symmetric_eta-0.1339
sys_1_site_0_shielding_symmetric_etaSP_0_operation_4_Scale_factor0.1247


The best fit solution

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

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()
plot 2 13C glycine 960Hz

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

Gallery generated by Sphinx-Gallery