diff --git a/modules/mol/alg/pymod/scoring.py b/modules/mol/alg/pymod/scoring.py
index 4444e463709703042b6395da3c866804cd971237..d415b0aad53c7f9327d2f7c47191a8e4d28e1940 100644
--- a/modules/mol/alg/pymod/scoring.py
+++ b/modules/mol/alg/pymod/scoring.py
@@ -381,6 +381,7 @@ class Scorer:
self._rigid_n_target_not_mapped = None
self._rigid_transform = None
+ self._gdt_window_sizes = [5, 7, 9, 12, 24, 36, 48]
self._gdt_05 = None
self._gdt_1 = None
self._gdt_2 = None
@@ -1399,7 +1400,13 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_05 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 0.5)
+ N = list()
+ for window_size in self._gdt_window_sizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 0.5)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_05 = float(n) / n_full
@@ -1417,7 +1424,13 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_1 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 1.0)
+ N = list()
+ for window_size in self._gdt_window_sizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 1.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_1 = float(n) / n_full
@@ -1436,7 +1449,13 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_2 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 2.0)
+ N = list()
+ for window_size in self._gdt_window_sizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 2.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_2 = float(n) / n_full
@@ -1454,7 +1473,13 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_4 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 4.0)
+ N = list()
+ for window_size in self._gdt_window_sizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 4.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_4 = float(n) / n_full
@@ -1471,7 +1496,13 @@ class Scorer:
:type: :class:`float`
"""
if self._gdt_8 is None:
- n, m = GDT(self.rigid_mapped_model_pos, self.rigid_mapped_target_pos, 7, 1000, 8.0)
+ N = list()
+ for window_size in self._gdt_window_sizes:
+ n = GDT(self.rigid_mapped_model_pos,
+ self.rigid_mapped_target_pos,
+ window_size, 1000, 8.0)[0]
+ N.append(n)
+ n = max(N)
n_full = len(self.rigid_mapped_target_pos) + self.rigid_n_target_not_mapped
if n_full > 0:
self._gdt_8 = float(n) / n_full
diff --git a/modules/mol/alg/src/gdt.cc b/modules/mol/alg/src/gdt.cc
index 93b6e8b173b51b04fb4c0e01c092f661b7507d04..8eedd7796241085c89ec8ed15d8f28ca168afd03 100644
--- a/modules/mol/alg/src/gdt.cc
+++ b/modules/mol/alg/src/gdt.cc
@@ -302,8 +302,126 @@ void GDT(const geom::Vec3List& mdl_pos, const geom::Vec3List& ref_pos,
int n_pos = mdl_pos.size();
+ // deal with special cases that don't produce valid transforms
+ if(n_pos == 1) {
+ transform = geom::Mat4::Identity();
+ transform.PasteTranslation(ref_pos[0] - mdl_pos[0]);
+ n_superposed = 1;
+ return;
+ }
+
+ if(n_pos == 2) {
+ Real mdl_d = geom::Distance(mdl_pos[0], mdl_pos[1]);
+ Real ref_d = geom::Distance(ref_pos[0], ref_pos[1]);
+ Real dd = std::abs(mdl_d - ref_d);
+ if(dd/2 <= distance_thresh) {
+ // the two can be superposed within specified distance threshold
+ // BUT: cannot construct valid transformation from two positions
+ // => Construct matrix with four positions
+ // Two are constructed starting from the center point +- some direction
+ // vector that is orthogonal to the vector connecting the original two
+ // points.
+ geom::Vec3 mdl_center = (mdl_pos[0] + mdl_pos[1])*0.5;
+ geom::Vec3 ref_center = (ref_pos[0] + ref_pos[1])*0.5;
+ Eigen::Matrix<double, 4, 3> eigen_mdl_pos = \
+ Eigen::Matrix<double, 4, 3>::Zero(4, 3);
+ Eigen::Matrix<double, 4, 3> eigen_ref_pos = \
+ Eigen::Matrix<double, 4, 3>::Zero(4, 3);
+ Eigen::Matrix<double,3,3> eigen_rot = \
+ Eigen::Matrix<double,3,3>::Identity();
+
+ geom::Vec3 mdl_dir = geom::Normalize(mdl_pos[1] - mdl_pos[0]);
+ geom::Vec3 ref_dir = geom::Normalize(ref_pos[1] - ref_pos[0]);
+ geom::Vec3 mdl_normal;
+ geom::Vec3 ref_normal;
+
+ // Use cross product to get some normal on mdl_dir
+ // The direction of the second vector doesn't really matter, but shouldnt
+ // be collinear with mdl_dir
+ if(mdl_dir[0] < 0.999) {
+ mdl_normal = geom::Cross(geom::Vec3(1,0,0), mdl_dir);
+ } else {
+ mdl_normal = geom::Cross(geom::Vec3(0,1,0), mdl_dir);
+ }
+
+ // same for ref_dir
+ if(ref_dir[0] < 0.999) {
+ ref_normal = geom::Cross(geom::Vec3(1,0,0), ref_dir);
+ } else {
+ ref_normal = geom::Cross(geom::Vec3(0,1,0), ref_dir);
+ }
+
+ eigen_mdl_pos(0, 0) = mdl_pos[0][0] - mdl_center[0];
+ eigen_mdl_pos(0, 1) = mdl_pos[0][1] - mdl_center[1];
+ eigen_mdl_pos(0, 2) = mdl_pos[0][2] - mdl_center[2];
+ eigen_mdl_pos(1, 0) = mdl_pos[1][0] - mdl_center[0];
+ eigen_mdl_pos(1, 1) = mdl_pos[1][1] - mdl_center[1];
+ eigen_mdl_pos(1, 2) = mdl_pos[1][2] - mdl_center[2];
+ eigen_mdl_pos(2, 0) = mdl_normal[0];
+ eigen_mdl_pos(2, 1) = mdl_normal[1];
+ eigen_mdl_pos(2, 2) = mdl_normal[2];
+ eigen_mdl_pos(3, 0) = -mdl_normal[0];
+ eigen_mdl_pos(3, 1) = -mdl_normal[1];
+ eigen_mdl_pos(3, 2) = -mdl_normal[2];
+ eigen_ref_pos(0, 0) = ref_pos[0][0] - ref_center[0];
+ eigen_ref_pos(0, 1) = ref_pos[0][1] - ref_center[1];
+ eigen_ref_pos(0, 2) = ref_pos[0][2] - ref_center[2];
+ eigen_ref_pos(1, 0) = ref_pos[1][0] - ref_center[0];
+ eigen_ref_pos(1, 1) = ref_pos[1][1] - ref_center[1];
+ eigen_ref_pos(1, 2) = ref_pos[1][2] - ref_center[2];
+ eigen_ref_pos(2, 0) = ref_normal[0];
+ eigen_ref_pos(2, 1) = ref_normal[1];
+ eigen_ref_pos(2, 2) = ref_normal[2];
+ eigen_ref_pos(3, 0) = -ref_normal[0];
+ eigen_ref_pos(3, 1) = -ref_normal[1];
+ eigen_ref_pos(3, 2) = -ref_normal[2];
+
+ Real tmp; // no need to store RMSD
+ TheobaldRMSD(eigen_mdl_pos, eigen_ref_pos, tmp, eigen_rot);
+ transform = geom::Mat4();
+ transform(0, 0) = eigen_rot(0, 0);
+ transform(0, 1) = eigen_rot(0, 1);
+ transform(0, 2) = eigen_rot(0, 2);
+ transform(1, 0) = eigen_rot(1, 0);
+ transform(1, 1) = eigen_rot(1, 1);
+ transform(1, 2) = eigen_rot(1, 2);
+ transform(2, 0) = eigen_rot(2, 0);
+ transform(2, 1) = eigen_rot(2, 1);
+ transform(2, 2) = eigen_rot(2, 2);
+
+ // there are three transformation to be applied to reach ref_pos from
+ // mdl_pos:
+ // 1: shift mdl_pos to center
+ // 2: apply estimated rotation
+ // 3: shift onto average of ref_pos
+ Eigen::Matrix<double,1,3> eigen_avg_mdl = Eigen::Matrix<double,1,3>::Zero();
+ Eigen::Matrix<double,1,3> eigen_avg_ref = Eigen::Matrix<double,1,3>::Zero();
+ eigen_avg_mdl(0,0) = mdl_center[0];
+ eigen_avg_mdl(0,1) = mdl_center[1];
+ eigen_avg_mdl(0,2) = mdl_center[2];
+ eigen_avg_ref(0,0) = ref_center[0];
+ eigen_avg_ref(0,1) = ref_center[1];
+ eigen_avg_ref(0,2) = ref_center[2];
+ Eigen::Matrix<double,1,3> translation = eigen_rot *
+ (-eigen_avg_mdl.transpose()) +
+ eigen_avg_ref.transpose();
+ transform(0, 3) = translation(0, 0);
+ transform(1, 3) = translation(0, 1);
+ transform(2, 3) = translation(0, 2);
+ n_superposed = 2;
+ } else {
+ // the two cannot be superposed within specified distance threshold
+ // => just set n_superposed to 1 and generate the same transformation
+ // as in n_pos == 1
+ transform = geom::Mat4::Identity();
+ transform.PasteTranslation(ref_pos[0] - mdl_pos[0]);
+ n_superposed = 1;
+ }
+ return;
+ }
+
if(window_size > n_pos) {
- throw ost::Error("Window size in GDT algorithm is larger than positions");
+ window_size = n_pos;
}
Eigen::Matrix<double, Eigen::Dynamic, 3> eigen_mdl_pos = \