Modified the distance function in RDF (#2013)

* using capped_distance function in rdf

* updated docs, added cutoff in determine methods

* extended RDF benchmarks to cover range of atoms

benchmarks/benchmarks/analysis/rdf.py

@@ -17,21 +17,23 @@ class SimpleRdfBench(object):
     """
 
     params = ([20,75,200],
-              [[0,5], [0,15], [0,20]])
+              [[0,5], [0,15], [0,20]],
+              [1, 100, 1000, 10000])
 
     param_names = ['nbins',
-                   'range_val']
+                   'range_val',
+                   'natoms']
 
-    def setup(self, nbins, range_val):
+    def setup(self, nbins, range_val, natoms):
 
         self.sel_str = 'name OW'
 
         self.u = MDAnalysis.Universe(TPR, XTC)
 
         try:
-            self.sel = self.u.select_atoms(self.sel_str)[:200]
+            self.sel = self.u.select_atoms(self.sel_str)[:natoms]
         except AttributeError:
-            self.sel = self.u.selectAtoms(self.sel_str)[:200]
+            self.sel = self.u.selectAtoms(self.sel_str)[:natoms]
 
         # do not include initialization of the
         # InterRDF object in the benchmark itself

package/CHANGELOG

@@ -21,10 +21,13 @@ The rules for this file:
 
 Enhancements
 
+
   * Replaced multiple apply (_apply_distmat, _apply_kdtree)
     methods in distance based selections with
     lib.distances.capped_distance for automatic selection of
     distance evaluation method. (PR #2035)
+  * Modified analysis.rdf.InterRDF to use lib.distances.capped_distance
+    to handle rdf calculations for large systems. (PR #2013)
   * Added return_distances argument in lib.distances.capped_distances
     to evaluate and return distances only when required. Modified
     the optimization rules in lib.distances._determine_method for

package/MDAnalysis/analysis/rdf.py

@@ -187,26 +187,23 @@ def _prepare(self):
 
         # Need to know average volume
         self.volume = 0.0
+        # Set the max range to filter the search radius
+        self._maxrange = self.rdf_settings['range'][1]
 
-        # Allocate a results array which we will reuse
-        self._result = np.zeros((len(self.g1), len(self.g2)), dtype=np.float64)
-        # If provided exclusions, create a mask of _result which
-        # lets us take these out
-        if self._exclusion_block is not None:
-            self._exclusion_mask = blocks_of(self._result,
-                                             *self._exclusion_block)
-            self._maxrange = self.rdf_settings['range'][1] + 1.0
-        else:
-            self._exclusion_mask = None
 
     def _single_frame(self):
-        distances.distance_array(self.g1.positions, self.g2.positions,
-                                 box=self.u.dimensions, result=self._result)
+        pairs, dist = distances.capped_distance(self.g1.positions,
+                                                self.g2.positions,
+                                                self._maxrange,
+                                                box=self.u.dimensions)
         # Maybe exclude same molecule distances
-        if self._exclusion_mask is not None:
-            self._exclusion_mask[:] = self._maxrange
+        if self._exclusion_block is not None:
+            idxA, idxB = pairs[:, 0]//self._exclusion_block[0], pairs[:, 1]//self._exclusion_block[1]
+            mask = np.where(idxA != idxB)[0]
+            dist = dist[mask]
+
 
-        count = np.histogram(self._result, **self.rdf_settings)[0]
+        count = np.histogram(dist, **self.rdf_settings)[0]
         self.count += count
 
         self.volume += self._ts.volume
@@ -370,15 +367,12 @@ def _conclude(self):
 
     def get_cdf(self):
         """Calculate the cumulative distribution functions (CDF) for all sites.
-
         Note that this is the actual count within a given radius, i.e.,
         :math:`N(r)`.
-
         Returns
         -------
               cdf : list
                       list of arrays with the same structure as :attr:`rdf`
-
         """
         # Calculate cumulative distribution function
         # Empty list to restore CDF
@@ -391,4 +385,4 @@ def get_cdf(self):
         # self.cdf is a list of cdf between pairs of AtomGroups in ags
         self.cdf = cdf
 
-        return cdf
+        return cdf