CoCl₂.2D₂O, ²H (I=1) Shifting-d echo

²H (I=1) 2D NMR CSA-Quad 1st order correlation spectrum.

The following is an example of fitting static shifting-d echo NMR correlation spectrum of \(\text{NiCl}_2\cdot 2\text{D}_2\text{O}\) crystalline solid. The spectrum used here is from Walder et al. 1.

import numpy as np
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 Method2D
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

filename = "https://sandbox.zenodo.org/record/830903/files/CoCl2.2D2O.csdf"
experiment = cp.load(filename)

# standard deviation of noise from the dataset
sigma = 11.11578

# 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.
max_amp = experiment.max()
levels = (np.arange(24) + 1) * max_amp / 25  # contours are drawn at these levels.
options = dict(levels=levels, linewidths=0.5)  # plot options

plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.contour(experiment, colors="k", **options)
ax.set_xlim(2500, -1500)
ax.set_ylim(2000, -2200)
plt.grid()
plt.tight_layout()
plt.show()

Create a fitting model

Guess model

Create a guess list of spin systems.

site = Site(
    isotope="2H",
    isotropic_chemical_shift=200,  # in ppm
    shielding_symmetric={
        "zeta": -1300,  # in ppm
        "eta": 0.2,
        "alpha": np.pi,  # in rads
        "beta": np.pi / 2,  # in rads
        "gamma": np.pi / 2,  # in rads
    },
    quadrupolar={"Cq": 110e3, "eta": 0.83},  # Cq in Hz
)

spin_systems = [SpinSystem(sites=[site])]

Method

Use the generic 2D method, Method2D, to generate a shifting-d echo method.

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

