diff --git a/package/CHANGELOG b/package/CHANGELOG
index 83bff74002..cdefee0ac9 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -14,7 +14,8 @@ The rules for this file:
 
 
 -------------------------------------------------------------------------------
-??/??/?? IAlibay, ChiahsinChu, RMeli
+??/??/?? IAlibay, ChiahsinChu, RMeli, talagayev
+
 
  * 2.9.0
 
@@ -23,6 +24,8 @@ Fixes
 Enhancements
  * Add check and warning for empty (all zero) coordinates in RDKit converter (PR #4824)
  * Added `precision` for XYZWriter (Issue #4775, PR #4771)
+ * Enables parallelization of `analysis.atomicdistances.AtomicDistances`
+   (Issue #4819)
 
 Changes
 
diff --git a/package/MDAnalysis/analysis/atomicdistances.py b/package/MDAnalysis/analysis/atomicdistances.py
index 1860d3285b..08ab33f71b 100644
--- a/package/MDAnalysis/analysis/atomicdistances.py
+++ b/package/MDAnalysis/analysis/atomicdistances.py
@@ -114,6 +114,7 @@
 
 import logging
 from .base import AnalysisBase
+from .results import Results
 
 logger = logging.getLogger("MDAnalysis.analysis.atomicdistances")
 
@@ -145,6 +146,10 @@ class AtomicDistances(AnalysisBase):
 
 
     .. versionadded:: 2.5.0
+
+    .. versionchanged:: 2.9.0
+        Enables parallelization through
+        the use of `self.results.distances`.
     """
 
     def __init__(self, ag1, ag2, pbc=True, **kwargs):
@@ -163,14 +168,21 @@ def __init__(self, ag1, ag2, pbc=True, **kwargs):
         self._ag1 = ag1
         self._ag2 = ag2
         self._pbc = pbc
+        self.results = Results()
 
     def _prepare(self):
         # initialize NumPy array of frames x distances for results
-        self.results = np.zeros((self.n_frames, self._ag1.atoms.n_atoms))
+        self.results.distances = np.zeros((
+            self.n_frames, self._ag1.atoms.n_atoms
+        ))
 
     def _single_frame(self):
         # if PBCs considered, get box size
         box = self._ag1.dimensions if self._pbc else None
-        self.results[self._frame_index] = calc_bonds(self._ag1.positions,
-                                                     self._ag2.positions,
-                                                     box)
+        self.results.distances[self._frame_index] = calc_bonds(
+            self._ag1.positions, self._ag2.positions, box
+        )
+
+    def _conclude(self):
+        # adjust self.results to self.results.distances
+        self.results = self.results.distances
diff --git a/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py b/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py
index 443173cff7..38dead5f24 100644
--- a/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py
+++ b/testsuite/MDAnalysisTests/analysis/test_atomicdistances.py
@@ -118,8 +118,9 @@ def test_ad_pairwise_dist(self, ad_ag1, ad_ag2,
                               expected_dist):
         '''Ensure that pairwise distances between atoms are
         correctly calculated without PBCs.'''
-        pairwise_no_pbc = (ad.AtomicDistances(ad_ag1, ad_ag2,
-                                              pbc=False).run())
+        pairwise_no_pbc = ad.AtomicDistances(
+            ad_ag1, ad_ag2, pbc=False
+        ).run()
         actual = pairwise_no_pbc.results
 
         # compare with expected values from dist()
@@ -129,7 +130,7 @@ def test_ad_pairwise_dist_pbc(self, ad_ag1, ad_ag2,
                                   expected_pbc_dist):
         '''Ensure that pairwise distances between atoms are
         correctly calculated with PBCs.'''
-        pairwise_pbc = (ad.AtomicDistances(ad_ag1, ad_ag2).run())
+        pairwise_pbc = ad.AtomicDistances(ad_ag1, ad_ag2).run()
         actual = pairwise_pbc.results
 
         # compare with expected values from dist()