Merge pull request #668 from brucefan1983/correct_velocity

correct velocity for different groups

doc/gpumd/input_parameters/correct_velocity.rst

@@ -4,21 +4,32 @@
 :attr:`correct_velocity`
 ========================
 
-This keyword allows one to enforce zero angular momentum.
-This can be useful as the stochastic nature of most thermostats and barostats available in :program:`gpumd` can cause the conservation of angular momentum to be violated.
+This keyword allows one to enforce zero linear and angular momenta at regular intervals.
 
 Syntax
 ------
-* This keyword only has one parameter, which is the interval between which the angular momentum is set to zero::
+* This keyword can have one or two parameters, which are the interval between which the linear and angular momenta are set to zero and an optional group method::
   
-    correct_velocity <interval>
+    correct_velocity <interval> [<group_method>]
 
+* :attr:`interval` must be larger than or equal to 10. A value between 10 and 100 should be good.
 
-Example
--------
+* The :attr:`group_method` must be defined in :attr:`model.xyz`.
+
+Example 1
+---------
 
 The command::
 
     correct_velocity 10
 
-implies that the angular momentum is strictly set to zero every 10 steps.
+implies that the linear and angular momenta are set to zero every 10 steps.
+
+Example 2
+---------
+
+The command::
+
+    correct_velocity 50 3
+
+implies that the linear and angular momenta for the individual groups in group method 3 are set to zero every 50 steps.

src/force/nep3.cu

@@ -194,8 +194,8 @@ NEP3::NEP3(const char* file_potential, const int num_atoms)
     paramb.use_typewise_cutoff_zbl = get_int_from_token(tokens[6], __FILE__, __LINE__);
   }
 #ifdef USE_TABLE
