²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 numpy as np
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)

# 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()

Estimate noise statistics from the dataset

coords = experiment.dimensions[0].coordinates
noise_region = np.where(np.logical_and(coords > -250e-6, coords < -210e-6))
noise_data = experiment[noise_region]

plt.figure(figsize=(3.75, 2.5))
ax = plt.subplot(projection="csdm")
ax.plot(noise_data, label="noise")
plt.title("Noise section")
plt.axis("off")
plt.tight_layout()
plt.show()

noise_mean, sigma = experiment[noise_region].mean(), experiment[noise_region].std()
noise_mean, sigma

(<Quantity -0.1204676>, <Quantity 0.35856736>)

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
)

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=1400),
    ]
)
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()

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"]))

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Exponential_FWHM            60     -inf      inf     True     None
SP_0_operation_3_Scale_factor              1400     -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

opt = sim.optimize()  # Pre-compute transition pathways
minner = Minimizer(
    sf.LMFIT_min_function,
    params,
    fcn_args=(sim, processor, sigma),
    fcn_kws={"opt": opt},
)
result = minner.minimize()
result

Fit Result

Fit Statistics
fitting method	leastsq
# function evals	20
# data points	8192
# variables	4
chi-square	985575.891
reduced chi-square	120.368331
Akaike info crit.	39248.2376
Bayesian info crit.	39276.2812

Parameters
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	33250.8411	78.8622873	(0.24%)	30000.0	-inf	inf	True
sys_0_site_0_quadrupolar_eta	0.26590115	4.88125965	(1835.74%)	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	63.0830703	0.32289068	(0.51%)	60.0	-inf	inf	True
SP_0_operation_3_Scale_factor	1242.71354	4.78433047	(0.38%)	1400.0	-inf	inf	True

Correlations (unreported values are < 0.100)
Parameter1	Parameter 2	Correlation
SP_0_operation_1_Exponential_FWHM	SP_0_operation_3_Scale_factor	+0.6607
sys_0_site_0_quadrupolar_Cq	sys_0_site_0_quadrupolar_eta	+0.2384
sys_0_site_0_quadrupolar_Cq	SP_0_operation_3_Scale_factor	+0.2360
sys_0_site_0_quadrupolar_eta	SP_0_operation_3_Scale_factor	-0.1860

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()

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

Gallery generated by Sphinx-Gallery