²H MAS NMR of Methionine

The following is a least-squares fitting example of a \(^{2}\text{H}\) MAS NMR spectrum of Methionine. The 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
from mrsimulator.spin_system.tensors import SymmetricTensor

Import the dataset

host = "https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/"
filename = "2H methiodine MAS.csdf"
experiment = cp.load(host + filename)

# standard deviation of noise from the dataset
sigma = 0.3026282

# 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=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.set_xlim(600, -700)
plt.grid()
plt.tight_layout()
plt.show()
plot 4 2H quad

Create a fitting model

Spin System

H_2 = Site(
    isotope="2H",
    isotropic_chemical_shift=-57.12,  # in ppm,
    quadrupolar=SymmetricTensor(Cq=3e4, eta=0.0),  # Cq in Hz
)

spin_systems = [SpinSystem(sites=[H_2])]

Method

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

MAS = BlochDecaySpectrum(
    channels=["2H"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=4517.1,  # 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.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Exponential(FWHM="60 Hz"),
        sp.FFT(),
        sp.Scale(factor=140),
    ]
)
processed_dataset = processor.apply_operations(dataset=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.plot(processed_dataset, linewidth=2, alpha=0.6, label="Guess Spectrum")
ax.set_xlim(600, -700)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plot 4 2H quad

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

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Exponential_FWHM            60     -inf      inf     True     None
SP_0_operation_3_Scale_factor               140     -inf      inf     True     None
sys_0_abundance                             100        0      100    False      100
sys_0_site_0_isotropic_chemical_shift    -57.12     -inf      inf    False     None
sys_0_site_0_quadrupolar_Cq               3e+04     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta                  0        0        1     True     None
None

Solve the minimizer using LMFIT

Fit Statistics

fitting methodleastsq
# function evals21
# data points8192
# variables4
chi-square 1350041.14
reduced chi-square 164.880452
Akaike info crit. 41825.9668
Bayesian info crit. 41854.0105

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift -57.1200000 0.00000000 (0.00%) -57.12 -inf inf False
sys_0_site_0_quadrupolar_Cq 33652.7963 79.9919768 (0.24%) 30000.0 -inf inf True
sys_0_site_0_quadrupolar_eta 0.32118873 0.00436216 (1.36%) 0.0 0.00000000 1.00000000 True
sys_0_abundance 100.000000 0.00000000 (0.00%) 100 0.00000000 100.000000 False 100
SP_0_operation_1_Exponential_FWHM 62.9584689 0.32076792 (0.51%) 60.0 -inf inf True
SP_0_operation_3_Scale_factor 132.896559 0.50526319 (0.38%) 140.0 -inf inf True

Correlations (unreported correlations are < 0.100)

SP_0_operation_1_Exponential_FWHMSP_0_operation_3_Scale_factor0.6665
sys_0_site_0_quadrupolar_Cqsys_0_site_0_quadrupolar_eta-0.2857
sys_0_site_0_quadrupolar_CqSP_0_operation_3_Scale_factor0.2601
sys_0_site_0_quadrupolar_etaSP_0_operation_3_Scale_factor0.1075


The best fit solution

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

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
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, label="Best Fit")
ax.set_xlim(600, -700)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plot 4 2H quad

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

Gallery generated by Sphinx-Gallery