-  if (paramd.use_typewise_cutoff) {
-    PRINT_ERROR("Cannot use tabulated radial functions with typewise cutoff.")
+  if (paramb.use_typewise_cutoff) {
+    PRINT_INPUT_ERROR("Cannot use tabulated radial functions with typewise cutoff.");
   }
 #endif
 

src/force/nep3_multigpu.cu

@@ -196,8 +196,8 @@ NEP3_MULTIGPU::NEP3_MULTIGPU(
     paramb.use_typewise_cutoff_zbl = get_int_from_token(tokens[6], __FILE__, __LINE__);
   }
 #ifdef USE_TABLE
-  if (paramd.use_typewise_cutoff) {
-    PRINT_ERROR("Cannot use tabulated radial functions with typewise cutoff.")
+  if (paramb.use_typewise_cutoff) {
+    PRINT_INPUT_ERROR("Cannot use tabulated radial functions with typewise cutoff.");
   }
 #endif
 

src/main_gpumd/run.cu

@@ -282,6 +282,7 @@ void Run::perform_a_run()
 
     velocity.correct_velocity(
       step,
+      group,
       atom.cpu_mass,
       atom.position_per_atom,
       atom.cpu_position_per_atom,
@@ -398,7 +399,7 @@ void Run::parse_one_keyword(std::vector<std::string>& tokens)
   } else if (strcmp(param[0], "time_step") == 0) {
     parse_time_step(param, num_param);
   } else if (strcmp(param[0], "correct_velocity") == 0) {
-    parse_correct_velocity(param, num_param);
+    parse_correct_velocity(param, num_param, group);
   } else if (strcmp(param[0], "dump_thermo") == 0) {
     measure.dump_thermo.parse(param, num_param);
   } else if (strcmp(param[0], "dump_position") == 0) {
@@ -517,17 +518,40 @@ void Run::parse_velocity(const char** param, int num_param)
   }
 }
 
-void Run::parse_correct_velocity(const char** param, int num_param)
+void Run::parse_correct_velocity(const char** param, int num_param, const std::vector<Group>& group)
 {
-  if (num_param != 2) {
-    PRINT_INPUT_ERROR("correct_velocity should have 1 parameter.\n");
+  printf("Correct linear and angular momenta.\n");
+
+  if (num_param != 2 && num_param != 3) {
+    PRINT_INPUT_ERROR("correct_velocity should have 1 or 2 parameters.\n");
   }
   if (!is_valid_int(param[1], &velocity.velocity_correction_interval)) {
     PRINT_INPUT_ERROR("velocity correction interval should be an integer.\n");
   }
-  if (velocity.velocity_correction_interval <= 0) {
-    PRINT_INPUT_ERROR("velocity correction interval should be positive.\n");
+  if (velocity.velocity_correction_interval < 10) {
+    PRINT_INPUT_ERROR("velocity correction interval should >= 10.\n");
+  }
+
+  printf("    every %d steps.\n", velocity.velocity_correction_interval);
+
+  if (num_param == 3) {
+    if (!is_valid_int(param[2], &velocity.velocity_correction_group_method)) {
+      PRINT_INPUT_ERROR("velocity correction group method should be an integer.\n");
+    }
+    if (velocity.velocity_correction_group_method < 0) {
+      PRINT_INPUT_ERROR("grouping method should >= 0.\n");
+    }
+    if (velocity.velocity_correction_group_method >= group.size()) {
+      PRINT_INPUT_ERROR("grouping method should < maximum number of grouping methods.\n");
+    }
+  }
+
+  if (velocity.velocity_correction_group_method < 0) {
+    printf("    for the whole system.\n");
+  } else {
+    printf("    for individual groups in group method %d.\n", velocity.velocity_correction_group_method);
   }
+
   velocity.do_velocity_correction = true;
 }
 

src/main_gpumd/run.cuh

@@ -48,7 +48,7 @@ private:
   void parse_neighbor(const char** param, int num_param);
   void parse_velocity(const char** param, int num_param);
   void parse_change_box(const char** param, int num_param);
-  void parse_correct_velocity(const char** param, int num_param);
+  void parse_correct_velocity(const char** param, int num_param, const std::vector<Group>& group);
   void parse_time_step(const char** param, int num_param);
   void parse_run(const char** param, int num_param);
 

src/main_gpumd/velocity.cu

@@ -21,6 +21,7 @@ If DEBUG is on in the makefile, the velocities are the same from run to run.
 If DEBUG is off, the velocities are different in different runs.
 ------------------------------------------------------------------------------*/
 
+#include "model/group.cuh"
 #include "utilities/common.cuh"
 #include "utilities/gpu_vector.cuh"
 #include "velocity.cuh"
@@ -269,6 +270,7 @@ void Velocity::correct_velocity(
 
 void Velocity::correct_velocity(
   const int step,
+  const std::vector<Group>& group,
   const std::vector<double>& cpu_mass,
   GPU_Vector<double>& position_per_atom,
   std::vector<double>& cpu_position_per_atom,
@@ -279,7 +281,32 @@ void Velocity::correct_velocity(
     if ((step + 1) % velocity_correction_interval == 0) {
       position_per_atom.copy_to_host(cpu_position_per_atom.data());
       velocity_per_atom.copy_to_host(cpu_velocity_per_atom.data());
-      correct_velocity(cpu_mass, cpu_position_per_atom, cpu_velocity_per_atom);
+      if (velocity_correction_group_method < 0) {
+        correct_velocity(cpu_mass, cpu_position_per_atom, cpu_velocity_per_atom);
+      } else {
+        for (int g = 0; g < group[velocity_correction_group_method].number; ++g) {
+          int cpu_size = group[velocity_correction_group_method].cpu_size[g];
+          int cpu_size_sum = group[velocity_correction_group_method].cpu_size_sum[g];
+          std::vector<double> mass(cpu_size);
+          std::vector<double> position(cpu_size * 3);
+          std::vector<double> velocity(cpu_size * 3);
+          for (int m = 0; m < cpu_size; ++m) {
+            int n = group[velocity_correction_group_method].cpu_contents[cpu_size_sum + m];
+            mass[m] = cpu_mass[n];
+            for (int d = 0; d < 3; ++d) {
+              position[m + d * cpu_size] = cpu_position_per_atom[n + d * cpu_mass.size()];
+              velocity[m + d * cpu_size] = cpu_velocity_per_atom[n + d * cpu_mass.size()];
+            }
+          }
+          correct_velocity(mass, position, velocity);
+          for (int m = 0; m < cpu_size; ++m) {
+            int n = group[velocity_correction_group_method].cpu_contents[cpu_size_sum + m];
+            for (int d = 0; d < 3; ++d) {
+              cpu_velocity_per_atom[n + d * cpu_mass.size()] = velocity[m + d * cpu_size];
+            }
+          }
+        }
+      }
       velocity_per_atom.copy_from_host(cpu_velocity_per_atom.data());
     }
   }
@@ -324,4 +351,8 @@ void Velocity::initialize(
   velocity_per_atom.copy_from_host(cpu_velocity_per_atom.data());
 }
 
-void Velocity::finalize() { do_velocity_correction = false; }
+void Velocity::finalize()
+{
+  do_velocity_correction = false; 
+  velocity_correction_group_method = -1;
+}

src/main_gpumd/velocity.cuh

@@ -18,11 +18,14 @@
 #include "utilities/gpu_vector.cuh"
 #include <vector>
 
+class Group;
+
 class Velocity
 {
 public:
   bool do_velocity_correction = false;
   int velocity_correction_interval = 1000;
+  int velocity_correction_group_method = -1;
 
   void initialize(
     const bool has_velocity_in_xyz,
@@ -36,6 +39,7 @@ public:
 
   void correct_velocity(
     const int step,
+    const std::vector<Group>& group,
     const std::vector<double>& cpu_mass,
     GPU_Vector<double>& position_per_atom,
     std::vector<double>& cpu_position_per_atom,