Improve shock NEMD ensembles

doc/gpumd/input_parameters/dump_piston.rst → 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
+  dump_shock_nemd interval 1000 bin_size 20

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:

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.

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

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);
+}

src/integrate/ensemble_mirror.cuh → 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;
 };