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

# standard deviation of noise from the dataset
sigma = 3.258224

# 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()
plot 1 31p Na2PO4 static

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
)

# Optimize the script by pre-setting the transition pathways for each spin system from
# the method.
for sys in spin_systems:
    sys.transition_pathways = static1D.get_transition_pathways(sys)

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=4000),
    ]
)
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()
plot 1 31p Na2PO4 static

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

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM             3000     -inf      inf     True     None
SP_0_operation_3_Scale_factor              4000     -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

Fit Statistics

fitting methodleastsq
# function evals69
# data points2048
# variables5
chi-square 3405.48278
reduced chi-square 1.66690298
Akaike info crit. 1051.45512
Bayesian info crit. 1079.57821

Variables

name value standard error relative error initial value min max vary
sys_0_site_0_isotropic_chemical_shift 6.32228199 0.07291435 (1.15%) 5.0 -inf inf True
sys_0_site_0_shielding_symmetric_zeta -86.8324173 0.14007057 (0.16%) -80.0 -inf inf True
sys_0_site_0_shielding_symmetric_eta 0.18101165 0.02313119 (12.78%) 0.5 0.00000000 1.00000000 True
SP_0_operation_1_Gaussian_FWHM 5978.91425 101.488701 (1.70%) 3000.0 -inf inf True
SP_0_operation_3_Scale_factor 4626.51661 6.23103856 (0.13%) 4000.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_0_site_0_shielding_symmetric_etaSP_0_operation_1_Gaussian_FWHM-0.9747
sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_shielding_symmetric_zeta-0.5905
sys_0_site_0_shielding_symmetric_zetaSP_0_operation_1_Gaussian_FWHM0.3249
sys_0_site_0_shielding_symmetric_zetaSP_0_operation_3_Scale_factor-0.3075
sys_0_site_0_shielding_symmetric_zetasys_0_site_0_shielding_symmetric_eta-0.2963
sys_0_site_0_isotropic_chemical_shiftSP_0_operation_3_Scale_factor0.2785
sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_shielding_symmetric_eta-0.2489
sys_0_site_0_isotropic_chemical_shiftSP_0_operation_1_Gaussian_FWHM0.2161
SP_0_operation_1_Gaussian_FWHMSP_0_operation_3_Scale_factor0.1824
sys_0_site_0_shielding_symmetric_etaSP_0_operation_3_Scale_factor-0.1179


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()
plot 1 31p Na2PO4 static

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

Gallery generated by Sphinx-Gallery