¹³C MAS NMR of Glycine (CSA) multi-spectra fit

The following is a multi-dataset least-squares fitting example of \(^{13}\text{C}\) MAS NMR spectrum of Glycine spinning at 5 kHz, 1.94 kHz, and 960 Hz. Before trying multi-dataset fitting, we recommend that you first try individual fits. The experimental datasets are part of DMFIT [1] examples.

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
from mrsimulator.spin_system.tensors import SymmetricTensor

Import the datasets

Import the datasets and assign the standard deviation of noise for each dataset. Here, sigma1, sigma2, and sigma3 are the noise standard deviation for the dataset acquired at 5 kHz, 1.94 kHz, and 960 Hz spinning speeds, respectively.

host = "https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/"
filename1 = "13C MAS 5000Hz - Glycine.csdf"
filename2 = "13C MAS 1940Hz - Glycine.csdf"
filename3 = "13C MAS 960Hz - Glycine.csdf"

experiment1 = cp.load(host + filename1).real
experiment2 = cp.load(host + filename2).real
experiment3 = cp.load(host + filename3).real
experiments = [experiment1, experiment2, experiment3]

# standard deviation of noise from the dataset
sigma1 = 1.97637
sigma2 = 3.822249
sigma3 = 3.982936
sigmas = [sigma1, sigma2, sigma3]

# Convert the coordinates along each dimension from Hz to ppm for each dataset
fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"})
for i, experiment in enumerate(experiments):
    _ = [item.to("ppm", "nmr_frequency_ratio") for item in experiment.dimensions]

    # plot of the dataset.
    ax[i].plot(experiment, color="black", linewidth=0.5, label="Experiment")
    ax[i].set_xlim(280, -10)
    ax[i].grid()
plt.tight_layout()
plt.show()
plot 2 13C glycine multi spectra fit

Create a fitting model

Spin System: The objective of a multi-dataset fitting is to optimize the spin system parameters using multiple datasets. In this example, we create two single-site spin systems, which are then shared by three method objects.

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

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

Method: Create the three MAS method objects with respective MAS spinning speeds.

# Get the spectral dimension parameters from the respective experiment and setup the
# corresponding method.

# Method for dataset 1
spectral_dims1 = get_spectral_dimensions(experiment1)
MAS1 = BlochDecaySpectrum(
    channels=["13C"],
    magnetic_flux_density=7.05,  # in T
    rotor_frequency=5000,  # in Hz
    spectral_dimensions=spectral_dims1,
    experiment=experiment1,  # add experimental dataset 1
)

# Method for dataset 2
spectral_dims2 = get_spectral_dimensions(experiment2)
MAS2 = BlochDecaySpectrum(
    channels=["13C"],
    magnetic_flux_density=7.05,  # in T
    rotor_frequency=1940,  # in Hz
    spectral_dimensions=spectral_dims2,
    experiment=experiment2,  # add experimental dataset 2
)

# Method for dataset 3
spectral_dims3 = get_spectral_dimensions(experiment3)
MAS3 = BlochDecaySpectrum(
    channels=["13C"],
    magnetic_flux_density=7.05,  # in T
    rotor_frequency=960,  # in Hz
    spectral_dimensions=spectral_dims3,
    experiment=experiment3,  # add experimental dataset 3
)

Guess Model Spectrum

# Simulation
# ----------
# Add the spin systems and the three methods to the simulator object.
sim = Simulator(spin_systems=spin_systems, methods=[MAS1, MAS2, MAS3])
sim.config.decompose_spectrum = "spin_system"
sim.run()

# Post Simulation Processing
# --------------------------
# Add signal processing to simulation dataset from the three methods.

# Processor for dataset 1
processor1 = 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=10),  # dataset is scaled independently using scale factor.
    ]
)

# Processor for dataset 2
processor2 = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Exponential(FWHM="30 Hz", dv_index=0),  # spin system 0
        sp.apodization.Exponential(FWHM="300 Hz", dv_index=1),  # spin system 1
        sp.FFT(),
        sp.Scale(factor=100),  # dataset is scaled independently using scale factor.
    ]
)

# Processor for dataset 3
processor3 = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Exponential(FWHM="10 Hz", dv_index=0),  # spin system 0
        sp.apodization.Exponential(FWHM="150 Hz", dv_index=1),  # spin system 1
        sp.FFT(),
        sp.Scale(factor=50),  # dataset is scaled independently using scale factor.
    ]
)
processors = [processor1, processor2, processor3]

