| 1 |
|
#include "math/RMSD.hpp" |
| 2 |
|
#include "math/SVD.hpp" |
| 3 |
|
|
| 4 |
< |
using namespace oopse; |
| 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 |
< |
RotMat3x3d U) { |
| 11 |
< |
|
| 9 |
> |
Vector3d mov_to_ref) { |
| 10 |
> |
|
| 11 |
|
assert(mov.size() == ref_.size()); |
| 12 |
|
int n; |
| 13 |
|
int n_vec = ref_.size(); |
| 30 |
|
} |
| 31 |
|
|
| 32 |
|
/* initialize */ |
| 33 |
< |
Mat3x3d R = Mat3x3d(0.0); |
| 33 |
> |
Mat3x3d R(0.0); |
| 34 |
|
RealType E0 = 0.0; |
| 35 |
|
|
| 36 |
|
for (n=0; n < n_vec; n++) { |
| 51 |
|
} |
| 52 |
|
E0 *= 0.5; |
| 53 |
|
|
| 54 |
< |
RectMatrix<RealType, n_vec, 3> v; |
| 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 |
< |
SVD<RealType, n_vec, 3> svd = SVD<RealType, n_vec, 3>(R); |
| 72 |
< |
svd.getU(v); |
| 73 |
< |
svd.getSingularValues(s); |
| 62 |
< |
svd.getV(w); |
| 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()) < 0.0; |
| 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,0.0); |
| 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 |
< |
RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov) { |
| 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 |
|
|
