diff --git a/CHANGELOG b/CHANGELOG
index 89d9b03d..a9bc8829 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -14,7 +14,7 @@ The rules for this file:
 
 ------------------------------------------------------------------------------
 =======
-MM/DD/YYYY VOD555, lilyminium, orbeckst
+MM/DD/YYYY VOD555, lilyminium, orbeckst, yuxuanzhuang
 
   * 0.4.0
 
@@ -29,6 +29,7 @@ Fixes
     - threaded --> threads
   * fixed RDF functions (gave wrong results if step != 1) (#114)
   * fixed InterRDF_s function (gave wrong results if density=True) (#120)
+  * fixed Contact fails with uneven blocks. (#140)
 
 Changes
   * requires MDAnalysis >= 1.0.0 (#122)
diff --git a/pmda/contacts.py b/pmda/contacts.py
index 419afbef..510cc0e9 100644
--- a/pmda/contacts.py
+++ b/pmda/contacts.py
@@ -286,7 +286,7 @@ def _single_frame(self, ts, atomgroups):
         return y
 
     def _conclude(self):
-        self.timeseries = np.hstack(self._results)
+        self.timeseries = np.concatenate(self._results)
 
 
 def q1q2(atomgroup, radius=4.5):
diff --git a/pmda/parallel.py b/pmda/parallel.py
index 636c7d5b..361edc4d 100644
--- a/pmda/parallel.py
+++ b/pmda/parallel.py
@@ -153,7 +153,7 @@ def _conclude(self):
                # for each frame are stored in ``self._results`` in a per block
                # basis. Here those results should be moved and reshaped into a
                # sensible new variable.
-               self.results = np.hstack(self._results)
+               self.results = np.concatenate(self._results)
                # Apply normalisation and averaging to results here if wanted.
                self.results /= np.sum(self.results
 
diff --git a/pmda/test/test_contacts.py b/pmda/test/test_contacts.py
index f2fbc896..c201b6dc 100644
--- a/pmda/test/test_contacts.py
+++ b/pmda/test/test_contacts.py
@@ -63,12 +63,13 @@ def _run_Contacts(self,
                       start=None,
                       step=None,
                       stop=None,
+                      n_blocks=1,
                       **kwargs):
         acidic = universe.select_atoms(self.sel_acidic)
         basic = universe.select_atoms(self.sel_basic)
         return contacts.Contacts(
             (acidic, basic), (acidic, basic), radius=6.0, **kwargs).run(
-                start=start, stop=stop, step=step)
+                start=start, stop=stop, step=step, n_blocks=n_blocks)
 
     def test_startframe(self, universe):
         """test_startframe: TestContactAnalysis1: start frame set to 0 (resolution of
@@ -78,6 +79,11 @@ def test_startframe(self, universe):
         CA1 = self._run_Contacts(universe)
         assert len(CA1.timeseries) == universe.trajectory.n_frames
 
+    def test_uneven_blocks(self, universe):
+        """ Issue #140"""
+        CA1 = self._run_Contacts(universe, n_blocks=3)
+        assert len(CA1.timeseries) == universe.trajectory.n_frames
+
     def test_end_zero(self, universe):
         """test_end_zero: TestContactAnalysis1: stop frame 0 is not ignored"""
         CA1 = self._run_Contacts(universe, stop=0)
diff --git a/pmda/test/test_custom.py b/pmda/test/test_custom.py
index 69edd7c5..5c4373aa 100644
--- a/pmda/test/test_custom.py
+++ b/pmda/test/test_custom.py
@@ -35,6 +35,10 @@ def custom_function_tensor(ag):
     return ag.moment_of_inertia()
 
 
+def custom_function_scalar_tensor(ag):
+    return [ag.radius_of_gyration(), ag.moment_of_inertia()]
+
+
 @pytest.mark.parametrize('custom_function', [
     custom_function_vector,
     custom_function_scalar,
@@ -61,6 +65,28 @@ def test_AnalysisFromFunction(scheduler, universe, custom_function, step):
         assert_equal(results, ana.results)
 
 
+@pytest.mark.parametrize('n_blocks', [1, 3])
+def test_different_blocks(scheduler, universe, n_blocks):
+    #  This test will fail with np.hstack() #PR 88
+    ana1 = custom.AnalysisFromFunction(
+        custom_function_tensor, universe, universe.atoms).run(
+        n_blocks=n_blocks)
+    ana2 = custom.AnalysisFromFunction(
+        custom_function_tensor, universe, universe.atoms).run(
+        n_blocks=n_blocks)
+    ana3 = custom.AnalysisFromFunction(
+        custom_function_tensor, universe, universe.atoms).run(
+        n_blocks=n_blocks)
+
+    results = []
+    for ts in universe.trajectory:
+        results.append(custom_function_tensor(universe.atoms))
+    results = np.asarray(results)
+
+    for ana in (ana1, ana2, ana3):
+        assert_equal(results, ana.results)
+
+
 def custom_function_2(mobile, ref, ref2):
     return mobile.centroid() - ref.centroid() + 2 * ref2.centroid()
 
diff --git a/pmda/test/test_parallel.py b/pmda/test/test_parallel.py
index a0bed833..47b65297 100644
--- a/pmda/test/test_parallel.py
+++ b/pmda/test/test_parallel.py
@@ -56,7 +56,7 @@ def _prepare(self):
         pass
 
     def _conclude(self):
-        self.res = np.hstack(self._results)
+        self.res = np.concatenate(self._results)
 
     def _single_frame(self, ts, atomgroups):
         return ts.frame
diff --git a/pmda/test/test_rmsd.py b/pmda/test/test_rmsd.py
index 63442e05..9b8143e0 100644
--- a/pmda/test/test_rmsd.py
+++ b/pmda/test/test_rmsd.py
@@ -50,3 +50,10 @@ def test_rmsd_single_frame(self, universe, correct_values_frame_5):
         assert_almost_equal(RMSD1.rmsd, correct_values_frame_5, 4,
                             err_msg="error: rmsd profile should match " +
                             "test values")
+
+    @pytest.mark.parametrize('n_blocks', [1, 2, 3])
+    def test_rmsd_different_blocks(self, universe, n_blocks):
+        ca = universe.select_atoms("name CA")
+        universe.trajectory.rewind()
+        RMSD1 = RMSD(ca, ca).run(n_blocks=n_blocks)
+        assert_array_equal(RMSD1.rmsd.shape, (universe.trajectory.n_frames, 3))