実装して学ぶデータ処理システムその2 - 反射が起きる条件

スポット検出の話が間に入ってしまったが、第 1 回に引き続き、Kabsch, Wolfgang. "Integration, scaling, space-group assignment and post-refinement" Acta Crystallographica Section D (2010): 133-144 の 2.2 節で述べられている反射条件を説明する。

以下では、1つの反射 hkl に着目するので、添字に指数を明示しない。回転前( $\varphi = 0$ )の結晶方位における逆格子点は、逆格子の基底ベクトルを用いて $\mathbf{p_0^*} = h\mathbf{b_0^*} + k\mathbf{b_1^*} + l\mathbf{b_2^*}$ と表現できる。これを $\mathbf{m_2}$ 周りに $\varphi$ 回転させて $\mathbf{p^*}$ に移動し、回折が起きるようにしたい。

このときの様子を、回転軸 $\mathbf{m_2}$ に垂直な視線で眺めたのが下の図である。

結晶の回転により、 $\mathbf{p_0^*}$ は点線の軌道を移動して、 $\mathbf{p^*}$ において Ewald 球に接触する。これを、回転軸に平行な視線から眺める(上の図で、上空から見下ろす)と、 $\mathbf{p_0^*}$ の軌跡は $\mathbf{m_1}, \mathbf{m_3}$ 平面に平行な平面上の円になる。Ewald 球のこの高さでの断面は、 $|\mathbf{S_0}|$ 以下の半径を持つことに注意。

回折が起きるとき、Laue 方程式 $|\mathbf{S_0}| = |\mathbf{S}| = |\mathbf{S_0 + p^*|$ が成り立つ。両辺を二乗して、 $|\mathbf{S_0}|^2 = |\mathbf{S_0}|^2 + |\mathbf{p^*}|^2 + 2 \mathbf{S_0}\cdot\mathbf{p^*}$ だから、 $|\mathbf{p_0^*}|^2 + 2 \mathbf{S_0}\cdot\mathbf{p^*} = 0$ である。ここで、 $\mathbf{S_0}$ と $\mathbf{p^*}$ を goniostat system の基底を使って展開してやると、
$|\mathbf{p_0^*}|^2 + 2 (S_0^1\mathbf{m_1} + S_0^2\mathbf{m_2} + S_0^2\mathbf{m_3})(p_0^1\mathbf{m_1} + p_0^2\mathbf{m_2} + p_0^2\mathbf{m_3}) = 0$
となる。太字になっていない係数は $S_0^1 = \mathbf{S_0}\cdot\mathbf{m_1}$ のような成分であって、スカラーである。
$\mathbf{m_1}, \mathbf{m_2}, \mathbf{m_3}$ は正規直交系をなすので、
$|\mathbf{p^*}|^2 + 2 (S_0^1 p^1 + S_0^2 p^2 + S_0^3 p^3) = 0$
と簡単にできる。さらに、 $\mathbf{m_1} = \mathbf{m_2} \times \mathbf{S_0}$ という定義から $\mathbf{S_0}$ と $\mathbf{m_1}$ は直交するから、 $S_0^1 = 0$ である。また、回転は $\mathbf{m_2}$ 軸のまわりに起きるので、この軸方向の成分は変化せず、
$p_0^2 = p^2$
である。以上より、
$|\mathbf{p^*}|^2 + 2 (S_0^2 p^{2} + S_0^3 p^3) = 0$
となって、これを移項して整理すると、
$p^3 = -\frac{\frac{|\mathbf{p_0}|^2}{2} - S_0^2 p^{2}}{S_0^3}$
が得られる。最後の成分 $p^1$ は、回転によって $\mathbf{p}$ の長さが変わらないことから、
$p^1 = \pm\sqrt{|\mathbf{p_0^*}|^2 - (p^2)^2 - (p^3)^2}$
と求められる。

次に、回転角 $\varphi$ を求める。 $\mathbf{m_1}, \mathbf{m_3}$ 平面内に射影して考えればよい。
$\mathbf{p_0^*}$ を $\mathbf{p^*}$ に重ねるための角度だから、内積の定義から
$\cos\varphi = \frac{(p_0^{1}\mathbf{m_1}+p_0^{3}\mathbf{m_3})\cdot(p^{1}\mathbf{m_1}+p^{3}\mathbf{m_3})}{|p_0^{*1}\mathbf{m_1}+p_0^{*3}\mathbf{m_3}||p^{1}\mathbf{m_1}+p^3\mathbf{m_3}|} = \frac{p_0^{1}p^{1}+p_0^{3}p^{3}}{|\mathbf{p_0^*}|^2 - (p_0^2)^2}$
が分かる。sin については、(二次元の)外積の性質から
$\sin\varphi = \frac{(p_0^{1}\mathbf{m_1}+p_0^{3}\mathbf{m_3})\times(p^{1}\mathbf{m_1}+p^{3}\mathbf{m_3})}{|p_0^{1}\mathbf{m_1}+p_0^{3}\mathbf{m_3}||p^{1}\mathbf{m_1}+p^*3\mathbf{m_3}|} = \frac{p_0^{3}p^{1} - p_0^{1}p^{3}}{|\mathbf{p_0^*}|^2 - (p_0^2)^2}$
となる。

この sin と cos から回転角 $\varphi$ を求めることもできるが、実際の計算では三角関数のまま用いる。その理由は、φ の回折パラメータによる偏微分 - biochem_fan's note に書いた。

今回の計算は R に実装済みであるが、とりあえず理論だけ公開しておく。図を作るのが思ったより面倒で、コードの整理が間に合わなかった。反射が不可能な条件 (blind region) と合わせて、次回とする。