diff --git a/package/AUTHORS b/package/AUTHORS
index 028e3a3ef25..245e6bd7345 100644
--- a/package/AUTHORS
+++ b/package/AUTHORS
@@ -213,6 +213,7 @@ Chronological list of authors
   - Ahmed Salah Ghoneim
   - Alexander Schlaich
   - Josh Vermaas
+  - Xiaoxu Ruan
 
 External code
 -------------
diff --git a/package/CHANGELOG b/package/CHANGELOG
index 480bb224a6a..2d489056e09 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -16,7 +16,7 @@ The rules for this file:
   
 ??/??/?? IAlibay, pgbarletta, mglagolev, hmacdope, manuel.nuno.melo, chrispfae, 
          ooprathamm, MeetB7, BFedder, v-parmar, MoSchaeffler, jbarnoud, jandom,
-         xhgchen, jaclark5, DrDomenicoMarson, AHMED-salah00, schlaicha
+         xhgchen, jaclark5, DrDomenicoMarson, AHMED-salah00, schlaicha, SophiaRuan
  * 2.5.0
 
 Fixes
@@ -49,6 +49,8 @@ Enhancements
     and SegmentGroup. (PR #3953)
 
 Changes
+  * Add progress bars to track the progress of _conclude() functions
+    (_conclude_simple() and _conclude_fft()) in msd.py (Issue #4070, PR #4072)
   * As per NEP29 the minimum supported NumPy version has been raised to 1.21
     (note: in practice later versions of NumPy may be used depending on your
      architecture, operating system, or Python version) (PR #3983)
diff --git a/package/MDAnalysis/analysis/msd.py b/package/MDAnalysis/analysis/msd.py
index 6c3104a9ce6..c3d202e883a 100644
--- a/package/MDAnalysis/analysis/msd.py
+++ b/package/MDAnalysis/analysis/msd.py
@@ -250,6 +250,7 @@
 from ..due import due, Doi
 from .base import AnalysisBase
 from ..core import groups
+from tqdm import tqdm
 
 logger = logging.getLogger('MDAnalysis.analysis.msd')
 
@@ -387,7 +388,7 @@ def _conclude_simple(self):
         """
         lagtimes = np.arange(1, self.n_frames)
         positions = self._position_array.astype(np.float64)
-        for lag in lagtimes:
+        for lag in tqdm(lagtimes):
             disp = positions[:-lag, :, :] - positions[lag:, :, :]
             sqdist = np.square(disp).sum(axis=-1)
             self.results.msds_by_particle[lag, :] = np.mean(sqdist, axis=0)
@@ -411,7 +412,7 @@ def _conclude_fft(self):  # with FFT, np.float64 bit prescision required.
                 or set fft=False""")
 
         positions = self._position_array.astype(np.float64)
-        for n in range(self.n_particles):
+        for n in tqdm(range(self.n_particles)):
             self.results.msds_by_particle[:, n] = tidynamics.msd(
                 positions[:, n, :])
         self.results.timeseries = self.results.msds_by_particle.mean(axis=1)