diff --git a/doc/gpumd/input_parameters/dump_piston.rst b/doc/gpumd/input_parameters/dump_shock_nemd.rst
similarity index 64%
rename from doc/gpumd/input_parameters/dump_piston.rst
rename to doc/gpumd/input_parameters/dump_shock_nemd.rst
index 347953682..d743471a3 100644
--- a/doc/gpumd/input_parameters/dump_piston.rst
+++ b/doc/gpumd/input_parameters/dump_shock_nemd.rst
@@ -1,24 +1,23 @@
-.. _kw_dump_piston:
+.. _kw_dump_shock_nemd:
 .. index::
-   single: dump_piston (keyword in run.in)
+   single: dump_shock_nemd (keyword in run.in)
 
-:attr:`dump_piston`
-=====================
+:attr:`dump_shock_nemd`
+=======================
 
 In shock wave piston simulations, it's often crucial to compute thermo information at different regions, both before and after the shock wave passage.
 
-Piston simulations commonly involve millions of atoms. Dumping all the virial and velocity data for each atom can lead to excessively large output files, making data processing cumbersome. The `dump_piston` command addresses this by calculating spatial thermo information during the simulation.
+Piston simulations commonly involve millions of atoms. Dumping all the virial and velocity data for each atom can lead to excessively large output files, making data processing cumbersome. The `dump_shock_nemd` command addresses this by calculating spatial thermo information during the simulation.
 
-This feature calculates the spatial distribution of partical velocity (km/h), stress (GPa), temperature and density (g/cm3).
+This feature calculates the spatial distribution of partical velocity (km/h), stress (GPa), temperature and density (g/cm3) in *x*-direction.
 
 Syntax
 ------
 
 .. code::
 
-   dump_piston direction <direction> interval <time_interval> bin_size <size of each bin>
+   dump_shock_nemd interval <time_interval> bin_size <size of each bin>
 
-- The :attr:`direction` parameter specifies the shock wave direction (x, y, or z).
 - The :attr:`interval` parameter sets the output interval (number of steps).
 - The :attr:`bin_size` parameter, optional with a default value of 10, defines the thickness of each histogram bin, measured in Angstroms.
 
@@ -29,4 +28,4 @@ To output spatial thermo information every 1000 steps for a run in the x-directi
 
 .. code::
 
-  dump_piston direction x interval 1000 bin_size 20
\ No newline at end of file
+  dump_shock_nemd interval 1000 bin_size 20
\ No newline at end of file
diff --git a/doc/gpumd/input_parameters/ensemble.rst b/doc/gpumd/input_parameters/ensemble.rst
index 8378c2ac0..5559ecdb8 100644
--- a/doc/gpumd/input_parameters/ensemble.rst
+++ b/doc/gpumd/input_parameters/ensemble.rst
@@ -14,8 +14,7 @@ There are different categories of methods accessible via this keyword, which are
 * :ref:`integrators for path integral molecular dynamics simulations <kw_ensemble_pimd>`
 * :ref:`MSST integrator for simulating compressive shock wave <kw_ensemble_msst>`
 * :ref:`NPHug integrator for simulating compressive shock wave <kw_ensemble_nphug>`
-* :ref:`piston integrator for simulating compressive shock wave <kw_ensemble_piston>`
-* :ref:`momentum mirror integrator for simulating compressive shock wave <kw_ensemble_mirror>`
+* :ref:`NEMD for simulating compressive shock wave <kw_ensemble_shock_nemd>`
 
 
 .. _choice_of_parameters:
