前編から引き続いて分子の回転についてです。今回は分子を重ね合わせるための数式を導出したいと思います。
目的関数
四元数のイントロが終わってようやく本題に入ります。分子の3次元的な重ね合わせですが、動かす分子の座標ベクトルを\( \mathbf{r}_{\alpha}=(x_{\alpha}, y_{\alpha}, z_{\alpha}) \)として、重ねる先の参照分子の座標ベクトルを\( \mathbf{r}_{0\alpha}=(x_{0\alpha}, y_{0\alpha}, z_{0\alpha}) \)とします。ここで\( \alpha \)は原子のインデックスとします。分子を動かして参照分子に重ねるということは、2つの分子の座標の差を最小にするという問題になります。
まず分子座標の残差ベクトルは以下のようになります。
$$
\mathbf{\epsilon}_{\alpha} = R\left(\mathbf{r}_{\alpha} – \mathbf{g}\right) – \left(\mathbf{r}_{0\alpha}-\mathbf{g}_{0}\right)
$$
ここで\(R\)は回転行列で、\(\mathbf{g}\)と\(\mathbf{g}_{0}\)はそれぞれ動かす分子と参照分子の母核重心（中心）ベクトルとなります。 2つの分子の母核重心を原点に移動させることは簡単なことです。予め動かす分子と参照分子の母核重心は座標系の原点に一致しているものとすれば、 以後は純粋に回転のみを考えればよいことになります。 したがって以後の議論では動かす分子と参照分子の母核重心は原点にあるものとして
$$
\mathbf{g}=(0,0,0) \\
\mathbf{g}_{0}=(0,0,0)
$$
であることを仮定します。そのうえで上記の残差を四元数を使って表現すると、
$$
\epsilon_{\alpha} = q r_{\alpha} \bar{q} – r_{0\alpha}
$$
となります。この残差を母核を構成する原子に全てにおいて最小化したいので、以下のような関数を最小化することになります。
$$
m(q) = \displaystyle\sum_{\alpha} \bar{\epsilon_{\alpha}} \epsilon_{\alpha}
$$
この関数を最小化するような四元数\(q\)を求めれば良いのです。ただし、回転を表す四元数は大きさが1なので、\(q^2=\bar{q}q=1\)という条件も課すことになります。
最小化する四元数を求める
上記の残差四元数のままだと計算が面倒なので、少し変形します。\(\epsilon_{\alpha}\)の右から\(q\)を掛けたものを考えます。
$$
\epsilon_{\alpha}^{\prime}=\epsilon_{\alpha}q = q r_{\alpha} – r_{0\alpha}q
$$
そして目的関数は次のような形になります。
$$
\begin{align}
m^{\prime}(q) &= \displaystyle\sum_{\alpha} \bar{\epsilon_{\alpha}^{\prime}} \epsilon_{\alpha}^{\prime} \\
&= \overline{(q r_{\alpha} – r_{0\alpha}q)}(q r_{\alpha} – r_{0\alpha}q)
\end{align}
$$
ここで\(q^2=\bar{q}q=1\)という条件を使うと\(m^{\prime}(q) = m(q)\)となります。ぜひ手を動かして確認してみてください。
四元数\(q=a+bi+cj+dk\)を4次元ベクトル\(\mathbf{\tilde{q}}=(a,b,c,d)\)として、とある\(4\times4\)の行列\(K_{\alpha}\)を使って上記目的関数は、
$$
m^{\prime}(q) = \displaystyle\sum_{\alpha}
\mathbf{\tilde{q}}^{\mathrm{T}}K_{\alpha}^{\mathrm{T}}K_{\alpha}\mathbf{\tilde{q}}
= \displaystyle\sum_{\alpha} \mathbf{\tilde{q}}^{\mathrm{T}}M_{\alpha}\mathbf{\tilde{q}}
$$
と書き表すことができます。行列\(M_{\alpha}=K_{\alpha}^{\mathrm{T}}K_{\alpha}\)は対称行列となります。そして四元数\(q\)の大きさが1であるという条件は、
$$
\mathbf{\tilde{q}}^{\mathrm{T}}\mathbf{\tilde{q}} = 1
$$
となります。したがって、この条件の下に\(m^{\prime}(q)\)を最小化するということは、
$$
L=\displaystyle\sum_{\alpha} \mathbf{\tilde{q}}^{\mathrm{T}}M_{\alpha}\mathbf{\tilde{q}}
-\lambda \left( \mathbf{\tilde{q}}^{\mathrm{T}}\mathbf{\tilde{q}} – 1 \right)
$$
を最小化するということになります。ここで\(\lambda\)はラグランジュの未定乗数です。最終的に以下のような固有値問題に帰着します。
$$
\displaystyle\sum_{\alpha} M_{\alpha}\mathbf{\tilde{q}}
= M\mathbf{\tilde{q}} = \lambda \mathbf{\tilde{q}} \\
M=\displaystyle\sum_{\alpha} M_{\alpha}
$$
4次の正方対称行列\(M\)を対角化して固有値最小の固有ベクトル\(\mathbf{\tilde{q}}=(a, b, c, d)\)を求めれば、四元数\(q=a+bi+cj+dk\)から回転行列\(R\)が計算できるので、参照分子の母核構造にピッタリ重なるように分子を回転できるはずです。次は\(M_{\alpha}\)の具体的な表式を計算してみましょう。
対角化すべき行列を求める
\(M_{\alpha}\)を計算するために、先ずは以下の式を計算して\(K_{\alpha}\)みましょう。
$$
\epsilon_{\alpha}^{\prime} = q r_{\alpha} – r_{0\alpha}q
= K_{\alpha} \mathbf{\tilde{q}}
$$
四元数\(q, r_{\alpha}, r_{0\alpha}\)のベクトル表現はそれぞれ以下の通りです。
$$
\begin{align}
&q = (a, \mathbf{q}) = (a, (b, c, d)) \\
&r_{\alpha} = (0, \mathbf{r}_{\alpha}) = (0, (x_{\alpha}, y_{\alpha}, z_{\alpha})) \\
&r_{0\alpha} = (0, \mathbf{r}_{0\alpha}) = (0, (x_{0\alpha}, y_{0\alpha}, z_{0\alpha}))
\end{align}
$$
前編で紹介したベクトル表現での積の公式を使うと、
$$
\begin{align}
&q r_{\alpha} = (-\mathbf{q} \cdot \mathbf{r}_{\alpha}, a\mathbf{r}_{\alpha} + \mathbf{q} \times \mathbf{r}_{\alpha}) \\ &r_{0\alpha} q = (-\mathbf{r}_{0\alpha} \cdot \mathbf{q}, a\mathbf{r}_{0\alpha} + \mathbf{r}_{0\alpha} \times \mathbf{q})
\end{align}
$$
2つの差をとって、
$$
\begin{align}
q r_{\alpha} – r_{0\alpha} q
&= ((\mathbf{r}_{0\alpha}-\mathbf{r}_{\alpha}) \cdot \mathbf{q}, \quad
-a(\mathbf{r}_{0\alpha}-\mathbf{r}_{\alpha}) + \mathbf{q} \times (\mathbf{r}_{0\alpha}+\mathbf{r}_{\alpha}) )\\
&=\begin{pmatrix}
(x_{0\alpha}-x_{\alpha})b+(y_{0\alpha}-y_{\alpha})c+(z_{0\alpha}-z_{\alpha})d \\
-(x_{0\alpha}-x_{\alpha})a+(z_{0\alpha}+z_{\alpha})c-(y_{0\alpha}+y_{\alpha})d \\
-(y_{0\alpha}-y_{\alpha})a-(z_{0\alpha}+z_{\alpha})b+(x_{0\alpha}+x_{\alpha})d \\
-(z_{0\alpha}-z_{\alpha})a+(y_{0\alpha}+y_{\alpha})b-(x_{0\alpha}+x_{\alpha})c
\end{pmatrix} \\
&=\begin{pmatrix}
0 & x_{0\alpha}-x_{\alpha} & y_{0\alpha}-y_{\alpha} & z_{0\alpha}-z_{\alpha} \\
-(x_{0\alpha}-x_{\alpha}) & 0 & z_{0\alpha}+z_{\alpha} & -(y_{0\alpha}+y_{\alpha}) \\
-(y_{0\alpha}-y_{\alpha}) & -(z_{0\alpha}+z_{\alpha}) & 0 & x_{0\alpha}+x_{\alpha} \\
-(z_{0\alpha}-z_{\alpha}) & y_{0\alpha}+y_{\alpha} & -(x_{0\alpha}+x_{\alpha}) & 0 \\
\end{pmatrix}
\begin{pmatrix} a \\ b \\ c \\ d \end{pmatrix} \\
&= K_{\alpha} \mathbf{\tilde{q}}
\end{align}
$$
これで行列\(K_{\alpha}\)が求まりました。
いよいよ行列\(M_{\alpha}\)を計算します。
この計算は見かけほど面倒ではないので根性でやりましょう。
$$
\begin{align}
M_{\alpha} &= K_{\alpha}^{\mathrm{T}} K_{\alpha} \\
&=\begin{pmatrix}
0 & -(x_{0\alpha}-x_{\alpha}) & -(y_{0\alpha}-y_{\alpha}) & -(z_{0\alpha}-z_{\alpha}) \\
x_{0\alpha}-x_{\alpha} & 0 & -(z_{0\alpha}+z_{\alpha}) & y_{0\alpha}+y_{\alpha} \\
y_{0\alpha}-y_{\alpha} & z_{0\alpha}+z_{\alpha} & 0 & -(x_{0\alpha}+x_{\alpha}) \\
z_{0\alpha}-z_{\alpha} & -(y_{0\alpha}+y_{\alpha}) & x_{0\alpha}+x_{\alpha} & 0 \\
\end{pmatrix}
\begin{pmatrix}
0 & x_{0\alpha}-x_{\alpha} & y_{0\alpha}-y_{\alpha} & z_{0\alpha}-z_{\alpha} \\
-(x_{0\alpha}-x_{\alpha}) & 0 & z_{0\alpha}+z_{\alpha} & -(y_{0\alpha}+y_{\alpha}) \\
-(y_{0\alpha}-y_{\alpha}) & -(z_{0\alpha}+z_{\alpha}) & 0 & x_{0\alpha}+x_{\alpha} \\
-(z_{0\alpha}-z_{\alpha}) & y_{0\alpha}+y_{\alpha} & -(x_{0\alpha}+x_{\alpha}) & 0 \\
\end{pmatrix} \\
&=\begin{pmatrix}
m_{\alpha 11} & m_{\alpha 12} & m_{\alpha 13} & m_{\alpha 14} \\
m_{\alpha 21} & m_{\alpha 22} & m_{\alpha 23} & m_{\alpha 14} \\
m_{\alpha 31} & m_{\alpha 32} & m_{\alpha 33} & m_{\alpha 14} \\
m_{\alpha 41} & m_{\alpha 42} & m_{\alpha 43} & m_{\alpha 44} \\
\end{pmatrix}
\end{align}
$$
式が長くて行列の形で書けないので、行列成分毎に書き下します。
また、\(M_{\alpha}\)は対称行列なので下半分だけの記述とします。
$$
\begin{align}
m_{\alpha 11} &= x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2
-2x_{\alpha}x_{0\alpha}-2y_{\alpha}y_{0\alpha}-2z_{\alpha}z_{0\alpha} \\
m_{\alpha 21} &= 2(z_{\alpha}y_{0\alpha}-y_{\alpha}z_{0\alpha}) \\
m_{\alpha 22} &= x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2
-2x_{\alpha}x_{0\alpha}+2y_{\alpha}y_{0\alpha}+2z_{\alpha}z_{0\alpha} \\
m_{\alpha 31} &= 2(x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) \\
m_{\alpha 32} &= -2(x_{\alpha}y_{0\alpha}+y_{\alpha}x_{0\alpha}) \\
m_{\alpha 33} &= x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2
+2x_{\alpha}x_{0\alpha}-2y_{\alpha}y_{0\alpha}+2z_{\alpha}z_{0\alpha} \\
m_{\alpha 41} &= 2(y_{\alpha}x_{0\alpha}-x_{\alpha}y_{0\alpha}) \\
m_{\alpha 42} &= -2(x_{\alpha}z_{0\alpha}+z_{\alpha}x_{0\alpha}) \\
m_{\alpha 43} &= -2(y_{\alpha}z_{0\alpha}+z_{\alpha}y_{0\alpha}) \\
m_{\alpha 44} &= x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2
+2x_{\alpha}x_{0\alpha}+2y_{\alpha}y_{0\alpha}-2z_{\alpha}z_{0\alpha} \\
\end{align}
$$
行列\(M\)は重ね合わせに使う原子に渡って\(M_{\alpha}\)の和をとれば良いので、
$$
\scriptsize
\begin{align}
M &= \displaystyle\sum_{\alpha} M_{\alpha} \\
&=\begin{pmatrix}
\displaystyle\sum_{\alpha} \left( d_{\alpha}^2-2(x_{\alpha}x_{0\alpha}+y_{\alpha}y_{0\alpha}+z_{\alpha}z_{0\alpha}) \right)&
\displaystyle\sum_{\alpha} 2(z_{\alpha}y_{0\alpha}-y_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} 2(x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} 2(y_{\alpha}x_{0\alpha}-x_{\alpha}y_{0\alpha}) \\
\displaystyle\sum_{\alpha} 2(z_{\alpha}y_{0\alpha}-y_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} \left( d_{\alpha}^2+2(-x_{\alpha}x_{0\alpha}+y_{\alpha}y_{0\alpha}+z_{\alpha}z_{0\alpha}) \right)&
\displaystyle\sum_{\alpha} -2(x_{\alpha}y_{0\alpha}+y_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} -2(x_{\alpha}z_{0\alpha}+z_{\alpha}x_{0\alpha}) \\
\displaystyle\sum_{\alpha} 2(x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} -2(x_{\alpha}y_{0\alpha}+y_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} \left( d_{\alpha}^2+2(x_{\alpha}x_{0\alpha}-y_{\alpha}y_{0\alpha}+z_{\alpha}z_{0\alpha}) \right)&
\displaystyle\sum_{\alpha} -2(y_{\alpha}z_{0\alpha}+z_{\alpha}y_{0\alpha}) \\
\displaystyle\sum_{\alpha} 2(y_{\alpha}x_{0\alpha}-x_{\alpha}y_{0\alpha}) &
\displaystyle\sum_{\alpha} -2(x_{\alpha}z_{0\alpha}+z_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} -2(y_{\alpha}z_{0\alpha}+z_{\alpha}y_{0\alpha}) &
\displaystyle\sum_{\alpha} \left( d_{\alpha}^2+2(x_{\alpha}x_{0\alpha}+y_{\alpha}y_{0\alpha}-z_{\alpha}z_{0\alpha}) \right) \\
\end{pmatrix}
\\
d_{\alpha}^2&=x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2
\end{align}
$$
となります。
そして数値計算へ
分子の母核同士を重ね合わせる回転行列\(R\)を求めるには、上記の\(4 \times 4\)行列\(M\)を対角化して、最小固有値の固有ベクトルから四元数を構成すれば良いことは分かりました。行列\(M\)をそのままプログラミングして計算しても良いのですが、もう少しシンプルな形に変形したいと思います。行列をよりシンプルな形に変形できればプログラミングするときのミスが減るだろうし、数値誤差も軽減できる場合もあります。
行列\(M\)を以下のように変形します。
$$
M = 2M^{\prime}+\displaystyle\sum_{\alpha}(x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2)I
$$
ここで\(I\)は\(4 \times 4\)の単位行列です。そして\(M^{\prime}\)は
$$
M^{\prime} = \begin{pmatrix}
\displaystyle\sum_{\alpha} (-x_{\alpha}x_{0\alpha}-y_{\alpha}y_{0\alpha}-z_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} (z_{\alpha}y_{0\alpha}-y_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} (x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha})&
\displaystyle\sum_{\alpha} (y_{\alpha}x_{0\alpha}-x_{\alpha}y_{0\alpha}) \\
\displaystyle\sum_{\alpha} (z_{\alpha}y_{0\alpha}-y_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} (-x_{\alpha}x_{0\alpha}+y_{\alpha}y_{0\alpha}+z_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} (-x_{\alpha}y_{0\alpha}-y_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} (-x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) \\
\displaystyle\sum_{\alpha} (x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} (-x_{\alpha}y_{0\alpha}-y_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} (x_{\alpha}x_{0\alpha}-y_{\alpha}y_{0\alpha}+z_{\alpha}z_{0\alpha}) &
\displaystyle\sum_{\alpha} (-y_{\alpha}z_{0\alpha}-z_{\alpha}y_{0\alpha}) \\
\displaystyle\sum_{\alpha} (y_{\alpha}x_{0\alpha}-x_{\alpha}y_{0\alpha}) &
\displaystyle\sum_{\alpha} (-x_{\alpha}z_{0\alpha}-z_{\alpha}x_{0\alpha}) &
\displaystyle\sum_{\alpha} (-y_{\alpha}z_{0\alpha}-z_{\alpha}y_{0\alpha}) &
\displaystyle\sum_{\alpha} (x_{\alpha}x_{0\alpha}+y_{\alpha}y_{0\alpha}-z_{\alpha}z_{0\alpha}) \\
\end{pmatrix}
$$
となります。この行列\(M^{\prime}\)の固有ベクトル\(\mathbf{u}\)は行列\(M\)の固有ベクトルでもあります。
$$
\begin{align}
M \mathbf{u} &= 2M^{\prime} \mathbf{u} + \displaystyle\sum_{\alpha}(x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2)I \mathbf{u} \\
&= \left( 2\lambda^{\prime} + \displaystyle\sum_{\alpha}(x_{\alpha}^2+y_{\alpha}^2+z_{\alpha}^2+x_{0\alpha}^2+y_{0\alpha}^2+z_{0\alpha}^2) \right) \mathbf{u} \\
&= \lambda \mathbf{u}
\end{align}
$$
したがってプログラミングするときは行列\(M\)よりもシンプルな\(M^{\prime}\)を用いることにします。
次回は実際にC++プログラムを作ってみます。
Category: OpenBabel, プログラミング関連