分子構造の重ねあわせ (行列、最小二乗法など？)
すいません、化学っぽい内容の質問なのですが、
アルゴリズムは数学になると思うのでこちらで質問させていただきます。
また、説明不足を補うために一部プログラムっぽい表記をさせていただきますのでご了承ください。
現在、分子構造を比較するためのプログラムを作成しています。
そのプログラムに構造の「重ね合わせ」機能を追加したいのですが
そのアルゴリズムに最小二乗法などが必要らしく、悩んでいます。
(詳細は長くなるので下に書きます)
分かる方がいましたらご教示お願いします。
よろしくお願いします。
■重ね合わせについて
分子構造のデータは、各原子がそれぞれX, Y, Z座標の3つの値を持っています。
例えば、200原子からなる分子のデータ構造は以下のようになります。
Atom { double x; double y; double z; }
Atom atom = new Atom[200];
重ね合わせる2つの構造は原子数は同じとします。なので、分子構造2つをa1, a2とし、
対応する点の誤差をとっていくと、誤差の平均を求める計算は以下のようになります。
Atom a1 = new Atom[200];
Atom a2 = new Atom[200];
double d = 0;
for (i = 0; i < a1.length; i++)
d += sqrt(pow((a1[i].x-a2[i].x),2) + pow((a1[i].y-a2[i].y),2) + pow((a1[i].z-a2[i].z),2));
d = d / a1.length;
このdの値が最小値をとるようにa2の原子の座標を変更することを「重ね合わせ」と定義します。
■重ね合わせのアルゴリズムについて
基準となる構造に対し、もう一つの構造を並進や回転によって合わせる感じになると思います。
今回質問する箇所はこのアルゴリズムになります。
この際に行列計算や最小二乗法を使うことになるかと思われます。
■質問者の考え
途中まで、しかも間違っている可能性がありますが私の考えたアルゴリズムを追記しておきます。
(i) a1, a2の重心を求める
double x1 = y1 = z1 = x2 = y2 = z2 = 0;
for (i = 0; i < a1.length; i++) {
x1 += a1[i].x; y1 += a1[i].y; z1 += a1[i].z;
x2 += a2[i].x; y2 += a2[i].y; z2 += a2[i].z;
}
x1 = x1 / a1.length; y1 = y1 / a1.length; z1 = z1 / a1.length;
x2 = x2 / a1.length; y2 = y2 / a1.length; z2 = z2 / a1.length;
a1の重心(x1, y1, z1)とa2の重心(x2, y2, z2)が求まります。
(ii) 重心の誤差を求め、誤差をa2の各原子の座標に反映させる
xd = (x2 - x1); yd = (y2 - y1); zd = (z2 - z1);
for (i = 0; i < a1.length; i++) {
a2[i].x -= xd;
a2[i].y -= yd;
a2[i].z -= zd;
}
この方法で誤差はかなり小さくなります。しかし、最小値にはなってないと思われます。
なので、この方法が正しいのかよく分かりません。