diff --git a/doc/gpumd/input_parameters/ensemble_mirror.rst b/doc/gpumd/input_parameters/ensemble_mirror.rst
deleted file mode 100644
index 584a816d8..000000000
--- a/doc/gpumd/input_parameters/ensemble_mirror.rst
+++ /dev/null
@@ -1,33 +0,0 @@
-.. _mirror:
-.. _kw_ensemble_mirror:
-.. index::
-   single: mirror (keyword in run.in)
-   single: mirror integrator
-
-:attr:`ensemble` (mirror)
-=========================
-
-This keyword is employed to configure a momentum mirror shock wave simulation, where atoms are deflected by a moving momentum mirror to generate a shock wave.
-
-The direction of the shock wave is along the x axis.
-
-We recommand to use it with the :ref:`dump_piston keyword <kw_dump_piston>`.
-
-Syntax
-------
-
-The parameters are specified as follows::
-
-    ensemble mirror vp <vp> thickness <thickness>
-
-- :attr:`<vp>`: Indicates the velocity of the moving mirror in km/s.
-- :attr:`<thickness>`: Defines the thickness of the fixed region in the other side of the cell. Its unit is Ang. This keyword is optional with the default value 20.
-
-Example
---------
-
-.. code-block:: rst
-
-    ensemble mirror vp 5 thickness 30
-
-This command makes the mirror moving at 5 km/s in x direction. The thickness of the fixed region is 30. 
\ No newline at end of file
diff --git a/doc/gpumd/input_parameters/ensemble_piston.rst b/doc/gpumd/input_parameters/ensemble_piston.rst
deleted file mode 100644
index e25c7aa4f..000000000
--- a/doc/gpumd/input_parameters/ensemble_piston.rst
+++ /dev/null
@@ -1,33 +0,0 @@
-.. _piston:
-.. _kw_ensemble_piston:
-.. index::
-   single: piston (keyword in run.in)
-   single: piston integrator
-
-:attr:`ensemble` (piston)
-=========================
-
-This keyword is employed to configure a piston shock wave simulation, where a fixed wall of atoms is displaced at a specified velocity to generate a shock wave.
-
-We recommand to use it with the :ref:`dump_piston keyword <kw_dump_piston>`.
-
-Syntax
-------
-
-The parameters are specified as follows::
-
-    ensemble piston direction <direction> vp <vp> thickness <thickness>
-
-- :attr:`<direction>`: Specifies the direction of the shock wave (`x`, `y`, or `z`).
-- :attr:`<vp>`: Indicates the velocity of the moving piston in km/s.
-- :attr:`<thickness>`: Defines the thickness of the wall in Angstroms. Note that the other side of the piston is also fixed to prevent atoms from escaping. This keyword is optional with the default value 20.
-
-Example
---------
-
-.. code-block:: rst
-
-    ensemble piston direction x vp 5 thickness 30
-
-This command makes the piston moving at 5 km/s in x direction. The thickness of the wall is 30. 
-
diff --git a/doc/gpumd/input_parameters/ensemble_shock_nemd.rst b/doc/gpumd/input_parameters/ensemble_shock_nemd.rst
new file mode 100644
index 000000000..0573d625a
--- /dev/null
+++ b/doc/gpumd/input_parameters/ensemble_shock_nemd.rst
@@ -0,0 +1,44 @@
+.. _kw_ensemble_shock_nemd:
+
+:attr:`ensemble` (NEMD shock methods)
+=====================================
+
+In a typical NEMD shock simulation, a shock wave is generated by a moving *wall*.
+The *wall* can be simulated in multiple ways:
+
+- :attr:`wall_piston`: A fixed layer of atoms.
+- :attr:`wall_mirror`: A momentum mirror that reflects atoms.
+- :attr:`wall_harmonic`: A harmonic potential that pushes atoms away. It is softer than a momentum mirror.
+
+Any of these methods can generate a shock wave. While there are other possible methods, the above methods are usually sufficient.
+
+Please note that the shock wave propagates in the *x*-direction.
+Another wall is placed on the opposite side of the cell to prevent atoms from escaping.
+It is important **NOT** to use non-periodic boundary conditions in the *x*-direction, as GPUMD currently does not handle it well. A vacuum layer is automatically added on the other side of the cell.
+
+It is recommended to use this with the :ref:`dump_shock_nemd keyword <kw_dump_shock_nemd>` to get the spatial distribution of thermo information.
+
+Syntax
+------
+
+The parameters are specified as follows:
+
+.. code-block:: rst
+
+    ensemble wall_piston vp <vp> thickness <thickness>
+
+- :attr:`<vp>`: Indicates the velocity of the moving piston in km/s.
+- :attr:`<thickness>`: Defines the thickness of the wall in Angstroms. This keyword is optional with the default value of 20.
+
+.. code-block:: rst
+
+    ensemble wall_mirror vp <vp>
+
+- :attr:`<vp>`: Indicates the velocity of the moving piston in km/s.
+
+.. code-block:: rst
+
+    ensemble wall_harmonic vp <vp> k <k>
+
+- :attr:`<vp>`: Indicates the velocity of the moving piston in km/s.
+- :attr:`<k>`: Defines the strength of the harmonic wall in eV/A^2. This keyword is optional with the default value of 10.
\ No newline at end of file
diff --git a/doc/gpumd/input_parameters/index.rst b/doc/gpumd/input_parameters/index.rst
index 05e41f39d..27b797e47 100644
--- a/doc/gpumd/input_parameters/index.rst
+++ b/doc/gpumd/input_parameters/index.rst
@@ -27,8 +27,7 @@ Below you can find a listing of keywords for the ``run.in`` input file.
    ensemble_ti_as
    ensemble_ti_rs
    ensemble_ti
-   ensemble_piston
-   ensemble_mirror
+   ensemble_shock_nemd
    ensemble_msst
    ensemble_nphug
    fix
@@ -76,4 +75,4 @@ Below you can find a listing of keywords for the ``run.in`` input file.
    dump_restart
    dump_thermo
    dump_velocity
