| 1 |
#include "math/RMSD.hpp" |
| 2 |
#include "math/SVD.hpp" |
| 3 |
|
| 4 |
using namespace OpenMD; |
| 5 |
using namespace JAMA; |
| 6 |
|
| 7 |
RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov, |
| 8 |
Vector3d mov_com, |
| 9 |
Vector3d mov_to_ref) { |
| 10 |
|
| 11 |
assert(mov.size() == ref_.size()); |
| 12 |
int n; |
| 13 |
int n_vec = ref_.size(); |
| 14 |
|
| 15 |
/* calculate the centre of mass */ |
| 16 |
mov_com = V3Zero; |
| 17 |
|
| 18 |
for (n=0; n < n_vec; n++) { |
| 19 |
mov_com += mov[n]; |
| 20 |
} |
| 21 |
|
| 22 |
mov_com /= (RealType)n_vec; |
| 23 |
|
| 24 |
mov_to_ref = ref_com - mov_com; |
| 25 |
|
| 26 |
/* shift mov to center of mass */ |
| 27 |
|
| 28 |
for (n=0; n < n_vec; n++) { |
| 29 |
mov[n] -= mov_com; |
| 30 |
} |
| 31 |
|
| 32 |
/* initialize */ |
| 33 |
Mat3x3d R(0.0); |
| 34 |
RealType E0 = 0.0; |
| 35 |
|
| 36 |
for (n=0; n < n_vec; n++) { |
| 37 |
|
| 38 |
/* |
| 39 |
* E0 = 1/2 * sum(over n): y(n)*y(n) + x(n)*x(n) |
| 40 |
*/ |
| 41 |
E0 += dot(mov[n], mov[n]) + dot(ref_[n], ref_[n]); |
| 42 |
|
| 43 |
/* |
| 44 |
* correlation matrix R: |
| 45 |
* R(i,j) = sum(over n): y(n,i) * x(n,j) |
| 46 |
* where x(n) and y(n) are two vector sets |
| 47 |
*/ |
| 48 |
|
| 49 |
R += outProduct(mov[n], ref_[n]); |
| 50 |
|
| 51 |
} |
| 52 |
E0 *= 0.5; |
| 53 |
|
| 54 |
// Pack everything into Dynamic arrays for the SVD: |
| 55 |
DynamicRectMatrix<RealType> Rtmp(3,3,0.0); |
| 56 |
DynamicRectMatrix<RealType> vtmp(3, 3); |
| 57 |
DynamicVector<RealType> stmp(3); |
| 58 |
DynamicRectMatrix<RealType> wtmp(3, 3); |
| 59 |
|
| 60 |
Rtmp.setSubMatrix(0, 0, R); |
| 61 |
SVD<RealType> svd(Rtmp); |
| 62 |
|
| 63 |
svd.getU(vtmp); |
| 64 |
svd.getSingularValues(stmp); |
| 65 |
svd.getV(wtmp); |
| 66 |
|
| 67 |
Mat3x3d v; |
| 68 |
Vector3d s; |
| 69 |
Mat3x3d w; |
| 70 |
|
| 71 |
vtmp.getSubMatrix(0,0,v); |
| 72 |
stmp.getSubVector(0,s); |
| 73 |
wtmp.getSubMatrix(0,0,w); |
| 74 |
|
| 75 |
int is_reflection = (v.determinant() * w.determinant()) < RealType(0.0); |
| 76 |
if (is_reflection) |
| 77 |
s(2) = -s(2); |
| 78 |
|
| 79 |
RealType rmsd_sq = (E0 - 2.0 * s.sum() )/ (RealType)n_vec; |
| 80 |
rmsd_sq = max(rmsd_sq, RealType(0.0)); |
| 81 |
RealType rmsd = sqrt(rmsd_sq); |
| 82 |
return rmsd; |
| 83 |
} |
| 84 |
|
| 85 |
RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov, |
| 86 |
Vector3d mov_com, |
| 87 |
Vector3d mov_to_ref) { |
| 88 |
|
| 89 |
assert(mov.size() == ref_.size()); |
| 90 |
int n; |
| 91 |
int n_vec = ref_.size(); |
| 92 |
|
| 93 |
/* calculate the centre of mass */ |
| 94 |
mov_com = V3Zero; |
| 95 |
|
| 96 |
for (n=0; n < n_vec; n++) { |
| 97 |
mov_com += mov[n]; |
| 98 |
} |
| 99 |
|
| 100 |
mov_com /= (RealType)n_vec; |
| 101 |
|
| 102 |
mov_to_ref = ref_com - mov_com; |
| 103 |
|
| 104 |
/* shift mov to center of mass */ |
| 105 |
|
| 106 |
for (n=0; n < n_vec; n++) { |
| 107 |
mov[n] -= mov_com; |
| 108 |
} |
| 109 |
|
| 110 |
/* initialize */ |
| 111 |
Mat3x3d R(0.0); |
| 112 |
|
| 113 |
for (n=0; n < n_vec; n++) { |
| 114 |
|
| 115 |
/* |
| 116 |
* correlation matrix R: |
| 117 |
* R(i,j) = sum(over n): y(n,i) * x(n,j) |
| 118 |
* where x(n) and y(n) are two vector sets |
| 119 |
*/ |
| 120 |
|
| 121 |
R += outProduct(mov[n], ref_[n]); |
| 122 |
} |
| 123 |
|
| 124 |
// Pack everything into Dynamic arrays for the SVD: |
| 125 |
DynamicRectMatrix<RealType> Rtmp(3,3,0.0); |
| 126 |
DynamicRectMatrix<RealType> vtmp(3, 3); |
| 127 |
DynamicVector<RealType> stmp(3); |
| 128 |
DynamicRectMatrix<RealType> wtmp(3, 3); |
| 129 |
|
| 130 |
Rtmp.setSubMatrix(0, 0, R); |
| 131 |
SVD<RealType> svd(Rtmp); |
| 132 |
|
| 133 |
svd.getU(vtmp); |
| 134 |
svd.getSingularValues(stmp); |
| 135 |
svd.getV(wtmp); |
| 136 |
|
| 137 |
Mat3x3d v; |
| 138 |
Vector3d s; |
| 139 |
Mat3x3d w; |
| 140 |
|
| 141 |
vtmp.getSubMatrix(0,0,v); |
| 142 |
stmp.getSubVector(0,s); |
| 143 |
wtmp.getSubMatrix(0,0,w); |
| 144 |
|
| 145 |
int is_reflection = (v.determinant() * w.determinant()) < 0.0; |
| 146 |
if (is_reflection) |
| 147 |
s(2) = -s(2); |
| 148 |
|
| 149 |
RotMat3x3d U = v*w; |
| 150 |
return U; |
| 151 |
} |
| 152 |
|
