diff --git a/package/CHANGELOG b/package/CHANGELOG
index bc64f31b466..dacd361d753 100644
--- a/package/CHANGELOG
+++ b/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
diff --git a/package/MDAnalysis/analysis/rdf.py b/package/MDAnalysis/analysis/rdf.py
index 504c877b511..0d96b7ac106 100644
--- a/package/MDAnalysis/analysis/rdf.py
+++ b/package/MDAnalysis/analysis/rdf.py
@@ -235,6 +235,7 @@ def _conclude(self):
 
 class InterRDF_s(AnalysisBase):
     """Site-specific intermolecular pair distribution function
+
     Arguments
     ---------
     u : Universe
@@ -253,32 +254,48 @@ class InterRDF_s(AnalysisBase):
     step : int (optional)
           The step size through the trajectory in frames (default is
           every frame)
+
     Example
     -------
+
     First create the :class:`InterRDF_s` object, by supplying one Universe and
     one list of pairs of AtomGroups, then use the :meth:`~InterRDF_s.run`
     method::
+
       from MDAnalysisTests.datafiles import GRO_MEMPROT, XTC_MEMPROT
       u = mda.Universe(GRO_MEMPROT, XTC_MEMPROT)
+
       s1 = u.select_atoms('name ZND and resid 289')
       s2 = u.select_atoms('(name OD1 or name OD2) and resid 51 and sphzone 5.0 (resid 289)')
       s3 = u.select_atoms('name ZND and (resid 291 or resid 292)')
       s4 = u.select_atoms('(name OD1 or name OD2) and sphzone 5.0 (resid 291)')
       ags = [[s1, s2], [s3, s4]]
+
       rdf = InterRDF_s(u, ags)
       rdf.run()
+
     Results are available through the :attr:`bins` and :attr:`rdf` attributes::
+
       plt.plot(rdf.bins, rdf.rdf[0][0][0])
+
     (Which plots the rdf between the first atom in ``s1`` and the first atom in
     ``s2``)
+
     To generate the *cumulative distribution function* (cdf), use the
     :meth:`~InterRDF_s.get_cdf` method ::
+
       cdf = rdf.get_cdf()
+
     Results are available through the :attr:'cdf' attribute::
+
       plt.plot(rdf.bins, rdf.cdf[0][0][0])
+
     (Which plots the cdf between the first atom in ``s1`` and the first atom in
     ``s2``)
+
+
     .. versionadded:: 0.19.0
+
     """
     def __init__(self, u, ags,
                  nbins=75, range=(0.0, 15.0), density=True, **kwargs):