-   dump_piston
+   dump_shock_nemd
diff --git a/src/integrate/ensemble_wall_harmonic.cu b/src/integrate/ensemble_wall_harmonic.cu
new file mode 100644
index 000000000..a5bd102a8
--- /dev/null
+++ b/src/integrate/ensemble_wall_harmonic.cu
@@ -0,0 +1,166 @@
+/*
+    Copyright 2017 Zheyong Fan, Ville Vierimaa, Mikko Ervasti, and Ari Harju
+    This file is part of GPUMD.
+    GPUMD is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    GPUMD is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+    You should have received a copy of the GNU General Public License
+    along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include "ensemble_wall_harmonic.cuh"
+
+namespace
+{
+
+static __global__ void gpu_velocity_verlet(
+  const bool is_step1,
+  const int number_of_particles,
+  const double k,
+  const double wall_pos_left,
+  const double wall_pos_right,
+  const double g_time_step,
+  const double* g_mass,
+  double* g_x,
+  double* g_y,
+  double* g_z,
+  double* g_vx,
+  double* g_vy,
+  double* g_vz,
+  double* g_fx,
+  double* g_fy,
+  double* g_fz)
+{
+  const int i = blockIdx.x * blockDim.x + threadIdx.x;
+  if (i < number_of_particles) {
+    const double time_step = g_time_step;
+    const double time_step_half = time_step * 0.5;
+    double vx = g_vx[i];
+    double vy = g_vy[i];
+    double vz = g_vz[i];
+    const double mass_inv = 1.0 / g_mass[i];
+    if (!is_step1) {
+      if (g_x[i] < wall_pos_left)
+        g_fx[i] += k * (wall_pos_left - g_x[i]);
+      if (g_x[i] > wall_pos_right)
+        g_fx[i] += k * (wall_pos_right - g_x[i]);
+    }
+    const double ax = g_fx[i] * mass_inv;
+    const double ay = g_fy[i] * mass_inv;
+    const double az = g_fz[i] * mass_inv;
+
+    vx += ax * time_step_half;
+    vy += ay * time_step_half;
+    vz += az * time_step_half;
+    g_vx[i] = vx;
+    g_vy[i] = vy;
+    g_vz[i] = vz;
+
+    if (is_step1) {
+      g_x[i] += vx * time_step;
+      g_y[i] += vy * time_step;
+      g_z[i] += vz * time_step;
+    }
+  }
+}
+} // namespace
+
+Ensemble_wall_harmonic::Ensemble_wall_harmonic(const char** params, int num_params)
+{
+  int i = 2;
+  while (i < num_params) {
+    if (strcmp(params[i], "vp") == 0) {
+      if (!is_valid_real(params[i + 1], &vp))
+        PRINT_INPUT_ERROR("Wrong inputs for vp keyword.");
+      i += 2;
+    } else if (strcmp(params[i], "k") == 0) {
+      if (!is_valid_real(params[i + 1], &k))
+        PRINT_INPUT_ERROR("Wrong inputs for k keyword.");
+      i += 2;
+    } else {
+      PRINT_INPUT_ERROR("Unknown keyword.");
+    }
+  }
+  printf("Piston velocity: %f km/s.\n", vp);
+  vp = vp / 100 * TIME_UNIT_CONVERSION;
+}
+
+void Ensemble_wall_harmonic::init()
+{
+  wall_pos_left = 0;
+  wall_pos_right = box->cpu_h[0];
+  box->cpu_h[0] += 20;
+}
+
+Ensemble_wall_harmonic::~Ensemble_wall_harmonic(void) {}
+
+void Ensemble_wall_harmonic::compute1(
+  const double time_step,
+  const std::vector<Group>& group,
+  Box& box,
+  Atom& atoms,
+  GPU_Vector<double>& thermo)
+{
+  if (*current_step == 0)
+    init();
+  find_thermo(
+    false,
+    box.get_volume(),
+    group,
+    atoms.mass,
+    atoms.potential_per_atom,
+    atoms.velocity_per_atom,
+    atoms.virial_per_atom,
+    thermo);
+  int n = atoms.number_of_atoms;
+  gpu_velocity_verlet<<<(n - 1) / 128 + 1, 128>>>(
+    true,
+    n,
+    k,
+    wall_pos_left,
+    wall_pos_right,
+    time_step,
+    atoms.mass.data(),
+    atoms.position_per_atom.data(),
+    atoms.position_per_atom.data() + n,
+    atoms.position_per_atom.data() + 2 * n,
+    atoms.velocity_per_atom.data(),
+    atoms.velocity_per_atom.data() + n,
+    atoms.velocity_per_atom.data() + 2 * n,
+    atoms.force_per_atom.data(),
+    atoms.force_per_atom.data() + n,
+    atoms.force_per_atom.data() + 2 * n);
+}
+
+void Ensemble_wall_harmonic::compute2(
+  const double time_step,
+  const std::vector<Group>& group,
+  Box& box,
+  Atom& atoms,
+  GPU_Vector<double>& thermo)
+{
+  int n = atoms.number_of_atoms;
+  wall_pos_left += time_step * vp;
+  gpu_velocity_verlet<<<(n - 1) / 128 + 1, 128>>>(
+    false,
+    n,
+    k,
+    wall_pos_left,
+    wall_pos_right,
+    time_step,
+    atoms.mass.data(),
+    atoms.position_per_atom.data(),
+    atoms.position_per_atom.data() + n,
+    atoms.position_per_atom.data() + 2 * n,
+    atoms.velocity_per_atom.data(),
+    atoms.velocity_per_atom.data() + n,
+    atoms.velocity_per_atom.data() + 2 * n,
+    atoms.force_per_atom.data(),
+    atoms.force_per_atom.data() + n,
+    atoms.force_per_atom.data() + 2 * n);
+}
diff --git a/src/integrate/ensemble_mirror.cuh b/src/integrate/ensemble_wall_harmonic.cuh
similarity index 85%
rename from src/integrate/ensemble_mirror.cuh
rename to src/integrate/ensemble_wall_harmonic.cuh
index d74e17846..57229c8f2 100644
--- a/src/integrate/ensemble_mirror.cuh
+++ b/src/integrate/ensemble_wall_harmonic.cuh
@@ -21,11 +21,11 @@
 #include "utilities/read_file.cuh"
 #include <math.h>
 
-class Ensemble_mirror : public Ensemble
+class Ensemble_wall_harmonic : public Ensemble
 {
 public:
-  Ensemble_mirror(const char** params, int num_params);
-  virtual ~Ensemble_mirror(void);
+  Ensemble_wall_harmonic(const char** params, int num_params);
+  virtual ~Ensemble_wall_harmonic(void);
 
   virtual void compute1(
     const double time_step,
@@ -44,9 +44,8 @@ public:
   void init();
 
 protected:
-  double thickness = 20;
-  double mirror_pos = 0;
+  double wall_pos_left, wall_pos_right;
   double vp;
-  GPU_Vector<bool> gpu_right_wall_list;
+  double k = 10;
   std::vector<double> thermo_cpu;
 };
\ No newline at end of file
diff --git a/src/integrate/ensemble_piston.cu b/src/integrate/ensemble_wall_mirror.cu
similarity index 56%
rename from src/integrate/ensemble_piston.cu
rename to src/integrate/ensemble_wall_mirror.cu
index 35cdf4315..c45052634 100644
--- a/src/integrate/ensemble_piston.cu
+++ b/src/integrate/ensemble_wall_mirror.cu
@@ -13,34 +13,17 @@
     along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include "ensemble_piston.cuh"
+#include "ensemble_wall_mirror.cuh"
 
 namespace
 {
-static __global__ void gpu_find_wall(
-  int number_of_atoms,
-  double left,
-  double right,
-  bool* left_wall_list,
-  bool* right_wall_list,
-  double* position)
-{
-  const int i = blockIdx.x * blockDim.x + threadIdx.x;
-  if (i < number_of_atoms) {
-    left_wall_list[i] = (position[i] < left);
-    right_wall_list[i] = (position[i] > right);
-  }
-}
-
 static __global__ void gpu_velocity_verlet(
   const bool is_step1,
   const int number_of_particles,
-  const double move_velocity_x,
-  const double move_velocity_y,
-  const double move_velocity_z,
+  const double mirror_vel_left,
+  const double mirror_pos_left,
+  const double mirror_pos_right,
   const double g_time_step,
-  bool* left_wall_list,
-  bool* right_wall_list,
   const double* g_mass,
   double* g_x,
   double* g_y,
@@ -63,57 +46,36 @@ static __global__ void gpu_velocity_verlet(
     const double ax = g_fx[i] * mass_inv;
     const double ay = g_fy[i] * mass_inv;
     const double az = g_fz[i] * mass_inv;
-    if (right_wall_list[i]) {
-      vx = 0;
-      vy = 0;
-      vz = 0;
-      g_vx[i] = 0.0;
-      g_vy[i] = 0.0;
-      g_vz[i] = 0.0;
-    } else if (left_wall_list[i]) {
-      vx = move_velocity_x;
-      vy = move_velocity_y;
-      vz = move_velocity_z;
-      g_vx[i] = 0.0;
-      g_vy[i] = 0.0;
-      g_vz[i] = 0.0;
-    } else {
-      vx += ax * time_step_half;
-      vy += ay * time_step_half;
-      vz += az * time_step_half;
-      g_vx[i] = vx;
-      g_vy[i] = vy;
-      g_vz[i] = vz;
-    }
+
+    vx += ax * time_step_half;
+    vy += ay * time_step_half;
+    vz += az * time_step_half;
+    g_vx[i] = vx;
+    g_vy[i] = vy;
+    g_vz[i] = vz;
 
     if (is_step1) {
       g_x[i] += vx * time_step;
       g_y[i] += vy * time_step;
       g_z[i] += vz * time_step;
     }
+
+    if (g_x[i] < mirror_pos_left) {
+      g_x[i] = 2 * mirror_pos_left - g_x[i];
+      g_vx[i] = 2 * mirror_vel_left - g_vx[i];
+    } else if (g_x[i] > mirror_pos_right) {
+      g_x[i] = 2 * mirror_pos_right - g_x[i];
+      g_vx[i] = -g_vx[i];
+    }
   }
 }
 } // namespace
 
-Ensemble_piston::Ensemble_piston(const char** params, int num_params)
+Ensemble_wall_mirror::Ensemble_wall_mirror(const char** params, int num_params)
 {
   int i = 2;
   while (i < num_params) {
-    if (strcmp(params[i], "direction") == 0) {
-      if (strcmp(params[i + 1], "x") == 0)
-        direction = 0;
-      else if (strcmp(params[i + 1], "y") == 0)
-        direction = 1;
-      else if (strcmp(params[i + 1], "z") == 0)
-        direction = 2;
-      else
-        PRINT_INPUT_ERROR("Direction should be x or y or z.");
-      i += 2;
-    } else if (strcmp(params[i], "thickness") == 0) {
-      if (!is_valid_real(params[i + 1], &thickness))
-        PRINT_INPUT_ERROR("Wrong inputs for thickness keyword.");
-      i += 2;
-    } else if (strcmp(params[i], "vp") == 0) {
+    if (strcmp(params[i], "vp") == 0) {
       if (!is_valid_real(params[i + 1], &vp))
         PRINT_INPUT_ERROR("Wrong inputs for vp keyword.");
       i += 2;
@@ -121,34 +83,20 @@ Ensemble_piston::Ensemble_piston(const char** params, int num_params)
       PRINT_INPUT_ERROR("Unknown keyword.");
     }
   }
-  if (direction == 0)
-    vp_x = vp / 100 * TIME_UNIT_CONVERSION;
-  else if (direction == 1)
-    vp_y = vp / 100 * TIME_UNIT_CONVERSION;
-  else if (direction == 2)
-    vp_z = vp / 100 * TIME_UNIT_CONVERSION;
-  printf("Shock wave direction: %d.\n", direction);
   printf("Piston velocity: %f km/s.\n", vp);
-  printf("The thickness of fixed wall: %f Ang.\n", thickness);
+  vp = vp / 100 * TIME_UNIT_CONVERSION;
 }
 
-void Ensemble_piston::init()
+void Ensemble_wall_mirror::init()
 {
-  int N = atom->number_of_atoms;
-  gpu_left_wall_list.resize(N, false);
-  gpu_right_wall_list.resize(N, false);
-  gpu_find_wall<<<(N - 1) / 128 + 1, 128>>>(
-    N,
-    thickness,
-    box->cpu_h[direction] - thickness,
-    gpu_left_wall_list.data(),
-    gpu_right_wall_list.data(),
-    atom->position_per_atom.data() + N * direction);
+  mirror_pos_left = 0;
+  mirror_pos_right = box->cpu_h[0];
+  box->cpu_h[0] += 20;
 }
 
-Ensemble_piston::~Ensemble_piston(void) {}
+Ensemble_wall_mirror::~Ensemble_wall_mirror(void) {}
 
-void Ensemble_piston::compute1(
+void Ensemble_wall_mirror::compute1(
   const double time_step,
   const std::vector<Group>& group,
   Box& box,
@@ -170,12 +118,10 @@ void Ensemble_piston::compute1(
   gpu_velocity_verlet<<<(n - 1) / 128 + 1, 128>>>(
     true,
     n,
-    vp_x,
-    vp_y,
-    vp_z,
+    vp,
+    mirror_pos_left,
+    mirror_pos_right,
     time_step,
-    gpu_left_wall_list.data(),
-    gpu_right_wall_list.data(),
     atoms.mass.data(),
     atoms.position_per_atom.data(),
     atoms.position_per_atom.data() + n,
@@ -188,7 +134,7 @@ void Ensemble_piston::compute1(
     atoms.force_per_atom.data() + 2 * n);
 }
 
-void Ensemble_piston::compute2(
+void Ensemble_wall_mirror::compute2(
   const double time_step,
   const std::vector<Group>& group,
   Box& box,
@@ -196,15 +142,15 @@ void Ensemble_piston::compute2(
   GPU_Vector<double>& thermo)
 {
   int n = atoms.number_of_atoms;
+
+  mirror_pos_left += time_step * vp;
   gpu_velocity_verlet<<<(n - 1) / 128 + 1, 128>>>(
     false,
     n,
-    vp_x,
-    vp_y,
-    vp_z,
+    vp,
+    mirror_pos_left,
+    mirror_pos_right,
     time_step,
-    gpu_left_wall_list.data(),
-    gpu_right_wall_list.data(),
     atoms.mass.data(),
     atoms.position_per_atom.data(),
     atoms.position_per_atom.data() + n,
@@ -215,4 +161,4 @@ void Ensemble_piston::compute2(
     atoms.force_per_atom.data(),
     atoms.force_per_atom.data() + n,
     atoms.force_per_atom.data() + 2 * n);
-}
\ No newline at end of file
+}
diff --git a/src/integrate/ensemble_wall_mirror.cuh b/src/integrate/ensemble_wall_mirror.cuh
new file mode 100644
index 000000000..6fa3fcdd1
--- /dev/null
+++ b/src/integrate/ensemble_wall_mirror.cuh
@@ -0,0 +1,50 @@
+/*
+    Copyright 2017 Zheyong Fan, Ville Vierimaa, Mikko Ervasti, and Ari Harju
+    This file is part of GPUMD.
+    GPUMD is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    GPUMD is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+    You should have received a copy of the GNU General Public License
+    along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#pragma once
+#include "ensemble.cuh"
+#include "model/box.cuh"
+#include "utilities/common.cuh"
+#include "utilities/error.cuh"
+#include "utilities/read_file.cuh"
+#include <math.h>
+
+class Ensemble_wall_mirror : public Ensemble
+{
+public:
+  Ensemble_wall_mirror(const char** params, int num_params);
+  virtual ~Ensemble_wall_mirror(void);
+
+  virtual void compute1(
+    const double time_step,
+    const std::vector<Group>& group,
+    Box& box,
+    Atom& atoms,
+    GPU_Vector<double>& thermo);
+
+  virtual void compute2(
+    const double time_step,
+    const std::vector<Group>& group,
+    Box& box,
+    Atom& atoms,
+    GPU_Vector<double>& thermo);
+
+  void init();
+
+protected:
+  double mirror_pos_left, mirror_pos_right;
+  double vp;
+  std::vector<double> thermo_cpu;
+};
\ No newline at end of file
diff --git a/src/integrate/ensemble_mirror.cu b/src/integrate/ensemble_wall_piston.cu
similarity index 76%
rename from src/integrate/ensemble_mirror.cu
rename to src/integrate/ensemble_wall_piston.cu
index 509fd4048..4def28f6f 100644
--- a/src/integrate/ensemble_mirror.cu
+++ b/src/integrate/ensemble_wall_piston.cu
@@ -13,25 +13,32 @@
     along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include "ensemble_mirror.cuh"
+#include "ensemble_wall_piston.cuh"
 
 namespace
 {
-static __global__ void
-gpu_find_wall(int number_of_atoms, double right, bool* right_wall_list, double* position)
+static __global__ void gpu_find_wall(
+  int number_of_atoms,
+  double wall_pos_left,
+  double wall_pos_right,
+  bool* wall_list_left,
+  bool* wall_list_right,
+  double* position)
 {
   const int i = blockIdx.x * blockDim.x + threadIdx.x;
-  if (i < number_of_atoms)
-    right_wall_list[i] = (position[i] > right);
+  if (i < number_of_atoms) {
+    wall_list_left[i] = (position[i] < wall_pos_left);
+    wall_list_right[i] = (position[i] > wall_pos_right);
+  }
 }
 
 static __global__ void gpu_velocity_verlet(
   const bool is_step1,
   const int number_of_particles,
-  const double mirror_vel,
-  const double mirror_pos,
-  const double g_time_step,
+  const double vp,
+  bool* left_wall_list,
   bool* right_wall_list,
+  const double g_time_step,
   const double* g_mass,
   double* g_x,
   double* g_y,
@@ -54,66 +61,72 @@ static __global__ void gpu_velocity_verlet(
     const double ax = g_fx[i] * mass_inv;
     const double ay = g_fy[i] * mass_inv;
     const double az = g_fz[i] * mass_inv;
+
     if (right_wall_list[i]) {
       vx = 0;
       vy = 0;
       vz = 0;
-      g_vx[i] = 0.0;
-      g_vy[i] = 0.0;
-      g_vz[i] = 0.0;
+    } else if (left_wall_list[i]) {
+      vx = vp;
+      vy = 0;
+      vz = 0;
     } else {
       vx += ax * time_step_half;
       vy += ay * time_step_half;
       vz += az * time_step_half;
-      g_vx[i] = vx;
-      g_vy[i] = vy;
-      g_vz[i] = vz;
     }
+    g_vx[i] = vx;
+    g_vy[i] = vy;
+    g_vz[i] = vz;
 
     if (is_step1) {
       g_x[i] += vx * time_step;
       g_y[i] += vy * time_step;
       g_z[i] += vz * time_step;
-    } else if (g_x[i] < mirror_pos) {
-      g_x[i] = 2 * mirror_pos - g_x[i];
-      g_vx[i] = 2 * mirror_vel - g_vx[i];
     }
   }
 }
 } // namespace
 
-Ensemble_mirror::Ensemble_mirror(const char** params, int num_params)
+Ensemble_wall_piston::Ensemble_wall_piston(const char** params, int num_params)
 {
   int i = 2;
   while (i < num_params) {
-    if (strcmp(params[i], "thickness") == 0) {
-      if (!is_valid_real(params[i + 1], &thickness))
-        PRINT_INPUT_ERROR("Wrong inputs for thickness keyword.");
-      i += 2;
-    } else if (strcmp(params[i], "vp") == 0) {
+    if (strcmp(params[i], "vp") == 0) {
       if (!is_valid_real(params[i + 1], &vp))
         PRINT_INPUT_ERROR("Wrong inputs for vp keyword.");
       i += 2;
+    } else if (strcmp(params[i], "thickness") == 0) {
+      if (!is_valid_real(params[i + 1], &thickness))
+        PRINT_INPUT_ERROR("Wrong inputs for thickness keyword.");
+      i += 2;
     } else {
       PRINT_INPUT_ERROR("Unknown keyword.");
     }
   }
-  printf("Piston velocity: %f km/s.\n", vp);
-  printf("The thickness of fixed wall on the right side: %f Ang.\n", thickness);
   vp = vp / 100 * TIME_UNIT_CONVERSION;
+  printf("Piston velocity: %f km/s.\n", vp);
+  printf("The thickness of fixed wall: %f Ang.\n", thickness);
 }
 
-void Ensemble_mirror::init()
+void Ensemble_wall_piston::init()
 {
   int N = atom->number_of_atoms;
+  gpu_left_wall_list.resize(N, false);
   gpu_right_wall_list.resize(N, false);
   gpu_find_wall<<<(N - 1) / 128 + 1, 128>>>(
-    N, box->cpu_h[0] - thickness, gpu_right_wall_list.data(), atom->position_per_atom.data());
+    N,
+    thickness,
+    box->cpu_h[0] - thickness,
+    gpu_left_wall_list.data(),
+    gpu_right_wall_list.data(),
+    atom->position_per_atom.data());
+  box->cpu_h[0] += 20;
 }
 
-Ensemble_mirror::~Ensemble_mirror(void) {}
+Ensemble_wall_piston::~Ensemble_wall_piston(void) {}
 
-void Ensemble_mirror::compute1(
+void Ensemble_wall_piston::compute1(
   const double time_step,
   const std::vector<Group>& group,
   Box& box,
@@ -136,9 +149,9 @@ void Ensemble_mirror::compute1(
     true,
     n,
     vp,
-    mirror_pos,
-    time_step,
+    gpu_left_wall_list.data(),
     gpu_right_wall_list.data(),
+    time_step,
     atoms.mass.data(),
     atoms.position_per_atom.data(),
     atoms.position_per_atom.data() + n,
@@ -151,7 +164,7 @@ void Ensemble_mirror::compute1(
     atoms.force_per_atom.data() + 2 * n);
 }
 
-void Ensemble_mirror::compute2(
+void Ensemble_wall_piston::compute2(
   const double time_step,
   const std::vector<Group>& group,
   Box& box,
@@ -159,14 +172,13 @@ void Ensemble_mirror::compute2(
   GPU_Vector<double>& thermo)
 {
   int n = atoms.number_of_atoms;
-  mirror_pos += time_step * vp;
   gpu_velocity_verlet<<<(n - 1) / 128 + 1, 128>>>(
     false,
     n,
     vp,
-    mirror_pos,
-    time_step,
+    gpu_left_wall_list.data(),
     gpu_right_wall_list.data(),
+    time_step,
     atoms.mass.data(),
     atoms.position_per_atom.data(),
     atoms.position_per_atom.data() + n,
@@ -177,4 +189,4 @@ void Ensemble_mirror::compute2(
     atoms.force_per_atom.data(),
     atoms.force_per_atom.data() + n,
     atoms.force_per_atom.data() + 2 * n);
-}
+}
\ No newline at end of file
diff --git a/src/integrate/ensemble_piston.cuh b/src/integrate/ensemble_wall_piston.cuh
similarity index 87%
rename from src/integrate/ensemble_piston.cuh
rename to src/integrate/ensemble_wall_piston.cuh
index b0d0f37a8..51d0c8f8d 100644
--- a/src/integrate/ensemble_piston.cuh
+++ b/src/integrate/ensemble_wall_piston.cuh
@@ -21,11 +21,11 @@
 #include "utilities/read_file.cuh"
 #include <math.h>
 
-class Ensemble_piston : public Ensemble
+class Ensemble_wall_piston : public Ensemble
 {
 public:
-  Ensemble_piston(const char** params, int num_params);
-  virtual ~Ensemble_piston(void);
+  Ensemble_wall_piston(const char** params, int num_params);
+  virtual ~Ensemble_wall_piston(void);
 
   virtual void compute1(
     const double time_step,
@@ -44,9 +44,8 @@ public:
   void init();
 
 protected:
-  int direction;
   double thickness = 20;
-  double vp, vp_x = 0, vp_y = 0, vp_z = 0;
+  double vp;
   GPU_Vector<bool> gpu_left_wall_list, gpu_right_wall_list;
   std::vector<double> thermo_cpu;
 };
\ No newline at end of file
diff --git a/src/integrate/integrate.cu b/src/integrate/integrate.cu
index 4236588e5..59a8b4a0d 100644
--- a/src/integrate/integrate.cu
+++ b/src/integrate/integrate.cu
@@ -21,7 +21,6 @@ The driver class for the various integrators.
 #include "ensemble_bdp.cuh"
 #include "ensemble_ber.cuh"
 #include "ensemble_lan.cuh"
-#include "ensemble_mirror.cuh"
 #include "ensemble_msst.cuh"
 #include "ensemble_mttk.cuh"
 #include "ensemble_nhc.cuh"
@@ -29,11 +28,13 @@ The driver class for the various integrators.
 #include "ensemble_npt_scr.cuh"
 #include "ensemble_nve.cuh"
 #include "ensemble_pimd.cuh"
-#include "ensemble_piston.cuh"
 #include "ensemble_ti.cuh"
 #include "ensemble_ti_as.cuh"
 #include "ensemble_ti_rs.cuh"
 #include "ensemble_ti_spring.cuh"
+#include "ensemble_wall_harmonic.cuh"
+#include "ensemble_wall_mirror.cuh"
+#include "ensemble_wall_piston.cuh"
 #include "integrate.cuh"
 #include "model/atom.cuh"
 #include "utilities/common.cuh"
@@ -135,6 +136,8 @@ void Integrate::initialize(
       break;
     case -9: // ti_as
       break;
+    case -10:
+      break;
     case 21: // heat-NHC
       ensemble.reset(new Ensemble_NHC(
         type,
@@ -406,24 +409,27 @@ void Integrate::parse_ensemble(
   } else if (strcmp(param[1], "ti_spring") == 0) {
     type = -2;
     ensemble.reset(new Ensemble_TI_Spring(param, num_param));
-  } else if (strcmp(param[1], "piston") == 0) {
+  } else if (strcmp(param[1], "wall_piston") == 0) {
     type = -4;
-    ensemble.reset(new Ensemble_piston(param, num_param));
+    ensemble.reset(new Ensemble_wall_piston(param, num_param));
   } else if (strcmp(param[1], "nphug") == 0) {
     type = -5;
     ensemble.reset(new Ensemble_NPHug(param, num_param));
   } else if (strcmp(param[1], "ti") == 0) {
     type = -6;
     ensemble.reset(new Ensemble_TI(param, num_param));
-  } else if (strcmp(param[1], "mirror") == 0) {
+  } else if (strcmp(param[1], "wall_mirror") == 0) {
     type = -7;
-    ensemble.reset(new Ensemble_mirror(param, num_param));
+    ensemble.reset(new Ensemble_wall_mirror(param, num_param));
   } else if (strcmp(param[1], "ti_rs") == 0) {
     type = -8;
     ensemble.reset(new Ensemble_TI_RS(param, num_param));
   } else if (strcmp(param[1], "ti_as") == 0) {
     type = -9;
     ensemble.reset(new Ensemble_TI_AS(param, num_param));
+  } else if (strcmp(param[1], "wall_harmonic") == 0) {
+    type = -10;
+    ensemble.reset(new Ensemble_wall_harmonic(param, num_param));
   } else {
     PRINT_INPUT_ERROR("Invalid ensemble type.");
   }
@@ -862,6 +868,8 @@ void Integrate::parse_ensemble(
       break;
     case -9:
       break;
+    case -10:
+      break;
     case 21:
       printf("Integrate with heating and cooling for this run.\n");
       printf("    choose the Nose-Hoover chain method.\n");
diff --git a/src/main_gpumd/run.cu b/src/main_gpumd/run.cu
index c406d2d17..ee266fe3e 100644
--- a/src/main_gpumd/run.cu
+++ b/src/main_gpumd/run.cu
@@ -421,8 +421,8 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
     measure.dump_beads.parse(param, num_param);
   } else if (strcmp(param[0], "dump_observer") == 0) {
     measure.dump_observer.parse(param, num_param);
-  } else if (strcmp(param[0], "dump_piston") == 0) {
-    measure.dump_piston.parse(param, num_param);
+  } else if (strcmp(param[0], "dump_shock_nemd") == 0) {
+    measure.dump_shock_nemd.parse(param, num_param);
   } else if (strcmp(param[0], "dump_dipole") == 0) {
     measure.dump_dipole.parse(param, num_param);
   } else if (strcmp(param[0], "dump_polarizability") == 0) {
diff --git a/src/measure/dump_piston.cu b/src/measure/dump_shock_nemd.cu
similarity index 87%
rename from src/measure/dump_piston.cu
rename to src/measure/dump_shock_nemd.cu
index 058419029..f42c51bb0 100644
--- a/src/measure/dump_piston.cu
+++ b/src/measure/dump_shock_nemd.cu
@@ -13,7 +13,7 @@
     along with GPUMD.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include "dump_piston.cuh"
+#include "dump_shock_nemd.cuh"
 #include <cstring>
 
 namespace
@@ -139,24 +139,14 @@ void write_to_file(FILE* file, double* array, int n)
 
 } // namespace
 
-void Dump_Piston::parse(const char** param, int num_param)
+void Dump_Shock_NEMD::parse(const char** param, int num_param)
 {
   dump_ = true;
 
   printf("Dump spatial histogram thermo information for piston shock wave simulation, ");
   int i = 1;
   while (i < num_param) {
-    if (strcmp(param[i], "direction") == 0) {
-      if (strcmp(param[i + 1], "x") == 0)
-        direction = 0;
-      else if (strcmp(param[i + 1], "y") == 0)
-        direction = 1;
-      else if (strcmp(param[i + 1], "z") == 0)
-        direction = 2;
-      else
-        PRINT_INPUT_ERROR("Direction should be x or y or z.");
-      i += 2;
-    } else if (strcmp(param[i], "interval") == 0) {
+    if (strcmp(param[i], "interval") == 0) {
       if (!is_valid_int(param[i + 1], &dump_interval_))
         PRINT_INPUT_ERROR("Dump interval should be an integer.");
       i += 2;
@@ -168,19 +158,9 @@ void Dump_Piston::parse(const char** param, int num_param)
       PRINT_INPUT_ERROR("Unknown keyword.");
     }
   }
-
-  if (strcmp(param[2], "x") == 0) {
-    direction = 0;
-  } else if (strcmp(param[2], "y") == 0) {
-    direction = 1;
-  } else if (strcmp(param[2], "z") == 0) {
-    direction = 2;
-  } else
-    PRINT_INPUT_ERROR("Direction should be x or y or z.");
-  printf("in %d direction.\n", direction);
 }
 
-void Dump_Piston::preprocess(Atom& atom, Box& box)
+void Dump_Shock_NEMD::preprocess(Atom& atom, Box& box)
 {
   if (!dump_)
     return;
@@ -199,12 +179,10 @@ void Dump_Piston::preprocess(Atom& atom, Box& box)
   pyy_file = my_fopen("pyy_hist.txt", "w");
   pzz_file = my_fopen("pzz_hist.txt", "w");
   density_file = my_fopen("density_hist.txt", "w");
-  com_vx_file = my_fopen("com_vx_hist.txt", "w");
-  com_vy_file = my_fopen("com_vy_hist.txt", "w");
-  com_vz_file = my_fopen("com_vz_hist.txt", "w");
+  com_vx_file = my_fopen("vp_hist.txt", "w");
 }
 
-void Dump_Piston::process(Atom& atom, Box& box, const int step)
+void Dump_Shock_NEMD::process(Atom& atom, Box& box, const int step)
 {
   if (!dump_ || step % dump_interval_ != 0)
     return;
@@ -295,11 +273,9 @@ void Dump_Piston::process(Atom& atom, Box& box, const int step)
   write_to_file(pzz_file, cpu_pzz.data(), bins);
   write_to_file(density_file, cpu_density.data(), bins);
   write_to_file(com_vx_file, cpu_com_vx.data(), bins);
-  write_to_file(com_vy_file, cpu_com_vy.data(), bins);
-  write_to_file(com_vz_file, cpu_com_vz.data(), bins);
 }
 
-void Dump_Piston::postprocess()
+void Dump_Shock_NEMD::postprocess()
 {
   if (!dump_)
     return;
@@ -310,7 +286,5 @@ void Dump_Piston::postprocess()
   fclose(pzz_file);
   fclose(density_file);
   fclose(com_vx_file);
-  fclose(com_vy_file);
-  fclose(com_vz_file);
   dump_ = false;
 }
\ No newline at end of file
diff --git a/src/measure/dump_piston.cuh b/src/measure/dump_shock_nemd.cuh
similarity index 94%
rename from src/measure/dump_piston.cuh
rename to src/measure/dump_shock_nemd.cuh
index 91271b6c8..bd077a589 100644
--- a/src/measure/dump_piston.cuh
+++ b/src/measure/dump_shock_nemd.cuh
@@ -24,7 +24,7 @@
 #include <string>
 #include <vector>
 
-class Dump_Piston
+class Dump_Shock_NEMD
 {
 public:
   void parse(const char** param, int num_param);
@@ -36,12 +36,11 @@ private:
   bool dump_ = false;
   int n;
   int dump_interval_ = -1;
-  int direction;
+  int direction = 0;
   int bins;
   double slice_vol = 1;
   double avg_window = 10;
-  FILE *temp_file, *pxx_file, *pyy_file, *pzz_file, *density_file, *com_vx_file, *com_vy_file,
-    *com_vz_file;
+  FILE *temp_file, *pxx_file, *pyy_file, *pzz_file, *density_file, *com_vx_file;
   GPU_Vector<double> gpu_temp, gpu_pxx, gpu_pyy, gpu_pzz, gpu_density, gpu_com_vx, gpu_com_vy,
     gpu_com_vz, gpu_number;
   std::vector<double> cpu_temp, cpu_pxx, cpu_pyy, cpu_pzz, cpu_density, cpu_com_vx, cpu_com_vy,
diff --git a/src/measure/measure.cu b/src/measure/measure.cu
index 6920ea3fb..a5d195950 100644
--- a/src/measure/measure.cu
+++ b/src/measure/measure.cu
@@ -55,7 +55,7 @@ void Measure::initialize(
   dump_exyz.preprocess(number_of_atoms);
   dump_beads.preprocess(number_of_atoms, atom.number_of_beads);
   dump_observer.preprocess(number_of_atoms, number_of_potentials, force);
-  dump_piston.preprocess(atom, box);
+  dump_shock_nemd.preprocess(atom, box);
   dump_dipole.preprocess(number_of_atoms, number_of_potentials, force);
   dump_polarizability.preprocess(number_of_atoms, number_of_potentials, force);
   active.preprocess(number_of_atoms, number_of_potentials, force);
@@ -85,7 +85,7 @@ void Measure::finalize(
   dump_exyz.postprocess();
   dump_beads.postprocess();
   dump_observer.postprocess();
-  dump_piston.postprocess();
+  dump_shock_nemd.postprocess();
   dump_dipole.postprocess();
   dump_polarizability.postprocess();
   active.postprocess();
@@ -204,7 +204,7 @@ void Measure::process(
     step, temperature, box.get_volume(), hnemd.fe, atom.velocity_per_atom, atom.virial_per_atom);
 
   lsqt.process(atom, box, step);
-  dump_piston.process(atom, box, step);
+  dump_shock_nemd.process(atom, box, step);
 
 #ifdef USE_NETCDF
   dump_netcdf.process(
diff --git a/src/measure/measure.cuh b/src/measure/measure.cuh
index c3ee2ce51..bf62380da 100644
--- a/src/measure/measure.cuh
+++ b/src/measure/measure.cuh
@@ -22,10 +22,10 @@
 #include "dump_exyz.cuh"
 #include "dump_force.cuh"
 #include "dump_observer.cuh"
-#include "dump_piston.cuh"
 #include "dump_polarizability.cuh"
 #include "dump_position.cuh"
 #include "dump_restart.cuh"
+#include "dump_shock_nemd.cuh"
 #include "dump_thermo.cuh"
 #include "dump_velocity.cuh"
 #include "force/force.cuh"
@@ -109,7 +109,7 @@ public:
   Dump_EXYZ dump_exyz;
   Dump_Beads dump_beads;
   Dump_Observer dump_observer;
-  Dump_Piston dump_piston;
+  Dump_Shock_NEMD dump_shock_nemd;
   Dump_Dipole dump_dipole;
   Dump_Polarizability dump_polarizability;
   Active active;