1#include "math/RMSD.hpp"
8RealType RMSD::calculate_rmsd(std::vector<Vector3d> mov,
Vector3d mov_com,
10 assert(mov.size() == ref_.size());
12 int n_vec = ref_.size();
17 for (n = 0; n < n_vec; n++) {
21 mov_com /= (RealType)n_vec;
23 mov_to_ref = ref_com - mov_com;
27 for (n = 0; n < n_vec; n++) {
35 for (n = 0; n < n_vec; n++) {
39 E0 +=
dot(mov[n], mov[n]) +
dot(ref_[n], ref_[n]);
47 R += outProduct(mov[n], ref_[n]);
52 DynamicRectMatrix<RealType> Rtmp(3, 3, 0.0);
53 DynamicRectMatrix<RealType> vtmp(3, 3);
54 DynamicVector<RealType> stmp(3);
55 DynamicRectMatrix<RealType> wtmp(3, 3);
57 Rtmp.setSubMatrix(0, 0, R);
58 SVD<RealType> svd(Rtmp);
61 svd.getSingularValues(stmp);
68 vtmp.getSubMatrix(0, 0, v);
69 stmp.getSubVector(0, s);
70 wtmp.getSubMatrix(0, 0, w);
72 int is_reflection = (v.determinant() * w.
determinant()) < RealType(0.0);
73 if (is_reflection) s(2) = -s(2);
75 RealType rmsd_sq = (E0 - 2.0 * s.
sum()) / (RealType)n_vec;
76 rmsd_sq = max(rmsd_sq, RealType(0.0));
77 RealType rmsd = sqrt(rmsd_sq);
81RotMat3x3d RMSD::optimal_superposition(std::vector<Vector3d> mov,
83 assert(mov.size() == ref_.size());
85 int n_vec = ref_.size();
90 for (n = 0; n < n_vec; n++) {
94 mov_com /= (RealType)n_vec;
96 mov_to_ref = ref_com - mov_com;
100 for (n = 0; n < n_vec; n++) {
107 for (n = 0; n < n_vec; n++) {
114 R += outProduct(mov[n], ref_[n]);
118 DynamicRectMatrix<RealType> Rtmp(3, 3, 0.0);
119 DynamicRectMatrix<RealType> vtmp(3, 3);
120 DynamicVector<RealType> stmp(3);
121 DynamicRectMatrix<RealType> wtmp(3, 3);
123 Rtmp.setSubMatrix(0, 0, R);
124 SVD<RealType> svd(Rtmp);
127 svd.getSingularValues(stmp);
134 vtmp.getSubMatrix(0, 0, v);
135 stmp.getSubVector(0, s);
136 wtmp.getSubMatrix(0, 0, w);
138 int is_reflection = (v.determinant() * w.
determinant()) < 0.0;
139 if (is_reflection) s(2) = -s(2);
Real determinant() const
Returns the determinant of this matrix.
Real sum()
Returns the sum of all elements of this vector.
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.