processed_dataset = []
for i, proc in enumerate(processors):
    processed_dataset.append(
        proc.apply_operations(dataset=sim.methods[i].simulation).real
    )


# Plot of the guess Spectrum
# --------------------------

fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"})
for i, exp_dataset in enumerate(experiments):
    ax[i].plot(exp_dataset, color="black", linewidth=0.5, label="Experiment")
    ax[i].plot(processed_dataset[i], linewidth=2, alpha=0.6)
    ax[i].set_xlim(280, -10)
    ax[i].grid()

plt.legend()
plt.tight_layout()
plt.show()
plot 2 13C glycine multi spectra fit

Least-squares minimization with LMFIT

Use the make_LMFIT_params() for a quick setup of the fitting parameters. Note, the first two arguments of this function is the simulator object and a list of SignalProcessor objects, processors. The fitting parameters corresponding to the signal processor objects are generated using SP_i_operation_j_FunctionName_FunctionArg, where i is the ith signal processor within the list, j is the operation index of the ith processor, and FunctionName and FunctionArg are the operation function name and function argument, respectively.

params = sf.make_LMFIT_params(sim, processors, 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                10     -inf      inf     True     None
SP_1_operation_1_Exponential_FWHM            30     -inf      inf     True     None
SP_1_operation_2_Exponential_FWHM           300     -inf      inf     True     None
SP_1_operation_4_Scale_factor               100     -inf      inf     True     None
SP_2_operation_1_Exponential_FWHM            10     -inf      inf     True     None
SP_2_operation_2_Exponential_FWHM           150     -inf      inf     True     None
SP_2_operation_4_Scale_factor                50     -inf      inf     True     None
mth_0_rotor_frequency                      5000     4900     5100     True     None
mth_1_rotor_frequency                      1940     1840     2040     True     None
mth_2_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        60     -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, processors, sigmas))
result = minner.minimize()
result

Fit Statistics

