四極子核のMAS NMRスペクトル

研究で無機の固体NMRを測定することが多いのですが、四極子結合定数がMHzくらいある原子核では、MASで平均化されない四極子相互作用の2次の項の影響で特徴的なピークの線形となります。+/-1/2の遷移について、MASの回転速度が十分に早ければ周波数領域でこの線形を比較的容易にシュミレーションできます。線形は非対称性因子に対応して特徴的な鬼の角のようなピークを示します。

某所との研究で必要になったので、飯島さんに助言してもらいつつPythonでスペクトルのシミュレーションプログラムを作りました。11Bのケミカルシフトが0 ppmで四極子結合定数が2.5 MHz（ほぼ3配位Bの値）のピークが磁場強度でどのように線形が変わるか計算した例です。パウダーの計算にalphaとbetaをゴージャスにサンプリングして、スペクトルを作りました。CSの分布を1, 10, 20, 50, 100Hzで計算しています。四極子結合定数を共鳴周波数で割り算するので、磁場が強いとどんどん先鋭化しますが、全体的に高磁場シフトしてピークが分裂します。800 MHzくらいの磁石になると2次の四極相互作用に由来するピークの幅は数ppmになるので、ほぼ鬼の角は観測されないでしょう。計算のために作成したPythonコードも参考に示します。

#!/usr/bin/env python
# Author: T. Ohkubo (Chiba-u) 2020/08/29
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('classic')
plt_dic = {}
plt_dic['font.size'] = 10
plt_dic['axes.labelsize'] = 10
plt_dic['figure.figsize'] = [16, 10]
plt.rcParams.update(plt_dic)

p = {}
p["gamma"] = 85.84e6  # rad/T/s (11B)
p["B0"] = [200, 400, 600, 800]  # 1H freq in MHz
p["I"] = 3/2
p["Cq"] = 2.5e6  # Hz
p["eta"] = [0, 0.25, 0.5, 0.75, 1.0]
p["bf"] = 1  # Hz
p["nprime"] = 18   # debug:12 fine:18
p["xaxis"] = (3, -35, 400)   # left right steps
p["outfile"] = "Qmas12-bf{:03d}.png".format(p["bf"])


def func(alpha, beta, eta):
    a = np.cos(2*alpha)
    b = np.cos(beta)
    A = 7/2 - 7/3*eta*a + 7/18*eta**2*a**2
    B = -3 + 2/9*eta**2 + 8/3*eta*a - 7/9*eta**2*a**2
    C = 5/6 - 1/3*eta*a + 7/18*eta**2*a**2
    f = A*b**4 + B*b**2 + C
    return f


def CalcSpectrum(I, Cq, nu0, f):
    nuq = Cq * 3 / (2*I*(2*I+1) - 3/4)
    nu = -1/16 * nuq**2/nu0**2 * (I*(I+1) - 3/4)*f
    return nu * 1e6


def SphericalZWC(m):
    phi = (1 + np.sqrt(5))/2
    M = (phi**m - (-phi)**(-m))/np.sqrt(5)
    N = (phi**(m+2) - (-phi)**(-(m+2)))/np.sqrt(5)
    M, N = int(M), int(N)
    j = np.arange(N, dtype=int)
    alpha = 2*np.pi * (j*M % N)/N
    beta = np.arccos(1-(2*j+1)/N)
    alpha_ = np.tile(alpha, beta.shape[0])
    beta_ = np.repeat(beta, alpha.shape[0])
    return alpha_, beta_


def WindowFunc(x, y, bf):
    # bf in ppm
    gpnt = 100
    dx = x[1] - x[0]
    xg = np.arange(-gpnt, gpnt + 1) * dx
    yg = np.exp(-np.log(2) * (xg/bf)**2)
    v = np.outer(y, yg)
    yf = np.zeros(y.shape[0] + 2 * gpnt)
    for i in range(v.shape[0]):
        yf[i:i+yg.shape[0]] = yf[i:i+yg.shape[0]] + v[i, :]
    yf = yf[gpnt:-gpnt]
    return x, yf


alpha, beta = SphericalZWC(p["nprime"])
etadata = {}
for eta in p["eta"]:
    print("Calculating eta={}...".format(eta))
    etadata[eta] = func(alpha, beta, eta)

fig, ax = plt.subplots(ncols=len(p["eta"]), nrows=len(p["B0"]))
bins = np.linspace(p["xaxis"][1], p["xaxis"][0], p["xaxis"][2])
for i, B0 in enumerate(p["B0"]):
    nu0 = B0/267.513 * p["gamma"]
    bf = p["bf"] / nu0 * 1e6  # ppm
    for j, eta in enumerate(p["eta"]):
        s = CalcSpectrum(p["I"], p["Cq"], nu0, etadata[eta])
        hist, bins = np.histogram(s, bins=bins)
        x = bins[0:-1] + (bins[1] - bins[0]) * 0.5
        x, hist = WindowFunc(x, hist, bf)
        ax[i][j].plot(x, hist)
        label = "$\\eta={}$\n${:.0f}\ {{\\rm MHz}}$".format(eta, B0)
        ax[i][j].text(0.75, 0.75, label, transform=ax[i][j].transAxes,
                      bbox=dict(facecolor='gray', alpha=0.1))
        ax[i][j].axes.yaxis.set_visible(False)
        ax[i][j].set_xlim(p["xaxis"][0:2])
        if i == len(p["B0"]) - 1:
            ax[i][j].set_xlabel("ppm")
fig.tight_layout()
fig.subplots_adjust(wspace=0.1, hspace=0.2)
fig.savefig(p["outfile"])
plt.show()