¹³C 2D MAT NMR of L-Histidine

The following is an illustration for fitting 2D MAT/PASS datasets. The example dataset is a \(^{13}\text{C}\) 2D MAT spectrum of L-Histidine from Walder et al. [1]

import numpy as np
import csdmpy as cp
import matplotlib.pyplot as plt
from lmfit import Minimizer

from mrsimulator import Simulator
from mrsimulator.method.lib import SSB2D
from mrsimulator import signal_processor as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions
from mrsimulator.utils.collection import single_site_system_generator

Import the dataset

host = "https://ssnmr.org/sites/default/files/mrsimulator/"
filename = "1H13C_CPPASS_LHistidine.csdf"
mat_dataset = cp.load(host + filename)

# standard deviation of noise from the dataset
sigma = 0.4192854

# For the spectral fitting, we only focus on the real part of the complex dataset.
mat_dataset = mat_dataset.real

# Convert the coordinates along each dimension from Hz to ppm.
_ = [item.to("ppm", "nmr_frequency_ratio") for item in mat_dataset.dimensions]

When using the SSB2D method, ensure the horizontal dimension of the dataset is the isotropic dimension. Here, we apply an appropriate transpose operation to the dataset.

mat_dataset = mat_dataset.T  # transpose

# plot of the dataset.
max_amp = mat_dataset.max()
levels = (np.arange(24) + 1) * max_amp / 25  # contours are drawn at these levels.
options = dict(levels=levels, alpha=0.75, linewidths=0.5)  # plot options

plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_dataset, colors="k", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 LHistidine PASS

Create a fitting model

Guess model

Create a guess list of spin systems.

shifts = [120, 128, 135, 175, 55, 25]  # in ppm
zeta = [-70, -65, -60, -60, -10, -10]  # in ppm
eta = [0.8, 0.4, 0.9, 0.3, 0.0, 0.0]

spin_systems = single_site_system_generator(
    isotope="13C",
    isotropic_chemical_shift=shifts,
    shielding_symmetric={"zeta": zeta, "eta": eta},
    abundance=100 / 6,
)

Method

Create the SSB2D method.

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

