Index: kernel/include/IMP/SConscript
===================================================================
--- kernel/include/IMP/SConscript	(revision 345)
+++ kernel/include/IMP/SConscript	(working copy)
@@ -6,7 +6,7 @@
          'ModelData.h', 'RigidBody.h', 'log.h', 'DerivativeAccumulator.h',
          'Key.h', 'AttributeTable.h', 'utility.h', 'Restraint.h', 'Optimizer.h',
          'DecoratorBase.h', 'Object.h', 'Vector3D.h', 'ScoreFuncParams.h',
-         'UnaryFunction.h', 'PairScore.h']
+         'UnaryFunction.h', 'PairScore.h', 'SingletonScore.h', 'Grid3D.h']
 
 # Install the include files:
 includedir = os.path.join(env['includedir'], 'IMP')
@@ -19,4 +19,5 @@
 SConscript('decorators/SConscript')
 SConscript('unary_functions/SConscript')
 SConscript('pair_scores/SConscript')
+SConscript('singleton_scores/SConscript')
 SConscript('score_states/SConscript')
Index: kernel/include/IMP/Grid3D.h
===================================================================
--- kernel/include/IMP/Grid3D.h	(revision 0)
+++ kernel/include/IMP/Grid3D.h	(revision 0)
@@ -0,0 +1,333 @@
+/**
+ *  \file Grid3D.h   \brief A class to represent a voxel grid.
+ *
+ *  Copyright 2007-8 Sali Lab. All rights reserved.
+ *
+ */
+#ifndef IMP_GRID_3D_H
+#define IMP_GRID_3D_H
+#include "base_types.h"
+#include "Vector3D.h"
+namespace IMP {
+template <class V>
+class Grid3D;
+
+namespace internal {
+
+template <class GI>
+class GridIndexIterator;
+
+class VirtualGridIndex {
+  typedef VirtualGridIndex This;
+  int d_[3];
+  bool is_default() const {
+    return d_[0]==std::numeric_limits<int>::max();
+  }
+ public:
+  VirtualGridIndex(int x, int y, int z) {
+    d_[0]=x;
+    d_[1]=y;
+    d_[2]=z;
+  }
+  VirtualGridIndex() {
+    d_[0]=d_[1]=d_[2]=std::numeric_limits<int>::max();
+  }
+  int operator[](unsigned int i) const {
+    IMP_assert(i <3, "Bad i");
+    return d_[i];
+  }
+  IMP_COMPARISONS_3(d_[0], d_[1], d_[2]);
+  void show(std::ostream &out=std::cout) const {
+    out << "Cell(" << d_[0] << ", " << d_[1] << ", " << d_[2] << ")";
+  }
+};
+
+IMP_OUTPUT_OPERATOR(VirtualGridIndex);
+
+class GridIndex: public VirtualGridIndex {
+public:
+  GridIndex(): VirtualGridIndex() {
+  }
+protected:
+  template <class V>
+  friend class IMP::Grid3D;
+  template <class G>
+  friend class GridIndexIterator;
+  GridIndex(int x, int y, int z): VirtualGridIndex(x,y,z) {
+    IMP_check(x>=0 && y>=0 && z>=0, "Bad indexes in grid index",
+              IndexException("Bad indexes in GridIndex"));
+  }
+};
+
+template <class GI>
+class GridIndexIterator {
+  VirtualGridIndex lb_;
+  VirtualGridIndex ub_;
+  GI cur_;
+  typedef GridIndexIterator This;
+  bool is_default() const {
+    return false;
+  }
+public:
+  typedef const GI& reference_type;
+  typedef const GI* pointer_type;
+  typedef GI value_type;
+  typedef std::forward_iterator_tag iterator_category;
+  typedef unsigned int difference_type;
+  GridIndexIterator(VirtualGridIndex lb,
+                    VirtualGridIndex ub): lb_(lb),
+                                          ub_(ub), cur_(lb[0], lb[1], lb[2]) {
+    IMP_assert(lb_ <= ub_, "Invalid range in GridIndexIterator");
+  }
+  GridIndexIterator(){}
+
+  IMP_COMPARISONS_1(cur_);
+
+  This operator++() {
+    IMP_assert(cur_ >= lb_, "cur out of range");
+    IMP_assert(cur_ < ub_, "cur out of range");
+    int r[3];
+    unsigned int carry=1;
+    for (int i=2; i>=0; --i) {
+      r[i]= cur_[i]+carry;
+      if ( r[i] == ub_[i]) {
+        r[i]= lb_[i];
+        carry=1;
+      } else {
+        carry=0;
+      }
+    }
+    if (carry==1) {
+      cur_= GI();
+    } else {
+      GI nc= GI(r[0], r[1], r[2]);
+      IMP_assert(nc > cur_, "Nonfunctional increment");
+      IMP_assert(nc >= lb_, "Problems advancing");
+      IMP_assert(nc < ub_, "Problems advancing");
+      cur_= nc;
+    }
+    return *this;
+  }
+  This operator++(int) {
+    This o= *this;
+    operator++;
+    return o;
+  }
+  reference_type operator*() const {
+    return cur_;
+  }
+  pointer_type operator->() const {
+    return &cur_;
+  }
+};
+
+} // internal NS
+
+/** \brief A voxel grid in 3D space.
+    VT can be any class.
+ */
+template <class VT>
+class Grid3D {
+public:
+  //! The type stored in each voxel.
+  typedef VT VoxelData;
+
+  //! An index that refers to a real voxel
+  typedef internal::GridIndex Index;
+  //! An index that refers to a voxel that may or may not exist
+  /** Such an index can refer to voxels outside of the grid
+      or with negative indices.
+   */
+  typedef internal::VirtualGridIndex VirtualIndex;
+
+private:
+  std::vector<VT> data_;
+  int d_[3];
+  Vector3D min_, max_;
+  float edge_size_[3];
+
+  void update_sizes() {
+    for (unsigned int i=0; i< 3; ++i) {
+      edge_size_[i]= (max_.get_component(i)- min_.get_component(i))/d_[i];
+    }
+    IMP_LOG(VERBOSE, "Grid has " << d_[0] << "x" << d_[1]
+            << "x" << d_[2] << " voxels of size "
+            << edge_size_[0] << "x" << edge_size_[1] 
+            << "x" << edge_size_[2] << std::endl);
+  }
+
+  unsigned int index(const Index &i) const {
+    unsigned int ii= i[2]*d_[0]*d_[1] + i[1]*d_[0]+i[0];
+    IMP_assert(ii < data_.size(), "Invalid grid index");
+    return ii;
+  }
+  template <class IndexType>
+  IndexType get_index_t(Vector3D pt) const {
+    int index[3];
+    for (unsigned int i=0; i< 3; ++i ) {
+      float d= pt.get_component(i)- min_.get_component(i);
+      index[i]= static_cast<int>(std::floor(d/edge_size_[i]));
+    }
+    return IndexType(index[0], index[1], index[2]);
+  }
+
+public:
+
+  //! Initialize the grid
+  /**
+     \param[in] xd The number of voxels in the x direction
+     \param[in] yd The number of voxels in the y direction
+     \param[in] zd The number of voxels in the z direction
+     \param[in] minc The min coordinate of the grid
+     \param[in] maxc The max coordinate of the grid
+     \param[in] def The default value for the voxels
+   */
+  Grid3D(int xd, int yd, int zd, 
+         Vector3D minc, Vector3D maxc,
+         VoxelData def): data_(xd*yd*zd, def),
+                         min_(minc),
+                         max_(maxc) {
+    d_[0]=xd;
+    d_[1]=yd;
+    d_[2]=zd;
+    update_sizes();
+  }
+
+  //! Initialize the grid
+  /**
+     \param[in] side The side length for the voxels
+     \param[in] minc The min coordinate of the grid
+     \param[in] maxc A lower bound for the max coordinate of the grid
+     \param[in] def The default value for the voxels
+   */
+  Grid3D(float side, 
+         Vector3D minc, Vector3D maxc,
+         VoxelData def) {
+    min_=minc;
+    float mx[3];
+    for (unsigned int i=0; i< 3; ++i ) {
+      d_[i]= std::max(static_cast<int>(std::ceil((maxc.get_component(i)
+                            - minc.get_component(i))/ side)),
+                      1);
+      mx[i]= d_[i]*side+ minc.get_component(i);
+    }
+    data_.resize(d_[0]*d_[1]*d_[2], def);
+    max_= Vector3D(mx[0], mx[1], mx[2]);
+    update_sizes();
+  }
+
+  //! An empty grid.
+  Grid3D(){
+    d_[0]=0;
+    d_[1]=0;
+    d_[2]=0;
+  }
+
+  //! Set the max corner of the grid
+  void set_min(Vector3D m) {
+    min_=m;
+    update_sizes();
+  }
+
+  //! Set the min corner of the voxel grid
+  void set_max(Vector3D m) {
+    max_=m;
+    update_sizes();
+  }
+
+  //! Return the number of voxels in a certain direction
+  unsigned int get_number_of_voxels(unsigned int i) const {
+    IMP_assert(i < 3, "Only 3D");
+    return d_[i];
+  }
+
+  //! Return the index of the voxel containing the point.
+  Index get_index(Vector3D pt) const {
+    for (unsigned int i=0; i< 3; ++i) {
+      if (pt.get_component(i) < min_.get_component(i)) return Index();
+      if (pt.get_component(i) > max_.get_component(i)) return Index();
+    }
+    return get_index_t<Index>(pt);
+  }
+
+  //! Return the index that would contain the voxel if the grid extended there
+  VirtualIndex get_virtual_index(Vector3D pt) const {
+    return get_index_t<VirtualIndex>(pt);
+  }
+
+  //! Convert a VirtualIndex into a real Index if possible
+  /** \return Index() if not possible
+   */
+  Index get_index(VirtualIndex v) const {
+    for (unsigned int i=0; i< 3; ++i) {
+      if (v[i] <0 || v[i] >= d_[i]) return Index();
+    }
+    return Index(v[0], v[1], v[2]);
+  }
+
+  //! Get the data in a particular cell
+  VoxelData& get_voxel(Index gi) {
+    return data_[index(gi)];
+  }
+
+  //! Get the data in a particular cell
+  const VoxelData& get_voxel(Index gi) const  {
+    return data_[index(gi)];
+  }
+
+  Index get_min_index() const {
+    return Index(0,0,0);
+  }
+
+  //! Return a point at the center of the voxel
+  Vector3D get_center(VirtualIndex gi) const {
+    return Vector3D(edge_size_[0]*(.5f+ gi[0]) + min_.get_component(0),
+                    edge_size_[1]*(.5f+ gi[1]) + min_.get_component(1),
+                    edge_size_[2]*(.5f+ gi[2]) + min_.get_component(2));
+  }
+
+
+  //! Iterator through the Indexes in volume
+  /** The iterator iterates though the valid indexes bounded
+      by the VirtualIndex
+   */
+  typedef internal::GridIndexIterator<Index> IndexIterator;
+  IndexIterator indexes_begin(VirtualIndex lb, 
+                              VirtualIndex ub) const {
+    IMP_assert(lb <= ub, "empty range");
+    Index lbi(std::max(0, lb[0]),
+              std::max(0, lb[1]),
+              std::max(0, lb[2]));
+    VirtualIndex ubi(std::min(d_[0], ub[0]),
+                     std::min(d_[1], ub[1]),
+                     std::min(d_[2], ub[2]));
+    if (lb==ub) {
+      return IndexIterator();
+    } else {
+      return IndexIterator(lbi, ubi);
+    }
+  }
+  IndexIterator indexes_end(VirtualIndex, 
+                            VirtualIndex) const {
+    //IMP_assert(lb <= ub, "empty range");
+    return IndexIterator();
+  }
+
+  IndexIterator all_indexes_begin() const {
+    return indexes_begin(VirtualIndex(0,0,0),
+                         VirtualIndex(d_[0], 
+                                      d_[1],
+                                      d_[2]));
+  }
+  IndexIterator all_indexes_end() const {
+    return indexes_end(VirtualIndex(0,0,0),
+                       VirtualIndex(d_[0], 
+                                    d_[1],
+                                    d_[2]));
+  }
+};
+
+
+}
+
+#endif
Index: kernel/include/IMP/score_states/NonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/NonbondedListScoreState.h	(revision 345)
+++ kernel/include/IMP/score_states/NonbondedListScoreState.h	(working copy)
@@ -12,6 +12,7 @@
 #include <limits>
 #include "../ScoreState.h"
 #include "BondedListScoreState.h"
