diff --git a/modules/index.rst b/modules/index.rst
index 8670525746c11f543886ed9108306aa756b2c6fa..147274735b6617024bf345b5d88c6ce4b031c8fa 100644
--- a/modules/index.rst
+++ b/modules/index.rst
@@ -11,6 +11,7 @@ OpenStructure documentation
base/base
geom/geom
mol/base/mol
+ mol/alg/molalg
conop/conop
img/base/img
img/alg/alg
diff --git a/modules/mol/alg/molalg.rst b/modules/mol/alg/molalg.rst
new file mode 100644
index 0000000000000000000000000000000000000000..00735f6def328ed899de4005fb2f803cd23f0206
--- /dev/null
+++ b/modules/mol/alg/molalg.rst
@@ -0,0 +1,25 @@
+:mod:`mol.alg <ost.mol.alg>` -- Algorithms for Structures
+================================================================================
+
+.. module:: ost.mol.alg
+ :synopsis: Algorithms operating on molecular structures
+
+.. function:: LocalDistTest(model, reference, tolerance, radius)
+
+ This function calculates the agreement of local contacts between the model and
+ the reference structure. The overlap is a number between zero and one, where
+ one indicates complete agreement, zero indicates no agreement at all. This
+ score is similar to the GDT, but does not require any superposition of the
+ model and the reference.
+
+ The distance of atom pairs in the reference structure that are closer than a
+ certain distance (radius) to each other is compared to the distance found in
+ the model. If the difference between these two distances is smaller than a
+ threshold value (tolerance), the model and the reference agree. Missing atoms
+ in the model are treated disagreement and thus lower the score.
+
+ For residue with symmetric sidechains (GLU, ASP, ARG, VAL, PHE, TYR), the
+ naming of the atoms is ambigous. For these residues, the overlap of both
+ possible solutions to the fixed atoms, that is, everything that is not
+ ambigous is calculated. The solution that gives higher overlap is then used to
+ calculate the actual overlap score.
\ No newline at end of file
diff --git a/modules/mol/alg/pymod/wrap_mol_alg.cc b/modules/mol/alg/pymod/wrap_mol_alg.cc
index ce0be429c6afb0afd5b6b55b2252cc79d446dbbd..240a029a42dd523e8a4dfab72be40fe77bd74fa7 100644
--- a/modules/mol/alg/pymod/wrap_mol_alg.cc
+++ b/modules/mol/alg/pymod/wrap_mol_alg.cc
@@ -22,6 +22,8 @@
*/
#include <boost/python.hpp>
#include <ost/config.hh>
+#include <ost/mol/alg/local_dist_test.hh>
+
using namespace boost::python;
void export_svdSuperPose();
@@ -36,4 +38,6 @@ BOOST_PYTHON_MODULE(_mol_alg)
#if OST_IMG_ENABLED
export_entity_to_density();
#endif
+
+ def("LocalDistTest", &ost::mol::alg::LocalDistTest);
}
diff --git a/modules/mol/alg/src/CMakeLists.txt b/modules/mol/alg/src/CMakeLists.txt
index 242eea3d8f3a4c80f29909b3ecad313df10753f7..17dc9acbf2a5d32ac17e62b0e50af18f35122b2f 100644
--- a/modules/mol/alg/src/CMakeLists.txt
+++ b/modules/mol/alg/src/CMakeLists.txt
@@ -2,11 +2,13 @@ set(OST_MOL_ALG_HEADERS
svd_superpose.hh
module_config.hh
sec_structure_segments.hh
+ local_dist_test.hh
)
set(OST_MOL_ALG_SOURCES
svd_superpose.cc
sec_structure_segments.cc
+ local_dist_test.cc
)
set(MOL_ALG_DEPS mol)
diff --git a/modules/mol/alg/src/local_dist_test.cc b/modules/mol/alg/src/local_dist_test.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a0c07bf58190170938e2a33340da4834837d042
--- /dev/null
+++ b/modules/mol/alg/src/local_dist_test.cc
@@ -0,0 +1,162 @@
+#include <ost/mol/mol.hh>
+#include "local_dist_test.hh"
+
+namespace ost { namespace mol { namespace alg {
+
+namespace {
+
+String swapped_name(const String& name)
+{
+ if (name=="OE1") return "OE2";
+ if (name=="OE2") return "OE1";
+
+ if (name=="OD1") return "OD2";
+ if (name=="OD2") return "OD1";
+
+ if (name=="CG1") return "CG2";
+ if (name=="CG2") return "CG1";
+
+ if (name=="CE1") return "CE2";
+ if (name=="CE2") return "CE1";
+
+ if (name=="CD1") return "CD2";
+ if (name=="CD2") return "CD1";
+
+ if (name=="NH1") return "NH2";
+ if (name=="NH2") return "NH1";
+
+ return name;
+}
+
+bool swappable(const String& rname, const String& aname)
+{
+ if (rname=="GLU") {
+ return (aname=="OE1" || aname=="OE2");
+ }
+ if (rname=="ASP") {
+ return (aname=="OD1" || aname=="OD2");
+ }
+ if (rname=="VAL") {
+ return (aname=="CG1" || aname=="CG2");
+ }
+ if (rname=="TYR" || rname=="PHE") {
+ return (aname=="CD1" || aname=="CD2" || aname=="CE1" || aname=="CE2");
+ }
+ if (rname=="ARG") {
+ return (aname=="NH1" || aname=="NH2");
+ }
+ return false;
+}
+
+std::pair<Real, Real> calc_overlap(ResidueView ref_res,
+ EntityView ref,
+ ChainView mdl_chain,
+ Real tol, Real max_dist, bool only_fixed,
+ bool swap)
+{
+ std::pair<Real, Real> overlap(0.0, 0.0);
+ AtomViewList ref_atoms=ref_res.GetAtomList();
+
+ ResidueView mdl_res=mdl_chain.FindResidue(ref_res.GetNumber());
+ for (AtomViewList::iterator ai=ref_atoms.begin(),
+ ae=ref_atoms.end(); ai!=ae; ++ai) {
+ if (ai->GetElement()=="H") { continue; }
+ String name=swap ? swapped_name(ai->GetName()) : ai->GetName();
+ AtomView av1=mdl_res ? mdl_res.FindAtom(name) : AtomView();
+ AtomViewList within=ref.FindWithin(ai->GetPos(), max_dist);
+ for (AtomViewList::iterator aj=within.begin(),
+ ae2=within.end(); aj!=ae2; ++aj) {
+ if (aj->GetElement()=="H") { continue; }
+ if (only_fixed) {
+ if (aj->GetResidue().GetNumber()==ref_res.GetNumber()) {
+ continue;
+ }
+ if (swappable(aj->GetResidue().GetName(), aj->GetName())) {
+ continue;
+ }
+ AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(),
+ aj->GetName());
+ overlap.second+=1.0;
+ if (!(av1 && av2)) {
+ continue;
+ }
+ Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+ Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
+ if (std::abs(mdl_dist-ref_dist)<tol) {
+ overlap.first+=1;
+ }
+ continue;
+ }
+ if (aj->GetResidue().GetNumber()>ref_res.GetNumber()) {
+ if (aj->GetResidue().GetNumber()<=ref_res.GetNumber()) {
+ continue;
+ }
+ AtomView av2=mdl_chain.FindAtom(aj->GetResidue().GetNumber(),
+ aj->GetName());
+ overlap.second+=1.0;
+ if (!(av1 && av2)) {
+ continue;
+ }
+ Real mdl_dist=geom::Length(av1.GetPos()-av2.GetPos());
+ Real ref_dist=geom::Length(ai->GetPos()-aj->GetPos());
+ if (std::abs(mdl_dist-ref_dist)<tol) {
+ overlap.first+=1;
+ }
+ }
+ }
+ }
+ return overlap;
+}
+
+}
+
+Real LocalDistTest(const EntityView& mdl, const EntityView& ref,
+ Real cutoff, Real max_dist)
+{
+ ResidueViewList ref_residues=ref.GetResidueList();
+ ChainView mdl_chain=mdl.GetChainList()[0];
+ // Residues with symmetric side-chains require special treatment as there are
+ // two possible ways to name the atoms. We calculate the overlap with the
+ // fixed atoms and take the solution that gives bigger scores.
+ XCSEditor edi=ref.GetHandle().RequestXCSEditor(BUFFERED_EDIT);
+ for (ResidueViewList::iterator
+ i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
+ const String rname=i->GetName();
+ ResidueView mdl_res=mdl_chain.FindResidue(i->GetNumber());
+ if (!mdl_res) {
+ continue;
+ }
+ if (!(rname=="GLU" || rname=="ASP" || rname=="VAL" || rname=="TYR" ||
+ rname=="PHE" || rname=="LYS" || rname=="ARG")) {
+ continue;
+ }
+
+ std::pair<Real, Real> ov1=calc_overlap(*i, ref, mdl_chain,
+ cutoff, max_dist, true,
+ false);
+ std::pair<Real, Real> ov2=calc_overlap(*i, ref, mdl_chain,
+ cutoff, max_dist, true,
+ true);
+ if (ov1.first/ov1.second<ov2.first/ov2.second) {
+ AtomViewList atoms=mdl_res.GetAtomList();
+ for (AtomViewList::iterator j=atoms.begin(),
+ e2=atoms.end(); j!=e2; ++j) {
+ if (swappable(rname, j->GetName())) {
+ edi.RenameAtom(j->GetHandle(), swapped_name(j->GetName()));
+ }
+ }
+
+ }
+ }
+ std::pair<Real, Real> total_ov(0.0, 0.0);
+ for (ResidueViewList::iterator
+ i=ref_residues.begin(), e=ref_residues.end(); i!=e; ++i) {
+ std::pair<Real, Real> ov1=calc_overlap(*i, ref, mdl_chain, cutoff,
+ max_dist, false, false);
+ total_ov.first+=ov1.first;
+ total_ov.second+=ov1.second;
+ }
+ return total_ov.first/total_ov.second;
+}
+
+}}}
diff --git a/modules/mol/alg/src/local_dist_test.hh b/modules/mol/alg/src/local_dist_test.hh
new file mode 100644
index 0000000000000000000000000000000000000000..f0ba3ebea20a6138a0cf14f2a34d24839abe3ef4
--- /dev/null
+++ b/modules/mol/alg/src/local_dist_test.hh
@@ -0,0 +1,15 @@
+#ifndef OST_MOL_ALG_LOCAL_DIST_TEST_HH
+#define OST_MOL_ALG_LOCAL_DIST_TEST_HH
+
+#include <ost/mol/entity_view.hh>
+#include <ost/mol/alg/module_config.hh>
+
+namespace ost { namespace mol { namespace alg {
+
+Real DLLEXPORT_OST_MOL_ALG LocalDistTest(const EntityView& mdl,
+ const EntityView& ref,
+ Real cutoff, Real max_dist);
+}}}
+
+#endif
+