²⁷Al MAS NMR of YAG (1st and 2nd order Quad)

The following is a quadrupolar lineshape fitting example for the 27Al MAS NMR of Yttrium aluminum garnet (YAG) crystal. 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, Site, SpinSystem
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 = "27Al Quad MAS YAG 400MHz.csdf"
experiment = cp.load(host + filename)

# standard deviation of noise from the dataset
sigma = 0.4895381

# 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=(8, 4))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.set_xlim(1200, -1200)
plt.grid()
plt.tight_layout()
plt.show()
plot 4 27Al YAG

Create a fitting model

Guess model

Create a guess list of spin systems.

Al_1 = Site(
    isotope="27Al",
    isotropic_chemical_shift=76,  # in ppm
    quadrupolar=SymmetricTensor(Cq=6e6, eta=0.0),  # Cq in Hz
)

Al_2 = Site(
    isotope="27Al",
    isotropic_chemical_shift=1,  # in ppm
    quadrupolar=SymmetricTensor(Cq=5e5, eta=0.3),  # Cq in Hz
)
spin_systems = [
    SpinSystem(sites=[Al_1], name="AlO4"),
    SpinSystem(sites=[Al_2], name="AlO6"),
]

Method

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

MAS = BlochDecaySpectrum(
    channels=["27Al"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=15250,  # 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.get_transition_pathways(sys)

Guess Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[MAS])
sim.config.decompose_spectrum = "spin_system"
sim.run()

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

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 4))
ax = plt.subplot(projection="csdm")
ax.plot(experiment, color="black", linewidth=0.5, label="Experiment")
ax.plot(processed_dataset, linewidth=2, alpha=0.6)
ax.set_xlim(1200, -1200)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plot 4 27Al YAG

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

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM              300     -inf      inf     True     None
SP_0_operation_3_Scale_factor                50     -inf      inf     True     None
mth_0_rotor_frequency                  1.525e+04 1.515e+04 1.535e+04     True     None
sys_0_abundance                              50        0      100     True     None
sys_0_site_0_isotropic_chemical_shift        76     -inf      inf     True     None
sys_0_site_0_quadrupolar_Cq               6e+06     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta                  0        0        1     True     None
sys_1_abundance                              50        0      100    False 100-sys_0_abundance
sys_1_site_0_isotropic_chemical_shift         1     -inf      inf     True     None
sys_1_site_0_quadrupolar_Cq               5e+05     -inf      inf     True     None
sys_1_site_0_quadrupolar_eta                0.3        0        1     True     None
None

Solve the minimizer using LMFIT

Fit Statistics

fitting methodleastsq
# function evals122
# data points4096
# variables10
chi-square 53126.2675
reduced chi-square 13.0020234
Akaike info crit. 10516.6578
Bayesian info crit. 10579.8355

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift 76.8355161 0.09433328 (0.12%) 76.0 -inf inf True
sys_0_site_0_quadrupolar_Cq 6023010.64 12064.0031 (0.20%) 6000000.0 -inf inf True
sys_0_site_0_quadrupolar_eta 0.04512660 0.00903697 (20.03%) 0.0 0.00000000 1.00000000 True
sys_0_abundance 47.8306865 0.45106200 (0.94%) 50.0 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift 0.84621252 0.00328035 (0.39%) 1.0 -inf inf True
sys_1_site_0_quadrupolar_Cq 587099.315 3851.74516 (0.66%) 500000.0 -inf inf True
sys_1_site_0_quadrupolar_eta 0.76658596 0.01589616 (2.07%) 0.3 0.00000000 1.00000000 True
sys_1_abundance 52.1693135 0.45106200 (0.86%) 50.0 0.00000000 100.000000 False 100-sys_0_abundance
mth_0_rotor_frequency 15251.4470 0.18079327 (0.00%) 15250.0 15150.0000 15350.0000 True
SP_0_operation_1_Gaussian_FWHM 319.998118 2.04120416 (0.64%) 300.0 -inf inf True
SP_0_operation_3_Scale_factor 142.027137 1.22767010 (0.86%) 50.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_0_abundanceSP_0_operation_3_Scale_factor0.8172
sys_1_site_0_quadrupolar_Cqsys_1_site_0_quadrupolar_eta-0.6812
sys_0_site_0_isotropic_chemical_shiftsys_0_site_0_quadrupolar_Cq0.6016
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_Cq0.5502
sys_1_site_0_isotropic_chemical_shiftmth_0_rotor_frequency-0.4901
sys_1_site_0_quadrupolar_Cqmth_0_rotor_frequency-0.3477
sys_1_site_0_isotropic_chemical_shiftsys_1_site_0_quadrupolar_eta-0.2255
sys_1_site_0_quadrupolar_etamth_0_rotor_frequency0.1754
sys_0_site_0_quadrupolar_etasys_0_abundance0.1709
sys_1_site_0_quadrupolar_CqSP_0_operation_1_Gaussian_FWHM-0.1640
sys_0_abundancesys_1_site_0_quadrupolar_Cq-0.1549
SP_0_operation_1_Gaussian_FWHMSP_0_operation_3_Scale_factor0.1487
sys_0_site_0_quadrupolar_etaSP_0_operation_3_Scale_factor0.1444
sys_0_abundancesys_1_site_0_quadrupolar_eta0.1328
sys_1_site_0_quadrupolar_CqSP_0_operation_3_Scale_factor0.1280
sys_1_site_0_quadrupolar_etaSP_0_operation_3_Scale_factor-0.1180
sys_0_abundanceSP_0_operation_1_Gaussian_FWHM-0.1096
sys_0_site_0_quadrupolar_CqSP_0_operation_3_Scale_factor0.1086
sys_0_site_0_quadrupolar_Cqsys_0_abundance0.1054


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=(8, 4))
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)
ax.set_xlim(1200, -1200)
plt.grid()
plt.legend()
plt.tight_layout()
plt.show()
plot 4 27Al YAG

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

Gallery generated by Sphinx-Gallery