correct velocity for different groups

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,