diff --git a/releasenotes/notes/parallel-probabilistic-measurement-mps-ec00a896a0322213.yaml b/releasenotes/notes/parallel-probabilistic-measurement-mps-ec00a896a0322213.yaml
new file mode 100644
index 0000000000..1ed6a7ba01
--- /dev/null
+++ b/releasenotes/notes/parallel-probabilistic-measurement-mps-ec00a896a0322213.yaml
@@ -0,0 +1,5 @@
+---
+upgrade:
+  - |
+    Improve performance in MPS measurement method ``matrix_product_state::sample_measure_all``
+    by parallelizing the shots.
diff --git a/src/simulators/matrix_product_state/matrix_product_state.hpp b/src/simulators/matrix_product_state/matrix_product_state.hpp
index b323d95a51..6621c0371e 100644
--- a/src/simulators/matrix_product_state/matrix_product_state.hpp
+++ b/src/simulators/matrix_product_state/matrix_product_state.hpp
@@ -807,7 +807,8 @@ std::vector<reg_t> State::sample_measure_all(uint_t shots, RngEngine &rng) {
   std::vector<reg_t> all_samples;
   all_samples.resize(shots);
 
-  for (uint_t i = 0; i < shots; i++) {
+#pragma omp parallel for if (getenv("PRL_PROB_MEAS"))
+  for (int_t i = 0; i < static_cast<int_t>(shots); i++) {
     auto single_result = qreg_.sample_measure(shots, rng);
     all_samples[i] = single_result;
   }
diff --git a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp
index e190577c7c..e6cacd0239 100644
--- a/src/simulators/matrix_product_state/matrix_product_state_internal.cpp
+++ b/src/simulators/matrix_product_state/matrix_product_state_internal.cpp
@@ -1642,8 +1642,11 @@ reg_t MPS::sample_measure(uint_t shots, RngEngine &rng) const {
   reg_t current_measure(num_qubits_);
   cmatrix_t mat;
   rvector_t rnds(num_qubits_);
-  for (uint_t i = 0; i < num_qubits_; ++i) {
-    rnds[i] = rng.rand(0., 1.);
+#pragma omp critical
+  {
+    for (uint_t i = 0; i < num_qubits_; ++i) {
+      rnds[i] = rng.rand(0., 1.);
+    }
   }
   for (uint_t i = 0; i < num_qubits_; i++) {
     current_measure[i] = sample_measure_single_qubit(i, prob, rnds[i], mat);
diff --git a/test/terra/backends/aer_simulator/test_measure.py b/test/terra/backends/aer_simulator/test_measure.py
index aaf40b598e..eda06dbf7c 100644
--- a/test/terra/backends/aer_simulator/test_measure.py
+++ b/test/terra/backends/aer_simulator/test_measure.py
@@ -25,6 +25,8 @@
 from test.terra.backends.simulator_test_case import SimulatorTestCase, supported_methods
 import numpy as np
 
+import os
+
 SUPPORTED_METHODS = [
     "automatic",
     "stabilizer",
@@ -322,6 +324,22 @@ def test_mps_measure_subset_alg_qv(self):
                 result1.get_counts(circuit), result2.get_counts(circuit), delta=0.1 * shots
             )
 
+            # Test also parallel version
+            os.environ["PRL_PROB_MEAS"] = "1"
+            result2_prl = backend.run(
+                circuit, shots=shots, mps_sample_measure_algorithm="mps_probabilities"
+            ).result()
+            self.assertTrue(getattr(result2_prl, "success", "True"))
+            del os.environ["PRL_PROB_MEAS"]  # Python 3.8 in Windows
+            # os.unsetenv("PRL_PROB_MEAS")  # SInce Python 3.9
+
+            self.assertDictAlmostEqual(
+                result1.get_counts(circuit), result2_prl.get_counts(circuit), delta=0.1 * shots
+            )
+            self.assertDictAlmostEqual(
+                result2.get_counts(circuit), result2_prl.get_counts(circuit), delta=0.1 * shots
+            )
+
     def test_mps_measure_with_limited_bond_dimension(self):
         """Test MPS measure with limited bond dimension,
         where the qubits are not in sorted order