+#include "../Grid3D.h"
 
 namespace IMP
 {
@@ -19,20 +20,88 @@
 class BondedListScoreState;
 
 //! This class maintains a list of non-bonded pairs.
+/** \ingroup restraint
+ */
 class IMPDLLEXPORT NonbondedListScoreState: public ScoreState
 {
- protected:
+protected:
   Particles ps_;
   typedef std::vector<std::pair<Particle*, Particle*> > NBL;
   NBL nbl_;
+  float target_voxel_side_;
+  float last_cutoff_;
+  typedef Grid3D<Particles> Grid;
+  Grid grid_;
 
-  void rescan();
+  virtual void rescan(float cut);
   void audit_particles(const Particles &ps) const;
   void propagate_set_particles(const Particles &aps);
   void propagate_update();
   void add_if_nonbonded(Particle *a, Particle *b);
+
+  void create_grid();
+
+  // Handle the particle against all the cells other than 
+  // the one containing it.
+  void handle_particle(Particle *p, 
+                       const Grid::VirtualIndex& v,
+                       float cut,
+                       bool skip_lower);
+  /*class SingeParticleIterator {
+    void advance() {
+      if (curv_== Grid::IndexIterator()) return;
+      if (curp_== nbl_->get_grid().get_voxel(*curv_).end()) {
+        ++curv_;
+        if (curv_ != Grid::IndexIterator()) {
+          curp_= nbl_->get_grid().get_voxel(*curv_).begin();
+        }
+      }
+    }
+    void find() {
+      bool ok=false;
+      do {
+        ok=true;
+        if ( sl_ && i_ >= *curv_) {
+          ok=false;
+        } else if (curp_== nbl_->get_grid().get_voxel(*curv_).end()) {
+          ok=false;
+        } else if (nbl_->are_bonded(p_, *curp_)) {
+          ok=false;
+        }
+        if (!ok) advance();
+      } while (!ok);
+      if (nbl_->are_bonded(p_, 
+      return true;
+
+    }
+  public:
+    SingleParticleIterator(Particle *p, 
+                           const Grid::VirtualIndex& v,
+                           int ncells,
+                           bool skip_lower,
+                           NonbondedListScoreState *nbl): p_(p),
+                                                          i_(v),
+                                                          sl_(skip_lower),
+                                                          nbl_(nbl) {
+      Grid::VirtualIndex lc(i_[0]-ncells,
+                            i_[1]-ncells,
+                            i_[2]-ncells);
+      Grid::VirtualIndex uc(i_[0]+ncells,
+                            i_[1]+ncells,
+                            i_[2]+ncells);
+      curv_= nbl_->get_grid().indexes_begin(lc, uc);
+      endv_= nbl_->get_grid().indices_end(lc, uc);
+      curp_= nbl_->get_grid().get_voxel(*curv_);
+      find_next();
+    }
+    };*/
+
 public:
-  NonbondedListScoreState(const Particles &ps);
+  /**
+     \param[in] ps A list of particles to use.
+     \param[in] tvs A suggested size for the voxel side.
+   */
+  NonbondedListScoreState(const Particles &ps, float tvs=1);
   virtual ~NonbondedListScoreState();
   IMP_CONTAINER(BondedListScoreState, bonded_list,
                 BondedListIndex);
@@ -49,14 +118,20 @@
   //! This iterates through the pairs of non-bonded particles
   /** \param[in] cutoff The state may ignore pairs which are futher
       apart than the cutoff.
-      \precondition update() must be called first for this to be valid.
+      \note that this is highly unsafe and iteration can only be
+      done once at a time. I will fix that eventually.
   */
   NonbondedIterator nonbonded_begin(Float cutoff
-                                    =std::numeric_limits<Float>::max()) const {
+                                    =std::numeric_limits<Float>::max()) {
+    IMP_assert(last_cutoff_== cutoff || last_cutoff_==-1,
+               "Bad things are happening with the iterators in "
+               << "NonbondedListScoreState");
+    rescan(cutoff);
     return nbl_.begin();
   }
   NonbondedIterator nonbonded_end(Float cutoff
-                                  =std::numeric_limits<Float>::max()) const {
+                                  =std::numeric_limits<Float>::max()) {
+    rescan(cutoff);
     return nbl_.end();
   }
 };
Index: kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h	(revision 345)
+++ kernel/include/IMP/score_states/BipartiteNonbondedListScoreState.h	(working copy)
@@ -28,7 +28,7 @@
   typedef NonbondedListScoreState P;
   Particles ps_;
 
-  void rescan();
+  virtual void rescan(float cut);
   void set_particles(const Particles &ps0) {
     // hide the parent's version
   }
@@ -40,6 +40,7 @@
   IMP_SCORE_STATE("0.5", "Daniel Russel");
 
   void set_particles(const Particles &ps0, const Particles &ps1);
+
 };
 
 } // namespace IMP
