²³Na MAS NMR of Nasicon

The following is a least-squares fitting example of a \(^{23}\text{Na}\) MAS NMR spectrum of Nasicon, \(\text{NaZr}_2(\text{PO}_4)_3\). 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, report_fit

from mrsimulator import Simulator, Site, SpinSystem
from mrsimulator.methods import BlochDecayCTSpectrum
from mrsimulator import signal_processing as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions

Import the dataset

host = "https://nmr.cemhti.cnrs-orleans.fr/Dmfit/Help/csdm/"
filename = "23Na QUAD MAS Nasicon.csdf"
experiment = cp.load(host + filename)

# standard deviation of noise from the dataset
sigma = 0.2368151

# 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(100, -150)
plt.grid()
plt.tight_layout()
plt.show()

Create a fitting model

Spin System

Na23 = Site(
    isotope="23Na",
    isotropic_chemical_shift=-20.0,  # in ppm
    quadrupolar={"Cq": 2.3e6, "eta": 0.03},  # Cq in Hz
)
spin_systems = [SpinSystem(sites=[Na23])]

Method

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

MAS_CT = BlochDecayCTSpectrum(
    channels=["23Na"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=15000,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=experiment,  # add the measurement to the method.
)

# 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_CT.get_transition_pathways(sys)

Guess Model Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[MAS_CT])
sim.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        sp.IFFT(),
        sp.apodization.Exponential(FWHM="100 Hz"),
        sp.FFT(),
        sp.Scale(factor=200),
    ]
)
processed_data = processor.apply_operations(data=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_data, linewidth=2, alpha=0.6, label="Guess Spectrum")
ax.set_xlim(100, -150)
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)
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Exponential_FWHM           100     -inf      inf     True     None
SP_0_operation_3_Scale_factor               200     -inf      inf     True     None
sys_0_abundance                             100        0      100    False      100
sys_0_site_0_isotropic_chemical_shift       -20     -inf      inf     True     None
sys_0_site_0_quadrupolar_Cq             2.3e+06     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta               0.03        0        1     True     None
None

Solve the minimizer using LMFIT

minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma))
result = minner.minimize()
report_fit(result)

Out:

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 56
    # data points      = 4096
    # variables        = 5
    chi-square         = 50387.1068
    reduced chi-square = 12.3165746
    Akaike info crit   = 10289.8313
    Bayesian info crit = 10321.4201
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift: -21.4706356 +/- 0.00298717 (0.01%) (init = -20)
    sys_0_site_0_quadrupolar_Cq:            2252000.71 +/- 205.098216 (0.01%) (init = 2300000)
    sys_0_site_0_quadrupolar_eta:           0.00341673 +/- 6.9412e-04 (20.32%) (init = 0.03)
    sys_0_abundance:                        100.000000 +/- 0.00000000 (0.00%) == '100'
    SP_0_operation_1_Exponential_FWHM:      45.7664986 +/- 0.38445046 (0.84%) (init = 100)
    SP_0_operation_3_Scale_factor:          213.721815 +/- 0.23474326 (0.11%) (init = 200)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_eta)      = -0.877
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_Cq)       =  0.876
    C(sys_0_site_0_quadrupolar_eta, SP_0_operation_1_Exponential_FWHM)          = -0.722
    C(sys_0_site_0_quadrupolar_Cq, sys_0_site_0_quadrupolar_eta)                = -0.680
    C(sys_0_site_0_isotropic_chemical_shift, SP_0_operation_1_Exponential_FWHM) =  0.505
    C(SP_0_operation_1_Exponential_FWHM, SP_0_operation_3_Scale_factor)         =  0.470
    C(sys_0_site_0_quadrupolar_Cq, SP_0_operation_1_Exponential_FWHM)           =  0.384
    C(sys_0_site_0_quadrupolar_eta, SP_0_operation_3_Scale_factor)              = -0.240
    C(sys_0_site_0_isotropic_chemical_shift, SP_0_operation_3_Scale_factor)     =  0.187
    C(sys_0_site_0_quadrupolar_Cq, SP_0_operation_3_Scale_factor)               =  0.168

The best fit solution

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

# 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(100, -150)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()

1

D.Massiot, F.Fayon, M.Capron, I.King, S.Le Calvé, B.Alonso, J.O.Durand, B.Bujoli, Z.Gan, G.Hoatson, ‘Modelling one and two-dimensional solid-state NMR spectra.’, Magn. Reson. Chem. 40 70-76 (2002) DOI: 10.1002/mrc.984