| 1 |
#include "math/RMSD.hpp" |
| 2 |
#include "math/SVD.hpp" |
| 3 |
|
| 4 |
using namespace oopse; |
| 5 |
using namespace JAMA; |
| 6 |
|
| 7 |
RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov, |
| 8 |
Vector3d mov_com, |
| 9 |
Vector3d mov_to_ref, |
| 10 |
RotMat3x3d U) { |
| 11 |
|
| 12 |
assert(mov.size() == ref_.size()); |
| 13 |
int n; |
| 14 |
int n_vec = ref_.size(); |
| 15 |
|
| 16 |
/* calculate the centre of mass */ |
| 17 |
mov_com = V3Zero; |
| 18 |
|
| 19 |
for (n=0; n < n_vec; n++) { |
| 20 |
mov_com += mov[n]; |
| 21 |
} |
| 22 |
|
| 23 |
mov_com /= (RealType)n_vec; |
| 24 |
|
| 25 |
mov_to_ref = ref_com - mov_com; |
| 26 |
|
| 27 |
/* shift mov to center of mass */ |
| 28 |
|
| 29 |
for (n=0; n < n_vec; n++) { |
| 30 |
mov[n] -= mov_com; |
| 31 |
} |
| 32 |
|
| 33 |
/* initialize */ |
| 34 |
Mat3x3d R = Mat3x3d(0.0); |
| 35 |
RealType E0 = 0.0; |
| 36 |
|
| 37 |
for (n=0; n < n_vec; n++) { |
| 38 |
|
| 39 |
/* |
| 40 |
* E0 = 1/2 * sum(over n): y(n)*y(n) + x(n)*x(n) |
| 41 |
*/ |
| 42 |
E0 += dot(mov[n], mov[n]) + dot(ref_[n], ref_[n]); |
| 43 |
|
| 44 |
/* |
| 45 |
* correlation matrix R: |
| 46 |
* R(i,j) = sum(over n): y(n,i) * x(n,j) |
| 47 |
* where x(n) and y(n) are two vector sets |
| 48 |
*/ |
| 49 |
|
| 50 |
R += outProduct(mov[n], ref_[n]); |
| 51 |
|
| 52 |
} |
| 53 |
E0 *= 0.5; |
| 54 |
|
| 55 |
RectMatrix<RealType, n_vec, 3> v; |
| 56 |
Vector3d s; |
| 57 |
Mat3x3d w; |
| 58 |
|
| 59 |
SVD<RealType, n_vec, 3> svd = SVD<RealType, n_vec, 3>(R); |
| 60 |
svd.getU(v); |
| 61 |
svd.getSingularValues(s); |
| 62 |
svd.getV(w); |
| 63 |
|
| 64 |
int is_reflection = (v.determinant() * w.determinant()) < 0.0; |
| 65 |
if (is_reflection) |
| 66 |
s(2) = -s(2); |
| 67 |
|
| 68 |
RealType rmsd_sq = (E0 - 2.0 * s.sum() )/ (RealType)n_vec; |
| 69 |
rmsd_sq = max(rmsd_sq,0.0); |
| 70 |
RealType rmsd = sqrt(rmsd_sq); |
| 71 |
return rmsd; |
| 72 |
} |
| 73 |
|
| 74 |
|
| 75 |
RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov) { |
| 76 |
} |
| 77 |
|
