It’s !
どうもこんにちは。二回目の投稿です。今回はPythonでX線回折のシミュレーションをしてみました。コード全文(Jupyter Notebookを使っております)はGitHub¹
に載せています。間違っている点などがあれば、ご指摘いただけると幸いです。

X線回折(X-Ray Diffraction, XRD)

X線回折とは、X線が結晶中で回折する現象のことです。X線とは、波長が一定範囲内(1 nm ~ 1 pm程度)の電磁波のことで、回折とは、波の性質を持つものが、障害物を廻り込んで伝播することです。波は回折すると、互いに干渉することで強め合うところと相殺するところができ、縞状に分布するようになります。X線回折による分析では、回折によって生じた模様を解析することで対象の結晶構造を特定します。
ちなみに、回折現象はX線だけでなく、電子線や中性子線でも起き、それらを用いた結晶構造解析法も広く使われています。

X線回折による結晶構造解析の原理

X線回折によってできた模様では、回折による干渉によって強めあったところ(回折条件を満たしたところ)がスポットとして得られます。回折条件と結晶構造の関係をブラッグ条件といい、
$$2dsin{\theta} = n\lambda$$
で表されます。dは結晶中の原子がつくる面(結晶面)同士の間隔、$\theta$は結晶面と入射X線のなす角度、$\lambda$はX線の波長です。

原子散乱因子と結晶構造因子

先に説明したブラッグ条件では、X線を散乱する原子を点としていますが、実際にX線を散乱するのは、原子のまわりに分布している電子雲であり、散乱したX線の振幅は、(弾性散乱であるとすれば)散乱された位置での電子密度に比例します。そのため、原子によって散乱したX線の振幅fは、
$$f = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$$
となります。この式は、電子密度$\rho(\vec{r})$を散乱ベクトル$\vec{k}$でフーリエ変換することに相当します。散乱ベクトルとは、逆空間中に成分(h, k, l)で定義されるベクトルのことです。このfを原子散乱因子と言います。原子散乱因子を求めるには、電子密度が必要ですが、原子のまわりの電子密度は単純な場合を除いて解析的に求めることができないため、ハートリー-フォック法などで近似計算されます。各原子の原子散乱因子はすでに求められ、データベース化されているため、今回の計算ではそれらを用います。

さらに、結晶全体においても、同様の考え方を使うことができます。結晶中の電子密度をフーリエ変換することで結晶構造因子Fが得られます。
$$F(\vec{k}) = \int{\rho(\vec{r})e^{2π\vec{k}\vec{r}}} dV$$
式が原子散乱因子と一緒ですね。しかし、結晶中においては、原子位置の周辺にしか電子は分布していないので、これを原子の座標(x, y, z)と原子散乱因子fで書き換えることができ、
$$F(hkl) = \sum_{j}{f_jT_je^{2πi(hx_j + ky_j + lz_j)}}$$
となります。h, k, lは散乱ベクトルkの成分で、逆空間での格子点(逆格子点)の座標でもあります。
$T_j$は後ほど説明しますが、デバイ-ワラー因子というものです。

プログラム

ここからが実際の実装方法になります。
まず、逆格子点hklを列挙します。本来であれば、エワルト球上の点のみを挙げればいいのですが、X線の波長ではかなり幅広い範囲がエワルト球に載ってしまうので、なるだけ多く列挙します。GitHubにあるhkl.csvは列挙されたhklのファイルです。
結晶の単位格子の原子座標も必要になります。こちらは、単位格子内の原子だけで充分です。今回は、例として面心立方格子(fcc)の座標を用いています。GitHubではpos.csvという名前のファイルになっています。

次に、格子定数a, b, c, $\alpha$, $\beta$, $\gamma$から(これらはパラメータとしてあらかじめ持っておく)、実空間での格子テンソル$G$を組み立てます。
$$
G = \left(
\begin{matrix}
a^{2} & ab\cos\gamma & ac\cos\beta\\
ba\cos\gamma & b^{2} & bc\cos\alpha\\
ca\cos\beta & cb\cos\alpha & c^{2}
\end{matrix}
\right)
$$
次に、
$$G^{\ast} = G^{-1}$$
から、この格子テンソルの逆行列をとれば、逆空間での格子テンソル$G^{\ast}$が求まります。

G = [[a**2, a*b*np.cos(gamma), a*c*np.cos(beta)], [b*a*np.cos(gamma), b**2, b*c*np.cos(alpha)], [c*a*np.cos(beta), c*b*np.cos(alpha), c**2]]
invG = np.linalg.inv(G)