PASS = SSB2D(
    channels=["13C"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=1500,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=mat_dataset,  # 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 = PASS.get_transition_pathways(sys)

Guess Spectrum

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

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        # Lorentzian convolution along the isotropic dimensions.
        sp.FFT(dim_index=0),
        sp.apodization.Exponential(FWHM="50 Hz"),
        sp.IFFT(dim_index=0),
        sp.Scale(factor=212260),
    ]
)
processed_dataset = processor.apply_operations(dataset=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_dataset, colors="k", **options)
ax.contour(processed_dataset, colors="r", linestyles="--", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 LHistidine PASS

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            50     -inf      inf     True     None
SP_0_operation_3_Scale_factor          2.123e+05     -inf      inf     True     None
sys_0_abundance                           16.67        0      100     True     None
sys_0_site_0_isotropic_chemical_shift       120     -inf      inf     True     None
sys_0_site_0_shielding_symmetric_eta        0.8        0        1     True     None
sys_0_site_0_shielding_symmetric_zeta       -70     -inf      inf     True     None
sys_1_abundance                           16.67        0      100     True     None
sys_1_site_0_isotropic_chemical_shift       128     -inf      inf     True     None
sys_1_site_0_shielding_symmetric_eta        0.4        0        1     True     None
sys_1_site_0_shielding_symmetric_zeta       -65     -inf      inf     True     None
sys_2_abundance                           16.67        0      100     True     None
sys_2_site_0_isotropic_chemical_shift       135     -inf      inf     True     None
sys_2_site_0_shielding_symmetric_eta        0.9        0        1     True     None
sys_2_site_0_shielding_symmetric_zeta       -60     -inf      inf     True     None
sys_3_abundance                           16.67        0      100     True     None
sys_3_site_0_isotropic_chemical_shift       175     -inf      inf     True     None
sys_3_site_0_shielding_symmetric_eta        0.3        0        1     True     None
sys_3_site_0_shielding_symmetric_zeta       -60     -inf      inf     True     None
sys_4_abundance                           16.67        0      100     True     None
sys_4_site_0_isotropic_chemical_shift        55     -inf      inf     True     None
sys_4_site_0_shielding_symmetric_eta          0        0        1     True     None
sys_4_site_0_shielding_symmetric_zeta       -10     -inf      inf     True     None
sys_5_abundance                           16.67        0      100    False 100-sys_0_abundance-sys_1_abundance-sys_2_abundance-sys_3_abundance-sys_4_abundance
sys_5_site_0_isotropic_chemical_shift        25     -inf      inf     True     None
sys_5_site_0_shielding_symmetric_eta          0        0        1     True     None
sys_5_site_0_shielding_symmetric_zeta       -10     -inf      inf     True     None
None

Solve the minimizer using LMFIT

Fit Statistics

fitting methodleastsq
# function evals184
# data points32768
# variables25
chi-square 54900.1891
reduced chi-square 1.67670003
Akaike info crit. 16960.3971
Bayesian info crit. 17170.3273

Variables

name value standard error relative error initial value min max vary expression
sys_0_site_0_isotropic_chemical_shift 119.171472 0.00371241 (0.00%) 120.0 -inf inf True
sys_0_site_0_shielding_symmetric_zeta -72.1497225 0.32674232 (0.45%) -70.0 -inf inf True
sys_0_site_0_shielding_symmetric_eta 0.98433203 0.00761883 (0.77%) 0.8 0.00000000 1.00000000 True
sys_0_abundance 16.2121326 0.07759578 (0.48%) 16.666666666666664 0.00000000 100.000000 True
sys_1_site_0_isotropic_chemical_shift 128.196863 0.00311768 (0.00%) 128.0 -inf inf True
sys_1_site_0_shielding_symmetric_zeta -75.5752331 0.27374721 (0.36%) -65.0 -inf inf True
sys_1_site_0_shielding_symmetric_eta 0.94600110 0.00581075 (0.61%) 0.4 0.00000000 1.00000000 True
sys_1_abundance 20.4506825 0.07808314 (0.38%) 16.666666666666664 0.00000000 100.000000 True
sys_2_site_0_isotropic_chemical_shift 136.197001 0.00472752 (0.00%) 135.0 -inf inf True
sys_2_site_0_shielding_symmetric_zeta -86.2975277 0.38515604 (0.45%) -60.0 -inf inf True
sys_2_site_0_shielding_symmetric_eta 0.42650009 0.00812494 (1.91%) 0.9 0.00000000 1.00000000 True
sys_2_abundance 12.3786320 0.07827283 (0.63%) 16.666666666666664 0.00000000 100.000000 True
sys_3_site_0_isotropic_chemical_shift 172.996047 0.00304270 (0.00%) 175.0 -inf inf True
sys_3_site_0_shielding_symmetric_zeta -68.6805894 0.24913973 (0.36%) -60.0 -inf inf True
sys_3_site_0_shielding_symmetric_eta 0.97842942 0.00623753 (0.64%) 0.3 0.00000000 1.00000000 True
sys_3_abundance 19.3442016 0.07563897 (0.39%) 16.666666666666664 0.00000000 100.000000 True
sys_4_site_0_isotropic_chemical_shift 54.5174353 0.00147239 (0.00%) 55.0 -inf inf True
sys_4_site_0_shielding_symmetric_zeta -20.1121436 0.12766278 (0.63%) -10.0 -inf inf True
sys_4_site_0_shielding_symmetric_eta 0.41788671 0.03675663 (8.80%) 0.0 0.00000000 1.00000000 True
sys_4_abundance 18.1057090 0.05440404 (0.30%) 16.666666666666664 0.00000000 100.000000 True
sys_5_site_0_isotropic_chemical_shift 26.9910802 0.00162727 (0.01%) 25.0 -inf inf True
sys_5_site_0_shielding_symmetric_zeta -10.3456036 0.55449432 (5.36%) -10.0 -inf inf True
sys_5_site_0_shielding_symmetric_eta 0.87675181 0.27990294 (31.92%) 0.0 0.00000000 1.00000000 True
sys_5_abundance 13.5086423 0.05063896 (0.37%) 16.666666666666693 0.00000000 100.000000 False 100-sys_0_abundance-sys_1_abundance-sys_2_abundance-sys_3_abundance-sys_4_abundance
SP_0_operation_1_Exponential_FWHM 98.5827516 0.31355107 (0.32%) 50.0 -inf inf True
SP_0_operation_3_Scale_factor 414385.328 935.779129 (0.23%) 212260.0 -inf inf True

Correlations (unreported correlations are < 0.100)

sys_5_site_0_shielding_symmetric_zetasys_5_site_0_shielding_symmetric_eta0.9427
sys_4_site_0_shielding_symmetric_zetasys_4_site_0_shielding_symmetric_eta0.7190
SP_0_operation_1_Exponential_FWHMSP_0_operation_3_Scale_factor0.5643
sys_3_site_0_shielding_symmetric_zetasys_3_site_0_shielding_symmetric_eta0.4389
sys_1_site_0_shielding_symmetric_zetasys_1_site_0_shielding_symmetric_eta0.4377
sys_0_site_0_shielding_symmetric_zetasys_0_site_0_shielding_symmetric_eta0.4299
sys_2_site_0_shielding_symmetric_zetasys_2_site_0_shielding_symmetric_eta0.3401
sys_4_site_0_shielding_symmetric_zetasys_4_abundance-0.2906
sys_0_site_0_shielding_symmetric_zetasys_0_abundance-0.2905
sys_3_site_0_shielding_symmetric_zetasys_3_abundance-0.2842
sys_4_abundanceSP_0_operation_3_Scale_factor-0.2750
sys_0_abundancesys_1_abundance-0.2744
sys_1_site_0_shielding_symmetric_zetasys_1_abundance-0.2704
sys_1_abundancesys_2_abundance-0.2697
sys_1_abundancesys_3_abundance-0.2618
sys_0_abundancesys_3_abundance-0.2565
sys_2_abundancesys_3_abundance-0.2469
sys_0_abundancesys_2_abundance-0.2324
sys_2_site_0_shielding_symmetric_etasys_2_abundance0.2230
sys_2_site_0_shielding_symmetric_zetasys_2_abundance-0.2182
sys_2_abundancesys_4_abundance-0.2073
sys_1_site_0_isotropic_chemical_shiftsys_1_abundance0.1994
sys_0_abundancesys_4_abundance-0.1842
sys_4_site_0_isotropic_chemical_shiftSP_0_operation_1_Exponential_FWHM-0.1617
sys_2_abundanceSP_0_operation_3_Scale_factor0.1607
sys_1_site_0_isotropic_chemical_shiftSP_0_operation_1_Exponential_FWHM-0.1606
sys_4_site_0_shielding_symmetric_etasys_4_abundance0.1584
sys_3_abundancesys_4_abundance-0.1542
sys_1_abundancesys_4_abundance-0.1530
sys_3_site_0_shielding_symmetric_zetaSP_0_operation_3_Scale_factor-0.1177
sys_0_site_0_shielding_symmetric_zetaSP_0_operation_3_Scale_factor-0.1163
sys_1_site_0_shielding_symmetric_zetaSP_0_operation_3_Scale_factor-0.1140


The best fit solution

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

# Plot of the best fit solution
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_dataset, colors="k", **options)
ax.contour(best_fit, colors="r", linestyles="--", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()
plot 1 LHistidine PASS

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

Gallery generated by Sphinx-Gallery