shifting_d = Method2D(
    channels=["2H"],
    magnetic_flux_density=9.395,  # in T
    spectral_dimensions=[
        {
            **spectral_dims[0],
            "label": "Quadrupolar frequency",
            "events": [
                {
                    "rotor_frequency": 0,
                    "transition_query": {"P": [-1]},
                    "freq_contrib": ["Quad1_2"],
                }
            ],
        },
        {
            **spectral_dims[1],
            "label": "Paramagnetic shift",
            "events": [
                {
                    "rotor_frequency": 0,
                    "transition_query": {"P": [-1]},
                    "freq_contrib": ["Shielding1_0", "Shielding1_2"],
                }
            ],
        },
    ],
    experiment=experiment,  # also 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 = shifting_d.get_transition_pathways(sys)

Guess Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[shifting_d])
sim.config.integration_volume = "hemisphere"
sim.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        # Gaussian convolution along both dimensions.
        sp.IFFT(dim_index=0),
        sp.apodization.Gaussian(FWHM="10 kHz", dim_index=0),  # along dimension 0
        sp.FFT(dim_index=0),
        sp.Scale(factor=5e8),
    ]
)
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.contour(experiment, colors="k", **options)
ax.contour(processed_data, colors="r", linestyles="--", **options)
ax.set_xlim(2500, -1500)
ax.set_ylim(2000, -2200)
plt.grid()
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["sys_0_site_0_shielding_symmetric_alpha"].vary = False
params["sys_0_site_0_shielding_symmetric_beta"].vary = False
params["sys_0_site_0_shielding_symmetric_gamma"].vary = False
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                       Value      Min      Max     Vary     Expr
SP_0_operation_1_Gaussian_FWHM                10     -inf      inf     True     None
SP_0_operation_3_Scale_factor              5e+08     -inf      inf     True     None
sys_0_abundance                              100        0      100    False      100
sys_0_site_0_isotropic_chemical_shift        200     -inf      inf     True     None
sys_0_site_0_quadrupolar_Cq              1.1e+05     -inf      inf     True     None
sys_0_site_0_quadrupolar_eta                0.83        0        1     True     None
sys_0_site_0_shielding_symmetric_alpha     3.142     -inf      inf    False     None
sys_0_site_0_shielding_symmetric_beta      1.571     -inf      inf    False     None
sys_0_site_0_shielding_symmetric_eta         0.2        0        1     True     None
sys_0_site_0_shielding_symmetric_gamma     1.571     -inf      inf    False     None
sys_0_site_0_shielding_symmetric_zeta      -1300     -inf      inf     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   = 442
    # data points      = 65536
    # variables        = 7
    chi-square         = 309489.005
    reduced chi-square = 4.72293191
    Akaike info crit   = 101747.037
    Bayesian info crit = 101810.670
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift:   202.512829 +/- 0.89419451 (0.44%) (init = 200)
    sys_0_site_0_shielding_symmetric_zeta:  -1355.90491 +/- 1.48289502 (0.11%) (init = -1300)
    sys_0_site_0_shielding_symmetric_eta:    0.22077495 +/- 0.00550582 (2.49%) (init = 0.2)
    sys_0_site_0_shielding_symmetric_alpha:  3.141593 (fixed)
    sys_0_site_0_shielding_symmetric_beta:   1.570796 (fixed)
    sys_0_site_0_shielding_symmetric_gamma:  1.570796 (fixed)
    sys_0_site_0_quadrupolar_Cq:             114245.459 +/- 94.9713991 (0.08%) (init = 110000)
    sys_0_site_0_quadrupolar_eta:            0.99999667 +/- 0.00311737 (0.31%) (init = 0.83)
    sys_0_abundance:                         100.000000 +/- 0.00000000 (0.00%) == '100'
    SP_0_operation_1_Gaussian_FWHM:          37.8199216 +/- 0.18725437 (0.50%) (init = 10)
    SP_0_operation_3_Scale_factor:           9.3760e+08 +/- 1661029.24 (0.18%) (init = 5e+08)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_0_site_0_quadrupolar_Cq, sys_0_site_0_quadrupolar_eta)                    = -0.935
    C(sys_0_site_0_shielding_symmetric_eta, SP_0_operation_1_Gaussian_FWHM)         = -0.636
    C(SP_0_operation_1_Gaussian_FWHM, SP_0_operation_3_Scale_factor)                =  0.355
    C(sys_0_site_0_shielding_symmetric_eta, sys_0_site_0_quadrupolar_eta)           = -0.354
    C(sys_0_site_0_shielding_symmetric_eta, sys_0_site_0_quadrupolar_Cq)            =  0.347
    C(sys_0_site_0_quadrupolar_eta, SP_0_operation_1_Gaussian_FWHM)                 =  0.202
    C(sys_0_site_0_quadrupolar_Cq, SP_0_operation_1_Gaussian_FWHM)                  = -0.197
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_shielding_symmetric_zeta) =  0.153
    C(sys_0_site_0_shielding_symmetric_zeta, SP_0_operation_3_Scale_factor)         = -0.132
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_site_0_quadrupolar_Cq)           = -0.117
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_site_0_quadrupolar_eta)          =  0.115
    C(sys_0_site_0_isotropic_chemical_shift, sys_0_site_0_quadrupolar_Cq)           = -0.112
    C(sys_0_site_0_shielding_symmetric_zeta, SP_0_operation_1_Gaussian_FWHM)        =  0.106

The best fit solution

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

# Plot the spectrum
plt.figure(figsize=(4.25, 3.0))
ax = plt.subplot(projection="csdm")
ax.contour(experiment, colors="k", **options)
ax.contour(best_fit, colors="r", linestyles="--", **options)
ax.set_xlim(2500, -1500)
ax.set_ylim(2000, -2200)
plt.grid()
plt.tight_layout()
plt.show()

Image plots with residuals

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

fig, ax = plt.subplots(
    1, 3, sharey=True, figsize=(10, 3.0), subplot_kw={"projection": "csdm"}
)
vmax, vmin = experiment.max(), experiment.min()
for i, dat in enumerate([experiment, best_fit, residuals]):
    ax[i].imshow(dat, aspect="auto", cmap="gist_ncar_r", vmax=vmax, vmin=vmin)
    ax[i].set_xlim(2500, -1500)
ax[0].set_ylim(2000, -2200)
plt.tight_layout()
plt.show()

1

Walder B.J, Patterson A.M., Baltisberger J.H, and Grandinetti P.J Hydrogen motional disorder in crystalline iron group chloride dihydrates spectroscopy, J. Chem. Phys. (2018) 149, 084503. DOI: 10.1063/1.5037151