この逆格子テンソル$G^{\ast}$を使って、逆格子点hklでの反射Kに対する逆格子ベクトルの長さ
$$|\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}|$$
を求めます。
$$
|\vec{r}| = |ha^{\ast} + kb^{\ast} + lc^{\ast}| = h^{\mathrm{T}}G^{\ast}h$$

hkl = [h[i], k[i], l[i]]
thkl = np.transpose(hkl)
dk = 1/((np.matmul(np.matmul(invG, thkl), hkl))**(1/2))

|r|とブラッグ条件から、反射Kでの面間隔dkとブラッグ角θkが求まります。
$$d_k = \frac{1}{|\vec{r}|^2}$$
$$\theta_k = \arcsin(\frac{λ}{2d_k})$$

sinthetak = lamb/(2*dk)
thetak = np.rad2deg(np.arcsin(sinthetak))

原子散乱因子にはデータベースの値を使いますが、原子散乱因子はsinθ/λに依存する関数なので、データベースに載っているのはその係数だけです。係数から組み立てるには、
$$f(\frac{\sin\theta}{\lambda}) = \sum_{j=1}^{4}{a_j e^{-b_j(\frac{\sin\theta}{λ})^2}} + c$$
とします。$a_j$, $b_j$, cが係数です。プログラムでは、Na原子を想定してデータを用いました。

ここで、デバイ-ワラー因子について説明しなければなりません。
デバイ-ワラー因子とは、熱振動の影響を補正する因子です。等方的な熱振動のデバイ-ワラー因子は
$$T = e^{-B(\frac{\sin\theta}{\lambda})^2}$$
で表されます。Bは等方性原子変位パラメータというものです。
熱振動が等方的でない場合(異方性と言いますが)、より複雑で、$3\times3$対称行列$\beta$を用いて
$$T = e^{-h^{\mathrm{T}} β h}$$
となります。
今回は等方的な熱振動を考え、B = 2.0にしています。原子変位パラメータの値を探したのですが、見つからず、適当な値です。すいません。

原子散乱因子と結晶の単位格子の原子座標をもとに、結晶構造因子を計算します。
$$F_i(x, y, z, h_i, k_i, l_i, f_i) = f_i\sum_{j}{e^{2πi(hx_j + ky_j + lz_j)}}$$
見ての通り、結晶構造因子は複素数になります。

求めた結晶構造因子の絶対値の二乗(散乱振幅の絶対値の二乗なので、散乱強度になる)を縦軸、$2\theta$を横軸としてプロットすれば、X線回折で得られるようなプロットが得られます。入射X線の角度が$\theta$、散乱X線の角度が$\theta$なので、横軸を$2\theta$で取ることが習慣となっています。

結果

プロットすると以下のようになります。

ブラッグ角は小数点2桁までで丸めています。
このプロットでは、ブラッグ条件を満たす反射の角度でのみデルタ関数的に散乱強度が出ていますが、実際のX線回折では、散乱強度のピークは幅を持っています。そのピークの広がりを再現してみましょう。

デルタ関数に広がりを持たせたいときに、まず思いつくのは、ガウス関数(ガウシアン)による畳み込みでしょう。いわゆる、ガウシアンコンボリューションというやつです。
さっきのプロットに、ガウス関数を畳み込んだものが以下になります。分散は適当にそれっぽく見える値を入れています。

裾野に少し幅ができたのがわかると思います。X線回折の計測データを見たことがある人ならわかると思いますが、実際の物に比べて、このプロットは少し裾野が狭いように見えます。実はガウス関数は、粒子の熱振動などの挙動を表現するには適していますが、このようなスペクトルのピークの表現には適していません。
そこで、スペクトルのピークの表現に使われるローレンツ関数で畳み込みをしてみます。ローレンツ関数は以下のような関数です。
$$L = \frac{1}{1 + (\frac{x}{γ})^2}$$
$\gamma$は半値幅といい、ピーク高さの1/2でのピーク幅です。このパラメータで関数の広がり方が変わります。

def lorentzian(theta, gamma):
    lorentz = 1/(1 + (theta/gamma)**2)
    return lorentz

実際に畳み込んだものが以下です。またも半値幅の値はそれっぽく適当な値です。

