From 7436c4eba1ba0f2370d426fe20587d29089266df Mon Sep 17 00:00:00 2001
From: Pavel <pafonine@lbl.gov>
Date: Fri, 22 Sep 2023 16:55:59 -0700
Subject: [PATCH] Bug fix in phenix.douse (Reported by Rob).

---
 mmtbx/solvent/map_to_water.py |  4 +--
 mmtbx/utils/utils.h           | 60 +++++++++++++++++------------------
 mmtbx/utils/utils_ext.cpp     |  4 +--
 3 files changed, 33 insertions(+), 35 deletions(-)

diff --git a/mmtbx/solvent/map_to_water.py b/mmtbx/solvent/map_to_water.py
index fb57e0def3..cdcbf56bbc 100644
--- a/mmtbx/solvent/map_to_water.py
+++ b/mmtbx/solvent/map_to_water.py
@@ -502,8 +502,8 @@ def _filter_by_distance(self):
       "  distance (A), min: %4.2f max: %4.2f"%(self.dist_min, self.dist_max))
     self.ma.add("  start: %d"%self.xrs_water.scatterers().size())
     sel = mmtbx.utils.filter_water(
-      sites_frac_interaction = self.sites_frac_interaction,
-      sites_frac_other       = self.model.get_sites_frac().select(~self.interaction_selection),
+      interaction_selection  = self.interaction_selection,
+      sites_frac_other       = self.model.get_sites_frac(),
       sites_frac_water       = self.xrs_water.sites_frac(),
       dist_min               = self.dist_min,
       dist_max               = self.dist_max,
diff --git a/mmtbx/utils/utils.h b/mmtbx/utils/utils.h
index 49801ae38c..95b39ae772 100644
--- a/mmtbx/utils/utils.h
+++ b/mmtbx/utils/utils.h
@@ -141,7 +141,7 @@ class density_distribution_per_atom
 template <typename FloatType>
 af::shared<std::size_t>
   filter_water(
-    af::shared<vec3<FloatType> > const& sites_frac_interaction,
+    af::shared<bool> const& interaction_selection,
     af::shared<vec3<FloatType> > const& sites_frac_other,
     af::shared<vec3<FloatType> > const& sites_frac_water,
     FloatType const& dist_max,
@@ -149,62 +149,60 @@ af::shared<std::size_t>
     cctbx::uctbx::unit_cell const& unit_cell)
 {
   af::shared<std::size_t> result;
-  af::shared<std::size_t> result_;
   af::shared<std::size_t> first_shell;
   af::shared<std::size_t> second_shell;
+  MMTBX_ASSERT(interaction_selection.size()==sites_frac_other.size());
+  // Select water that a) do not clash with any of non-water atom and b) are
+  // within prescribed limits (first_shell).
+  // Set aside others non-clashing (second_shell).
   for(std::size_t i=0; i<sites_frac_water.size(); i+=1) {
     FloatType dist_closest = 1.e+9;
+    FloatType dist_closest_int = 1.e+9;
     cctbx::fractional<> sfw = sites_frac_water[i];
-    for(std::size_t j=0; j<sites_frac_interaction.size(); j+=1) {
-      cctbx::fractional<> sf = sites_frac_interaction[j];
-      FloatType dist = unit_cell.distance(sf, sfw);
-      if(dist < dist_closest) {
-        dist_closest = dist;
+    bool skip = false;
+    for(std::size_t j=0; j<sites_frac_other.size(); j+=1) {
+      cctbx::fractional<> sfo = sites_frac_other[j];
+      FloatType dist = unit_cell.distance(sfo, sfw);
+      if(dist<dist_min) { // clash detected
+        skip = true;
+        break;
       }
+      if(dist<dist_closest) dist_closest = dist;
+      if(dist<dist_closest_int && interaction_selection[j]) dist_closest_int = dist;
     }
-    if(dist_closest<=dist_max && dist_closest>=dist_min) {
+    if(skip) continue;
+    if(dist_closest_int<=dist_max &&
+       dist_closest_int>=dist_min &&
+       dist_closest    >=dist_min) {
       first_shell.push_back(i);
     }
     else {
       second_shell.push_back(i);
     }
   }
+  // Now check those set aside (second_shell) and fill into result
   for(std::size_t i=0; i<second_shell.size(); i+=1) {
     FloatType dist_closest = 1.e+9;
     cctbx::fractional<> sfi = sites_frac_water[second_shell[i]];
+    bool skip = false;
     for(std::size_t j=0; j<first_shell.size(); j+=1) {
       cctbx::fractional<> sfj = sites_frac_water[first_shell[j]];
       FloatType dist = unit_cell.distance(sfi, sfj);
-      if(dist < dist_closest) {
-        dist_closest = dist;
+      if(dist<dist_min) {
+        break; // clash detected
+        skip = true;
       }
+      if(dist < dist_closest) dist_closest = dist;
     }
+    if(skip) continue;
     if(dist_closest<=dist_max && dist_closest>=dist_min) {
-      result_.push_back(second_shell[i]);
+      result.push_back(second_shell[i]);
     }
   }
+  // Add first shell into result
   for(std::size_t i=0; i<first_shell.size(); i+=1) {
-    result_.push_back(first_shell[i]);
+    result.push_back(first_shell[i]);
   }
-
-
-  for(std::size_t i=0; i<result_.size(); i+=1) {
-    FloatType dist_closest = 1.e+9;
-    cctbx::fractional<> sfi = sites_frac_water[result_[i]];
-
-    for(std::size_t j=0; j<sites_frac_other.size(); j+=1) {
-      cctbx::fractional<> sfj = sites_frac_other[j];
-      FloatType dist = unit_cell.distance(sfi, sfj);
-      if(dist < dist_closest) {
-        dist_closest = dist;
-      }
-    }
-    if(dist_closest>=dist_min) {
-      result.push_back(result_[i]);
-    }
-  }
-
-
   return result;
 }
 
diff --git a/mmtbx/utils/utils_ext.cpp b/mmtbx/utils/utils_ext.cpp
index 35908701ed..d64744789f 100644
--- a/mmtbx/utils/utils_ext.cpp
+++ b/mmtbx/utils/utils_ext.cpp
@@ -56,13 +56,13 @@ namespace {
 
     def("filter_water",
          (af::shared<std::size_t>(*)
-               (af::shared<vec3<double> > const&,
+               (af::shared<bool> const&,
                 af::shared<vec3<double> > const&,
                 af::shared<vec3<double> > const&,
                 double const&,
                 double const&,
                 cctbx::uctbx::unit_cell const&)) filter_water,
-                  (arg("sites_frac_interaction"),
+                  (arg("interaction_selection"),
                    arg("sites_frac_other"),
                    arg("sites_frac_water"),
                    arg("dist_max"),