fitting methodleastsq
# function evals182
# data points12288
# variables19
chi-square 15689.6963
reduced chi-square 1.27880808
Akaike info crit. 3040.95414
Bayesian info crit. 3181.86533

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift 176.116887 5.4044e-04 (0.00%) 176.0 -inf inf True
sys_0_site_0_shielding_symmetric_zeta 72.3189698 0.15407965 (0.21%) 60.0 -inf inf True
sys_0_site_0_shielding_symmetric_eta 0.91811712 0.00396939 (0.43%) 0.6 0.00000000 1.00000000 True
sys_0_abundance 56.9212657 0.13659678 (0.24%) 50.0 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift 43.3550919 0.00384271 (0.01%) 43.0 -inf inf True
sys_1_site_0_shielding_symmetric_zeta 22.7833633 0.20784235 (0.91%) 30.0 -inf inf True
sys_1_site_0_shielding_symmetric_eta 0.44484360 0.05491572 (12.34%) 0.5 0.00000000 1.00000000 True
sys_1_abundance 43.0787343 0.13659678 (0.32%) 50.0 0.00000000 100.000000 False 100-sys_0_abundance
mth_0_rotor_frequency 5033.17026 0.41764162 (0.01%) 5000.0 4900.00000 5100.00000 True
mth_1_rotor_frequency 1940.76008 0.05775086 (0.00%) 1940.0 1840.00000 2040.00000 True
mth_2_rotor_frequency 962.164634 0.04765175 (0.00%) 960.0 860.000000 1060.00000 True
SP_0_operation_1_Exponential_FWHM 27.2770184 0.19323483 (0.71%) 20.0 -inf inf True
SP_0_operation_2_Exponential_FWHM 109.089988 1.31094156 (1.20%) 200.0 -inf inf True
SP_0_operation_4_Scale_factor 38.6045004 0.17629942 (0.46%) 10.0 -inf inf True
SP_1_operation_1_Exponential_FWHM 29.3880214 0.16504812 (0.56%) 30.0 -inf inf True
SP_1_operation_2_Exponential_FWHM 136.816220 0.99234027 (0.73%) 300.0 -inf inf True
SP_1_operation_4_Scale_factor 170.686810 0.57431206 (0.34%) 100.0 -inf inf True
SP_2_operation_1_Exponential_FWHM 34.4272669 0.30018680 (0.87%) 10.0 -inf inf True
SP_2_operation_2_Exponential_FWHM 170.832502 1.84941368 (1.08%) 150.0 -inf inf True
SP_2_operation_4_Scale_factor 170.583254 0.91131011 (0.53%) 50.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_1_site_0_shielding_symmetric_zetasys_1_site_0_shielding_symmetric_eta-0.6387
SP_0_operation_1_Exponential_FWHMSP_0_operation_4_Scale_factor0.5932
SP_2_operation_1_Exponential_FWHMSP_2_operation_4_Scale_factor0.5890
SP_1_operation_1_Exponential_FWHMSP_1_operation_4_Scale_factor0.5196
SP_1_operation_2_Exponential_FWHMSP_1_operation_4_Scale_factor0.4924
SP_0_operation_2_Exponential_FWHMSP_0_operation_4_Scale_factor0.4602
sys_0_abundanceSP_1_operation_2_Exponential_FWHM-0.4373
sys_0_site_0_shielding_symmetric_zetasys_0_site_0_shielding_symmetric_eta-0.4317
SP_2_operation_2_Exponential_FWHMSP_2_operation_4_Scale_factor0.3709
sys_0_abundanceSP_0_operation_2_Exponential_FWHM-0.3607
sys_0_abundanceSP_1_operation_1_Exponential_FWHM0.3026
sys_0_abundanceSP_0_operation_4_Scale_factor-0.2782
sys_0_abundancesys_1_site_0_shielding_symmetric_zeta-0.2019
SP_0_operation_1_Exponential_FWHMSP_0_operation_2_Exponential_FWHM0.1930
sys_0_abundanceSP_2_operation_2_Exponential_FWHM-0.1880
sys_0_abundanceSP_1_operation_4_Scale_factor-0.1870
SP_2_operation_1_Exponential_FWHMSP_2_operation_2_Exponential_FWHM0.1854
sys_1_site_0_shielding_symmetric_zetaSP_2_operation_4_Scale_factor0.1681
sys_0_abundanceSP_2_operation_1_Exponential_FWHM0.1672
SP_0_operation_2_Exponential_FWHMSP_1_operation_2_Exponential_FWHM0.1660
sys_0_abundanceSP_2_operation_4_Scale_factor-0.1613
sys_0_site_0_shielding_symmetric_zetasys_0_abundance0.1549
sys_0_abundanceSP_0_operation_1_Exponential_FWHM0.1530
sys_0_site_0_shielding_symmetric_zetaSP_1_operation_4_Scale_factor0.1457
SP_0_operation_4_Scale_factorSP_1_operation_2_Exponential_FWHM0.1315
SP_1_operation_2_Exponential_FWHMSP_2_operation_2_Exponential_FWHM0.1306
sys_1_site_0_shielding_symmetric_zetaSP_1_operation_4_Scale_factor0.1293
mth_1_rotor_frequencySP_1_operation_1_Exponential_FWHM0.1196
SP_1_operation_1_Exponential_FWHMSP_1_operation_2_Exponential_FWHM0.1142
SP_0_operation_2_Exponential_FWHMSP_1_operation_1_Exponential_FWHM-0.1128
sys_0_abundancesys_1_site_0_shielding_symmetric_eta-0.1090
sys_1_site_0_shielding_symmetric_etaSP_2_operation_4_Scale_factor0.1074
sys_1_site_0_shielding_symmetric_zetaSP_2_operation_2_Exponential_FWHM-0.1031
SP_1_operation_2_Exponential_FWHMSP_2_operation_1_Exponential_FWHM-0.1018


The best fit solution

all_best_fit = sf.bestfit(sim, processors)  # a list of best fit simulations
all_residuals = sf.residuals(sim, processors)  # a list of residuals

# Plot the spectrum
fig, ax = plt.subplots(1, 3, figsize=(12, 3), subplot_kw={"projection": "csdm"})
for i, proc in enumerate(processors):
    ax[i].plot(experiments[i], color="black", linewidth=0.5, label="Experiment")
    ax[i].plot(all_residuals[i].real, color="gray", linewidth=0.5, label="Residual")
    ax[i].plot(all_best_fit[i].real, linewidth=2, alpha=0.6)
    ax[i].set_xlim(280, -10)
    ax[i].grid()

plt.legend()
plt.tight_layout()
plt.show()
plot 2 13C glycine multi spectra fit

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

Gallery generated by Sphinx-Gallery