Note
Go to the end to download the full example code.
³¹P static NMR of crystalline Na2PO4 (CSA)¶
The following is a CSA static least-squares fitting example of a \(^{31}\text{P}\) MAS NMR spectrum of \(\text{Na}_{2}\text{PO}_{4}\). 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 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 = "31P Phophonate Static.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(200, -200)
plt.grid()
plt.tight_layout()
plt.show()
Estimate noise statistics from the dataset
coords = experiment.dimensions[0].coordinates
noise_region = np.where(coords < -110e-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 -1.9085335>, <Quantity 3.3310158>)
Create a fitting model¶
Spin System
P_31 = Site(
isotope="31P",
isotropic_chemical_shift=5.0, # in ppm,
shielding_symmetric=SymmetricTensor(zeta=-80, eta=0.5), # zeta in Hz
)
spin_systems = [SpinSystem(sites=[P_31])]
Method
# Get the spectral dimension parameters from the experiment.
spectral_dims = get_spectral_dimensions(experiment)
static1D = BlochDecaySpectrum(
channels=["31P"],
magnetic_flux_density=9.395, # in T
rotor_frequency=0, # in Hz
rotor_angle=0, # in rads
spectral_dimensions=spectral_dims,
experiment=experiment, # experimental dataset
)
Guess Model Spectrum
# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[static1D])
sim.run()
# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
operations=[
sp.IFFT(),
sp.apodization.Gaussian(FWHM="3000 Hz"),
sp.FFT(),
sp.Scale(factor=40000),
]
)
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(200, -200)
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.pop("sys_0_abundance")
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))
Name Value Min Max Vary Expr
SP_0_operation_1_Gaussian_FWHM 3000 -inf inf True None
SP_0_operation_3_Scale_factor 4e+04 -inf inf True None
sys_0_site_0_isotropic_chemical_shift 5 -inf inf True None
sys_0_site_0_shielding_symmetric_eta 0.5 0 1 True None
sys_0_site_0_shielding_symmetric_zeta -80 -inf inf 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
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(200, -200)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
Total running time of the script: (0 minutes 1.862 seconds)