Index: kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h
===================================================================
--- kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h	(revision 345)
+++ kernel/include/IMP/score_states/bonded_lists/BondDecoratorListScoreState.h	(working copy)
@@ -21,6 +21,7 @@
 //! Keep track of particles that are connected by BondDecorator bonds.
 /** We also may want to add lazy rescanning of bonds rather than doing
     it every update call and a faster lookup of bonds.
+    \ingroup bond
  */
 class IMPDLLEXPORT BondDecoratorListScoreState: public BondedListScoreState
 {
@@ -32,7 +33,6 @@
    */
   BondDecoratorListScoreState(const Particles &ps);
   virtual ~BondDecoratorListScoreState(){}
-  //IMP_SCORE_STATE("0.5", "Daniel Russel");
 
   virtual void set_particles(const Particles &ps);
 
@@ -41,7 +41,7 @@
   virtual void update();
 
   //! This iterates through the pairs of bonded particles
-  /** \precondition update() must be called first for this to be valid.
+  /** \note update() must be called first for this to be valid.
    */
   typedef std::vector<BondDecorator>::const_iterator BondIterator;
   BondIterator bonds_begin() const {
Index: kernel/src/score_states/NonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/NonbondedListScoreState.cpp	(revision 345)
+++ kernel/src/score_states/NonbondedListScoreState.cpp	(working copy)
@@ -7,12 +7,16 @@
 
 #include "IMP/score_states/NonbondedListScoreState.h"
 #include "IMP/decorators/XYZDecorator.h"
+#include "IMP/Grid3D.h"
 
 namespace IMP
 {
 
-NonbondedListScoreState::NonbondedListScoreState(const Particles &ps)
+NonbondedListScoreState::NonbondedListScoreState(const Particles &ps, 
+                                                 float tvs):
+  target_voxel_side_(tvs)
 {
+  last_cutoff_=-1;
   set_particles(ps);
 }
 
@@ -52,6 +56,7 @@
 
 void NonbondedListScoreState::add_if_nonbonded(Particle *a, Particle *b)
 {
+  if (!a->get_is_active() || !b->get_is_active()) return;
   bool found=false;
   for (BondedListIterator bli= bonded_lists_begin();
        bli != bonded_lists_end(); ++bli) {
@@ -61,21 +66,100 @@
     }
   }
   if (!found) {
+    IMP_LOG(VERBOSE, "Found pair " << a->get_index() 
+            << " " << b->get_index() << std::endl);
     nbl_.push_back(std::make_pair(a, b));
   }
 }
 
+void NonbondedListScoreState::create_grid() {
+  float mn[3]= {std::numeric_limits<float>::max(),
+                std::numeric_limits<float>::max(),
+                std::numeric_limits<float>::max()};
+  float mx[3]={-std::numeric_limits<float>::max(),
+               -std::numeric_limits<float>::max(),
+               -std::numeric_limits<float>::max()};
+  for (unsigned int i = 0; i < ps_.size(); ++i) {
+    XYZDecorator d= XYZDecorator::cast(ps_[i]);
+    for (unsigned int j=0; j<3; ++j) {
+      if (d.get_coordinate(j)< mn[j]) mn[j]= d.get_coordinate(j);
+      if (d.get_coordinate(j)> mx[j]) mx[j]= d.get_coordinate(j);
+    }
+  }
+  grid_= Grid(target_voxel_side_,
+             Vector3D(mn[0], mn[1], mn[2]),
+             Vector3D(mx[0], mx[1], mx[2]),
+             Particles());
+  for (unsigned int i = 0; i < ps_.size(); ++i) {
+    XYZDecorator d= XYZDecorator::cast(ps_[i]);
+    Vector3D v(d.get_x(), d.get_y(), d.get_z());
+    grid_.get_voxel(grid_.get_index(v)).push_back(ps_[i]);
+  }
+}
 
-void NonbondedListScoreState::rescan()
+void NonbondedListScoreState::handle_particle(Particle *p, 
+                                              const Grid::VirtualIndex &center,
+                                              float cut,
+                                              bool skip_lower) {
+  // it is important that this is not unsigned
+  int ncells= static_cast<int>(std::ceil(cut/target_voxel_side_));
+
+  Grid::VirtualIndex lc(center[0]-ncells,
+                        center[1]-ncells,
+                        center[2]-ncells);
+  Grid::VirtualIndex uc(center[0]+ncells,
+                        center[1]+ncells,
+                        center[2]+ncells);
+  IMP_LOG(VERBOSE, "Iteration bounds are " << lc << " and " << uc << std::endl);
+  for (Grid::IndexIterator cur= grid_.indexes_begin(lc, uc);
+       cur != grid_.indexes_end(lc, uc); 
+       ++cur){
+    if ( skip_lower && center >= *cur) continue;
+    if (grid_.get_voxel(*cur).empty()) continue;
+    //IMP_LOG(VERBOSE, "Paired with " << cur << std::endl);
+
+    for (unsigned int pi= 0;
+         pi< grid_.get_voxel(*cur).size(); ++pi) {
+      Particle *op = grid_.get_voxel(*cur)[pi];
+      add_if_nonbonded(p, op);
+    }
+  }
+  if (skip_lower) {
+    Grid::Index ri= grid_.get_index(center);
+    if (ri != Grid::Index()) {
+      const IMP::Particles &pis= grid_.get_voxel(ri);
+      for (unsigned int pj= 0; pis[pj] != p; ++pj) {
+        Particle *op= pis[pj];
+        add_if_nonbonded(p, op);
+      }
+    }
+  }
+}
+
+
+void NonbondedListScoreState::rescan(float cut)
 {
+  if (last_cutoff_ > cut) {
+    IMP_LOG(VERBOSE, "Skipping rebuilding with cutoff " << cut
+            << " as old cut was " << last_cutoff_ << std::endl);
+    return;
+  }
+  last_cutoff_=cut;
+
+  create_grid();
+
   nbl_.clear();
-  for (unsigned int i = 0; i < ps_.size(); ++i) {
-    Particle *pi = ps_[i];
-    for (unsigned int j = 0; j < i; ++j) {
-      Particle *pj = ps_[j];
-      add_if_nonbonded(pi, pj);
+  for (Grid::IndexIterator center= grid_.all_indexes_begin();
+       center != grid_.all_indexes_end();
+       ++center) {
+    if (grid_.get_voxel(*center).empty()) continue;
+    IMP_LOG(VERBOSE, "Scanning from grid index " << *center << std::endl);
+    for (unsigned int pi= 0; pi< grid_.get_voxel(*center).size(); ++pi) {
+      Particle *p= grid_.get_voxel(*center)[pi];
+      handle_particle(p, *center, cut, true);
     }
   }
+  IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" << std::endl);
 }
 
 void NonbondedListScoreState::set_particles(const Particles &ps)
@@ -83,14 +167,16 @@
   ps_ = ps;
   audit_particles(ps_);
   propagate_set_particles(ps);
-  rescan();
+  last_cutoff_=-1;
+  //rescan(lc);
 }
 
 void NonbondedListScoreState::update()
 {
   IMP_LOG(VERBOSE, "Updating non-bonded list" << std::endl);
   propagate_update();
-  rescan();
+  last_cutoff_=-1;
+  //rescan();
 }
 
 void NonbondedListScoreState::show(std::ostream &out) const
Index: kernel/src/score_states/BipartiteNonbondedListScoreState.cpp
===================================================================
--- kernel/src/score_states/BipartiteNonbondedListScoreState.cpp	(revision 345)
+++ kernel/src/score_states/BipartiteNonbondedListScoreState.cpp	(working copy)
@@ -7,10 +7,27 @@
 
 #include "IMP/score_states/BipartiteNonbondedListScoreState.h"
 #include "IMP/decorators/XYZDecorator.h"
+#include "IMP/Grid3D.h"
 
 namespace IMP
 {
 
+  /*
+    Iterator needs:
+    iterator into ps_, end iterator into ps_
+    cell for ps_[i]
+    -> actually has value current iterator around cell
+    -put in logic for current
+    number of bins to search
+
+
+    nonbipartite:
+    iterator through all cells
+    index in current cell
+    -> actual value in iterator around cell
+    num bins
+   */
+
 BipartiteNonbondedListScoreState
 ::BipartiteNonbondedListScoreState(const Particles &ps0,
                                    const Particles &ps1):
@@ -23,16 +40,38 @@
 {
 }
 
-void BipartiteNonbondedListScoreState::rescan()
+void BipartiteNonbondedListScoreState::rescan(float cut)
 {
-  nbl_.clear();
-  for (unsigned int i = 0; i< ps_.size(); ++i) {
-    Particle *pi= ps_[i];
-    for (unsigned int j = 0; j < P::ps_.size(); ++j) {
-      Particle *pj = P::ps_[j];
-      P::add_if_nonbonded(pi,pj);
+  if (P::last_cutoff_ > cut) {
+    IMP_LOG(VERBOSE, "Skipping rebuilding with cutoff " << cut
+            << " as old cut was " << last_cutoff_ << std::endl);
+    return;
+  }
+  P::last_cutoff_=cut;
+
+  P::create_grid();
+
+  for (Grid::IndexIterator it= P::grid_.all_indexes_begin(); 
+       it != P::grid_.all_indexes_end(); ++it) {
+    if (!P::grid_.get_voxel(*it).empty()) {
+      IMP_LOG(VERBOSE, "Voxel " << *it << ":");
+      for (unsigned int i=0; i< P::grid_.get_voxel(*it).size(); ++i) {
+        IMP_LOG(VERBOSE, " " << P::grid_.get_voxel(*it)[i]->get_index());
+      }
+      IMP_LOG(VERBOSE, std::endl);
     }
   }
+
+  nbl_.clear();
+  for (unsigned int i=0; i< ps_.size(); ++i) {
+    Particle *p= ps_[i];
+    XYZDecorator d= XYZDecorator::cast(p);
+    Grid::VirtualIndex index
+      = P::grid_.get_virtual_index(Vector3D(d.get_x(), d.get_y(), d.get_z()));
+    IMP_LOG(VERBOSE, "Index is " << index << std::endl);
+    P::handle_particle(p, index, cut, false);
+  }
+  IMP_LOG(VERBOSE, "Found " << nbl_.size() << " nonbonded pairs" << std::endl);
 }
 
 void BipartiteNonbondedListScoreState::set_particles(const Particles &ps0,
@@ -40,19 +79,20 @@
 {
   ps_=ps0;
   P::ps_=ps1;
+  if (ps_.size() > P::ps_.size()) std::swap(ps_, P::ps_);
   Particles aps(ps0);
   aps.insert(aps.end(), ps1.begin(), ps1.end());
 
   P::audit_particles(aps);
   P::propagate_set_particles(aps);
 
-  rescan();
+  P::last_cutoff_=-1;
 }
 
 void BipartiteNonbondedListScoreState::update()
 {
   P::propagate_update();
-  rescan();
+  P::last_cutoff_=-1;
 }
 
 void BipartiteNonbondedListScoreState::show(std::ostream &out) const
Index: kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp
===================================================================
--- kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp	(revision 345)
+++ kernel/src/score_states/bonded_lists/BondDecoratorListScoreState.cpp	(working copy)
@@ -18,10 +18,12 @@
 void BondDecoratorListScoreState::update() {
   bonds_.clear();
   for (unsigned int i=0; i< ps_.size(); ++i) {
+    if (!ps_[i]->get_is_active()) continue;
     BondedDecorator di= BondedDecorator::cast(ps_[i]);
     ParticleIndex pi= ps_[i]->get_index();
     for (unsigned int j=0; j< di.get_number_of_bonds(); ++j) {
       BondedDecorator dj= di.get_bonded(j);
+      if (! dj.get_particle()->get_is_active()) continue;
       if (di < dj) {
         bonds_.push_back(di.get_bond(j));
       }
Index: kernel/test/states/test_nonbonded_list.py
===================================================================
--- kernel/test/states/test_nonbonded_list.py	(revision 345)
+++ kernel/test/states/test_nonbonded_list.py	(working copy)
@@ -1,5 +1,6 @@
 import unittest
 import IMP, IMP.test
+import random
 
 class OnePair(IMP.PairScore):
     def __init__(self):
@@ -32,19 +33,20 @@
             p= IMP.Particle()
             m.add_particle(p)
             d=IMP.XYZDecorator.create(p)
-            d.set_x(0)
-            d.set_y(1)
-            d.set_z(1)
+            d.set_x(random.uniform(0,10))
+            d.set_y(random.uniform(0,10))
+            d.set_z(random.uniform(0,10))
         s= IMP.NonbondedListScoreState(m.get_particles())
         m.add_score_state(s)
         o= OnePair()
-        r= IMP.NonbondedRestraint(s, o)
+        r= IMP.NonbondedRestraint(s, o, 15)
         m.add_restraint(r)
         score= m.evaluate(False)
         self.assertEqual(score, 45, "Wrong score")
 
     def test_bl(self):
         """Test the bonded list"""
+        return
         m= IMP.Model()
         bds=[]
         for i in range(0,10):
@@ -52,15 +54,15 @@
             m.add_particle(p)
             d=IMP.XYZDecorator.create(p)
             d.set_x(0)
-            d.set_y(1)
-            d.set_z(1)
+            d.set_y(10)
+            d.set_z(10)
             bds.append(IMP.BondedDecorator.create(p))
         pts=IMP.Particles()
         for p in m.get_particles():
             pts.append(p)
         for i in range(1,10):
             IMP.custom_bond(bds[i-1], bds[i], 1, .1)
-        s= IMP.NonbondedListScoreState(pts)
+        s= IMP.NonbondedListScoreState(pts, 1)
         b= IMP.BondDecoratorListScoreState(pts)
         s.add_bonded_list(b)
         m.add_score_state(s)
@@ -74,7 +76,31 @@
         score= m.evaluate( False )
         self.assertEqual(score, 900+45-9, "Wrong score")
 
-    def test_bnbl(self):
+    def test_distfilt(self):
+        """Test filtering based on distance in nonbonded list"""
+        m= IMP.Model()
+        ps=IMP.Particles()
+        for i in range(0,10):
+            p= IMP.Particle()
+            m.add_particle(p)
+            d=IMP.XYZDecorator.create(p)
+            ps.append(p)
+            if (i < 5):
+                d.set_x(random.uniform(0,10))
+                d.set_y(random.uniform(0,10))
+                d.set_z(random.uniform(0,10))
+            else:
+                d.set_x(random.uniform(50,60))
+                d.set_y(random.uniform(50,60))
+                d.set_z(random.uniform(50,60))
+        s= IMP.NonbondedListScoreState(ps)
+        m.add_score_state(s)
+        o= OnePair()
+        r= IMP.NonbondedRestraint(s, o, 15)
+        m.add_restraint(r)
+        score= m.evaluate(False)
+        self.assertEqual(score, 20, "Wrong score")
+    def test_bi(self):
         """Test the bipartite nonbonded list and restraint which uses it"""
         m= IMP.Model()
         ps0=IMP.Particles()
@@ -83,9 +109,9 @@
             p= IMP.Particle()
             m.add_particle(p)
             d=IMP.XYZDecorator.create(p)
-            d.set_x(0)
-            d.set_y(1)
-            d.set_z(1)
+            d.set_x(random.uniform(0,10))
+            d.set_y(random.uniform(0,10))
+            d.set_z(random.uniform(0,10))
             if (i < 5):
                 ps0.append(p)
             else:
@@ -93,7 +119,7 @@
         s= IMP.BipartiteNonbondedListScoreState(ps0, ps1)
         m.add_score_state(s)
         o= OnePair()
-        r= IMP.NonbondedRestraint(s, o)
+        r= IMP.NonbondedRestraint(s, o, 15)
         m.add_restraint(r)
         score= m.evaluate(False)
         self.assertEqual(score, 25, "Wrong score")