先ほどのガウス関数よりも裾野が広くなりました。X線回折のピークっぽくなりましたね。
と言っても、ローレンツ関数だけでは正確ではありません。結晶内の原子が温度によって振動するために(たとえ絶対零度であっても零点振動があります)、ガウス関数的な挙動の熱振動もピークの広がりに影響を与えます。よってローレンツ関数とガウス関数を任意の重さで畳み込んだもので、より正確にピークを表現できます。

ローレンツ関数とガウス関数を畳み込んだものをフォークト関数と言います。フォークト関数は愚直にローレンツ関数とガウス関数を畳み込んで作るよりも、複素誤差関数を用いて作る方が早いです。筆者は勉強不足でよくわかっていませんが、なんか複素誤差関数の実部がフォークト関数っぽくなっているみたいです。²
ちょうどscipyに複素誤差関数を生成してくれる関数があるらしいので、今回はそれを使いました。

def voigt(theta, Hg, Hl):
    z = (theta + 1j*Hl)/(Hg * np.sqrt(2.0))
    w = scipy.special.wofz(z)
    v = (w.real)/(Hg * np.sqrt(2.0*np.pi))   
    return v

フォークト関数でピークの広がりを表現したものが以下です。

やってみたはいいものの、広がりをローレンツ関数のみで表現したものの方がそれっぽいですね。フォークト関数はパラメータとしてガウス関数部、ローレンツ関数部それぞれの半値幅が必要ですが、今回は適当に設定しています。半値幅を理論的に近づけることで、実際のピークの広がりに近づくかもしれません。

回折像

X線回折においては、上記のようなプロットの他に、回折像と言われる、X線の回折をそのまま写しとった像が得られます。回折像では逆格子が見て取れるので、今回はおまけとして作ってみようと思います。

実装としては、任意の大きさでnumpy.zeros()の配列を生成しておいて、反射の条件を満たすところのみ|F|^2の値を入れ、スポットっぽく見えるようにガウス関数を畳み込みます。

def gaussian(theta, lamb):
    gauss = (1/(2*np.pi*sigma))*np.exp(-(theta**2 + lamb**2)/(2*sigma))
    return gauss

回折像は入射X線の方向に対して垂直な逆格子点のみがスポットとして出てきます。よって、今までの条件に加えて、入射X線方向についても考慮しなければなりません。入射X線方向を[001]として、内積をとって0になる反射のみを取り出します。この処理は[001]だと楽ですが、入射方向が[110]や[111]などでは複雑な処理が必要でしょう。
Pythonでは、配列をそのままmatplotlib.pyplot.imshow()で画像にできますので、白黒のコントラストで画像にします。

t, l = np.meshgrid(np.arange(-1, 1, 0.1), np.arange(-1, 1, 0.1))
g = gaussian(t, l)
xrdimage = scipy.signal.convolve2d(I, g, boundary='wrap', mode='same')
plt.imshow(xrdimage, interpolation="nearest", cmap=plt.cm.gray)

便利ですね。fccでの像は下のようになりました。

それっぽいですね。さっきからそれっぽいしか言っていませんね。体心立方格子(bcc)でもやってみました。

fccより楽しい感じになりましたね。

まとめ

X線回折のシミュレーションをしてみました。回折理論が複雑ですが、コンピュータでする処理自体は簡単ですので、X線回折を理解するにはちょうどいい教材かもしれません。
実を言うと、精密なシミュレーションをするには、今回実装したものよりも多くの因子を考えなければなりません(ローレンツ・偏光因子や選択配向関数など)。加えて、今回適当に設定したパラメータ(原子変位パラメータ、フォークト関数の半値幅)も本来は正確な値を使うべきです。では、これらのパラメータは理論的に求めることができるのでしょうか。これらのパラメータは様々な要因が複雑に影響し合ってできているため、支配的な要因から類推はできても、解析的に求めるのは難しいでしょう。
上記のようなパラメータの正確な値を得るためには、実験値とシミュレーションをすり合わせ、誤差を減らしていく、という手法を用います。そこで、前回の記事³
で紹介した、非線形最小二乗法が活躍するわけです。

参考文献

今回の記事は、以下の記事や書籍を元に作成しています。勝手ながら紹介させていただきます。

• 粉末X線解析の実際　リートベルト法入門　(日本分析化学会X線分析研究懇談会　編　中井泉　泉富士夫　編著)

• pythonでフォークト関数(Voigt function)をプロットする方法

コード

• https://github.com/PicricAcid/XRD