Skipping the setup of the request list if the NL if stride is 1

plugins/pycv/src/PythonCVInterface.cpp

@@ -606,22 +606,7 @@ PythonCVInterface::PythonCVInterface(const ActionOptions&ao) try ://the catch on
 void PythonCVInterface::prepare() {
   try {
     if (nl) {
-      if (nl->getStride() > 0) {
-        if (firsttime || (getStep() % nl->getStride() == 0)) {
-          requestAtoms(nl->getFullAtomList());
-          invalidateList = true;
-          firsttime = false;
-        } else {
-          requestAtoms(nl->getReducedAtomList());
-          invalidateList = false;
-          if (getExchangeStep())
-            error("Neighbor lists should be updated on exchange steps - choose a "
-                  "NL_STRIDE which divides the exchange stride!");
-        }
-        if (getExchangeStep()) {
-          firsttime = true;
-        }
-      }
+      std::tie(firsttime,invalidateList) = nl->prepare(this,firsttime, invalidateList).get();
     }
     if (hasPrepare) {
       py::gil_scoped_acquire gil;

src/colvar/CoordinationBase.cpp

@@ -135,22 +135,7 @@ CoordinationBase::~CoordinationBase() {
 }
 
 void CoordinationBase::prepare() {
-  if(nl->getStride()>0) {
-    if(firsttime || (getStep()%nl->getStride()==0)) {
-      requestAtoms(nl->getFullAtomList());
-      invalidateList=true;
-      firsttime=false;
-    } else {
-      requestAtoms(nl->getReducedAtomList());
-      invalidateList=false;
-      if(getExchangeStep()) {
-        error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
-      }
-    }
-    if(getExchangeStep()) {
-      firsttime=true;
-    }
-  }
+  std::tie(firsttime,invalidateList) =nl->prepare(this,firsttime, invalidateList).get();
 }
 
 // calculator

src/core/Colvar.h

@@ -36,10 +36,12 @@ namespace PLMD {
 This is the abstract base class to use for implementing new collective variables, within it there is
 \ref AddingAColvar "information" as to how to go about implementing a new CV.
 */
-
+class NeighborList;
 class Colvar :
   public ActionAtomistic,
   public ActionWithValue {
+
+  friend class NeighborList;
 private:
 protected:
   void requestAtoms(const std::vector<AtomNumber> & a);

src/maze/Optimizer.cpp

@@ -387,28 +387,7 @@ Vector Optimizer::center_of_mass() const {
 }
 
 void Optimizer::prepare() {
-  if (neighbor_list_->getStride() > 0) {
-    if (first_time_ || (getStep() % neighbor_list_->getStride() == 0)) {
-      requestAtoms(neighbor_list_->getFullAtomList());
-
-      validate_list_ = true;
-      first_time_ = false;
-    } else {
-      requestAtoms(neighbor_list_->getReducedAtomList());
-
-      validate_list_ = false;
-
-      if (getExchangeStep()) {
-        plumed_merror(
-          "maze> Neighbor lists should be updated on exchange steps -- choose "
-          "an NL_STRIDE which divides the exchange stride.\n");
-      }
-    }
-
-    if (getExchangeStep()) {
-      first_time_ = true;
-    }
-  }
+  std::tie(first_time_,validate_list_) = neighbor_list_->prepare(this,first_time_, validate_list_).get();
 }
 
 double Optimizer::score() {

src/s2cm/S2ContactModel.cpp

@@ -248,22 +248,7 @@ S2ContactModel::~S2ContactModel() {
 }
 
 void S2ContactModel::prepare() {
-  if(nl->getStride()>0) {
-    if(firsttime || (getStep()%nl->getStride()==0)) {
-      requestAtoms(nl->getFullAtomList());
-      invalidateList=true;
-      firsttime=false;
-    } else {
-      requestAtoms(nl->getReducedAtomList());
-      invalidateList=false;
-      if(getExchangeStep()) {
-        error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
-      }
-    }
-    if(getExchangeStep()) {
-      firsttime=true;
-    }
-  }
+  std::tie(firsttime,invalidateList) = nl->prepare(this,firsttime,invalidateList).get();
 }
 
 // calculator

src/tools/NeighborList.cpp

@@ -27,6 +27,7 @@
 #include "Communicator.h"
 #include "OpenMP.h"
 #include "Tools.h"
+#include "core/Colvar.h"
 #include <string>
 #include "LinkCells.h"
 #include "View.h"
@@ -296,11 +297,17 @@ void NeighborList::update(const std::vector<Vector>& positions) {
       }
     }
   }
-
-  setRequestList();
+  if (stride_ >1) {
+    setRequestList();
+  } else {
+    reduced=true;
+  }
 }
 
 void NeighborList::setRequestList() {
+  // at time of adding the `if (stride_>1)` in `update()`
+  // this function is called only from `update()` and it is private
+  // so as now it is not necessary to add extra logic in this function
   requestlist_.clear();
   for(unsigned int i=0; i<size(); ++i) {
     requestlist_.push_back(fullatomlist_[neighbors_[i].first]);
@@ -311,25 +318,30 @@ void NeighborList::setRequestList() {
 }
 
 std::vector<AtomNumber>& NeighborList::getReducedAtomList() {
-  if(!reduced) {
-    for(unsigned int i=0; i<size(); ++i) {
-      AtomNumber index0=fullatomlist_[neighbors_[i].first];
-      AtomNumber index1=fullatomlist_[neighbors_[i].second];
+  if(stride_>1) {
+    if(!reduced) {
+      for(unsigned int i=0; i<size(); ++i) {
+        AtomNumber index0=fullatomlist_[neighbors_[i].first];
+        AtomNumber index1=fullatomlist_[neighbors_[i].second];
 // I exploit the fact that requestlist_ is an ordered vector
 // And I assume that index0 and index1 actually exists in the requestlist_ (see setRequestList())
 // so I can use lower_bond that uses binary seach instead of find
-      plumed_dbg_assert(std::is_sorted(requestlist_.begin(),requestlist_.end()));
-      auto p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index0);
-      plumed_dbg_assert(p!=requestlist_.end());
-      unsigned newindex0=p-requestlist_.begin();
-      p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index1);
-      plumed_dbg_assert(p!=requestlist_.end());
-      unsigned newindex1=p-requestlist_.begin();
-      neighbors_[i]=pairIDs(newindex0,newindex1);
+        plumed_dbg_assert(std::is_sorted(requestlist_.begin(),requestlist_.end()));
+        auto p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index0);
+        plumed_dbg_assert(p!=requestlist_.end());
+        unsigned newindex0=p-requestlist_.begin();
+        p = std::lower_bound(requestlist_.begin(), requestlist_.end(), index1);
+        plumed_dbg_assert(p!=requestlist_.end());
+        unsigned newindex1=p-requestlist_.begin();
+        neighbors_[i]=pairIDs(newindex0,newindex1);
+      }
     }
+    reduced=true;
+    return requestlist_;
+  } else {
+    //safeguard in case the stride is set to 0 or 1
+    return getFullAtomList();
   }
-  reduced=true;
-  return requestlist_;
 }
 
 unsigned NeighborList::getStride() const {
@@ -373,4 +385,28 @@ std::vector<unsigned> NeighborList::getNeighbors(const unsigned index) const {
   return neighbors;
 }
 
+NeighborList::preparestatus NeighborList::prepare(Colvar* const aa,
+    bool firsttime,
+    bool invalidateList) {
+  if(stride_>0) {
+    if(stride_==1) {
+      invalidateList=true;
+      firsttime=false;
+    } else if(firsttime || (aa->getStep()%stride_==0 )) {
+      aa->requestAtoms(getFullAtomList());
+      invalidateList=true;
+      firsttime=false;
+    } else {
+      aa->requestAtoms(getReducedAtomList());
+      invalidateList=false;
+      if(aa->getExchangeStep()) {
+        aa->error("Neighbor lists should be updated on exchange steps - choose a NL_STRIDE which divides the exchange stride!");
+      }
+    }
+    if(aa->getExchangeStep()) {
+      firsttime=true;
+    }
+  }
+  return {firsttime,invalidateList};
+}
 } // namespace PLMD

src/tools/NeighborList.h

@@ -26,11 +26,13 @@
 #include "AtomNumber.h"
 
 #include <vector>
+#include <tuple>
 
 namespace PLMD {
 
 class Pbc;
 class Communicator;
+class Colvar;
 
 /// \ingroup TOOLBOX
 /// A class that implements neighbor lists from two lists or a single list of atoms
@@ -84,10 +86,14 @@ class NeighborList {
                bool doCells=false);
   ~NeighborList();
 /// Return the list of all atoms. These are needed to rebuild the neighbor list.
+/// Please use the `prepare()` method instead of directly calling this outside the constructor of your action
   std::vector<PLMD::AtomNumber>& getFullAtomList();
 /// Update the indexes in the neighbor list to match the
 /// ordering in the new positions array
 /// and return the new list of atoms that must be requested to the main code
+///
+/// Please prefer calling the `prepare()` method
+/// this is kept public for backward compatibility
   std::vector<PLMD::AtomNumber>& getReducedAtomList();
 /// Update the neighbor list and prepare the new
 /// list of atoms that will be requested to the main code
@@ -108,6 +114,18 @@ class NeighborList {
   std::vector<unsigned> getNeighbors(unsigned i) const;
 /// Get the i-th pair of AtomNumbers from the neighbor list
   pairAtomNumbers getClosePairAtomNumber(unsigned i) const;
+  struct preparestatus {
+    bool firsttime;
+    bool invalidateList;
+    // this little method makes possible to write
+    // std::tie(firsttime,invalidateList) = nl->prepare(this,firsttime).get();
+    // and keep the named nature of this "return struct"
+    std::tuple<bool,bool>get() const {
+      return {firsttime,invalidateList};
+    }
+  };
+/// Returns if the neighborlist is invalidated for this step
+  preparestatus prepare(Colvar*,bool firsttime, bool invalidateList);
 };
 
 } // namespace PLMD