Merge pull request brucefan1983#1409 from luowh35/expand_group_method

Expand group method

Author: brucefan1983

doc/gpumd/input_parameters/fix.rst

@@ -9,13 +9,28 @@ This keyword can be used to fix (freeze) a group of atoms
 
 Syntax
 ------
-This keyword requires a single parameter which is the label of the group in which the atoms are to be fixed (velocities and forces are always set to zero such that the atoms in the group do not move).
+This keyword accepts 1 or 2 parameters.
+The atoms in the specified group will be fixed (velocities and forces are always set to zero such that the atoms do not move).
 The full command reads::
 
   fix <group_label>
+  fix <grouping_method> <group_label>
 
-Here, the :attr:`group_label` refers to the grouping method 0 defined in the :ref:`simulation model file <model_xyz>`.
+- If only :attr:`group_label` is given, grouping method 0 is used by default.
+- If :attr:`grouping_method` is also specified, the given grouping method defined in the :ref:`simulation model file <model_xyz>` will be used.
+
+Example
+-------
+
+Fix group 0 using the default grouping method 0::
+
+  fix 0
+
+Fix group 2 using grouping method 1::
+
+  fix 1 2
 
 Caveats
 -------
-This keyword is not propagating, which means that it only affects the simulation within the run it belongs to.
+* This keyword is not propagating, which means that it only affects the simulation within the run it belongs to.
+* When both :attr:`fix` and :ref:`move <kw_move>` are used, they must use the same grouping method.

doc/gpumd/input_parameters/move.rst

@@ -10,12 +10,15 @@ This keyword is used to move part of the system with a constant velocity.
 Syntax
 ------
 
-The :attr:`move` keyword requires 4 parameters::
+The :attr:`move` keyword accepts 4 or 5 parameters::
 
   move <moving_group_id> <velocity_x> <velocity_y> <velocity_z>
+  move <grouping_method> <moving_group_id> <velocity_x> <velocity_y> <velocity_z>
 
-Here, :attr:`moving_group_id` specifies the group id for the moving part, which should be defined in the grouping method 0.
-The next three parameters specify the moving velocity vector, in units of Ångstrom/fs.
+- If only :attr:`moving_group_id` and velocities are given, grouping method 0 is used by default.
+- If :attr:`grouping_method` is also specified, the given grouping method defined in the :ref:`simulation model file <model_xyz>` will be used.
+
+The last three parameters specify the moving velocity vector, in units of Ångstrom/fs.
 
 
 Example
@@ -28,12 +31,19 @@ One can first equilibrate the system and then move one group of atoms and at the
   run 1000000
 
   # production stage
-  ensemble nvt_scr 300 300 100
-  fix  0            # fix atoms in group 0
-  move 1 0.001 0 0  # move atoms in group 1, with a speed of 0.001 Ångstrom/fs in the x direction
+  ensemble nvt_ber 300 300 100
+  fix  0            # fix atoms in group 0 (default grouping method 0)
+  move 1 0.001 0 0  # move atoms in group 1, with a speed of 0.001 Ångstrom/fs in the x direction (default grouping method 0)
+  run  1000000
+
+Using a specific grouping method::
+
+  ensemble nvt_ber 300 300 100
+  fix  1 0            # fix group 0 in grouping method 1
+  move 1 1 0.001 0 0  # move group 1 in grouping method 1, with a speed of 0.001 Ångstrom/fs in the x direction
   run  1000000
 
 Caveats
 -------
-* One cannot use NPT when using this keyword.
-* Currently, the moving group must be defined in the grouping method 0. This might be extended in a future version.
+* One cannot use NPT when using this keyword. Only ``nvt_ber``, ``nvt_nhc``, ``nvt_bdp``, and ``heat_lan`` ensembles are supported.
+* When both :ref:`fix <kw_fix>` and :attr:`move` are used, they must use the same grouping method.

src/integrate/ensemble.cu

@@ -299,7 +299,7 @@ void Ensemble::velocity_verlet_v()
   int n = atom->number_of_atoms;
   const int* group_pointer;
   if (group->size())
-    group_pointer = (*group)[0].label.data();
+    group_pointer = (*group)[fixed_grouping_method].label.data();
   else
     group_pointer = 0;
 
@@ -324,7 +324,7 @@ void Ensemble::velocity_verlet_x()
   int n = atom->number_of_atoms;
   const int* group_pointer;
   if (group->size())
-    group_pointer = (*group)[0].label.data();
+    group_pointer = (*group)[fixed_grouping_method].label.data();
   else
     group_pointer = 0;
   gpu_velocity_verlet_x<<<(n - 1) / 128 + 1, 128>>>(
@@ -380,7 +380,7 @@ void Ensemble::velocity_verlet(
       move_velocity[0],
       move_velocity[1],
       move_velocity[2],
-      group[0].label.data(),
+      group[fixed_grouping_method].label.data(),
       time_step,
       mass.data(),
       position_per_atom.data(),
@@ -646,10 +646,10 @@ void Ensemble::find_thermo(
   const int number_of_atoms = mass.size();
   int num_atoms_for_temperature = number_of_atoms;
   if (fixed_group >= 0) {
-    num_atoms_for_temperature -= group[0].cpu_size[fixed_group];
+    num_atoms_for_temperature -= group[fixed_grouping_method].cpu_size[fixed_group];
   }
   if (move_group >= 0) {
-    num_atoms_for_temperature -= group[0].cpu_size[move_group];
+    num_atoms_for_temperature -= group[move_grouping_method].cpu_size[move_group];
   }
 
   gpu_find_thermo_instant_temperature<<<8, 1024>>>(

src/integrate/ensemble.cuh

@@ -76,6 +76,8 @@ public:
   int sink;
   int fixed_group = -1; // ID of the group in which the atoms will be fixed
   int move_group = -1;  // ID of the group in which the atoms will move with a constant velocity
+  int fixed_grouping_method = 0;
+  int move_grouping_method = 0;
   double move_velocity[3];
   double temperature; // target temperature at a specific time
   double delta_temperature;

src/integrate/ensemble_bao.cu

@@ -278,7 +278,7 @@ void Ensemble_BAO::operator_A(
     gpu_operator_A<<<(number_of_atoms - 1) / 128 + 1, 128>>>(
       number_of_atoms,
       fixed_group,
-      group[0].label.data(),
+      group[fixed_grouping_method].label.data(),
       time_step,
       mass.data(),
       position_per_atom.data(),
@@ -400,7 +400,7 @@ void Ensemble_BAO::operator_B(
     gpu_operator_B<<<(number_of_atoms - 1) / 128 + 1, 128>>>(
       number_of_atoms,
       fixed_group,
-      group[0].label.data(),
+      group[fixed_grouping_method].label.data(),
       time_step,
       mass.data(),
       position_per_atom.data(),

src/integrate/ensemble_bdp.cu

@@ -84,8 +84,8 @@ void Ensemble_BDP::integrate_nvt_bdp_2(
     false, time_step, group, mass, force_per_atom, position_per_atom, velocity_per_atom);
 
   // get thermo
-  int N_fixed = (fixed_group == -1) ? 0 : group[0].cpu_size[fixed_group];
-  N_fixed += (move_group == -1) ? 0 : group[0].cpu_size[move_group];
+  int N_fixed = (fixed_group == -1) ? 0 : group[fixed_grouping_method].cpu_size[fixed_group];
+  N_fixed += (move_group == -1) ? 0 : group[move_grouping_method].cpu_size[move_group];
   find_thermo(
     true, volume, group, mass, potential_per_atom, velocity_per_atom, virial_per_atom, thermo);
 

src/integrate/ensemble_npt_scr.cu

@@ -239,7 +239,7 @@ void Ensemble_NPT_SCR::compute2(
     atom.position_per_atom,
     atom.velocity_per_atom);
 
-  int N_fixed = (fixed_group == -1) ? 0 : group[0].cpu_size[fixed_group];
+  int N_fixed = (fixed_group == -1) ? 0 : group[fixed_grouping_method].cpu_size[fixed_group];
   find_thermo(
     true,
     box.get_volume(),

src/integrate/integrate.cu

@@ -57,6 +57,9 @@ void Integrate::initialize(
     if (fixed_group < 0) {
       PRINT_INPUT_ERROR("It is not allowed to have moving group but no fixed group.");
     }
+    if (fixed_grouping_method != move_grouping_method) {
+      PRINT_INPUT_ERROR("The fixed and moving groups must use the same grouping method.");
+    }
     if (move_group == fixed_group) {
       PRINT_INPUT_ERROR("The fixed and moving groups cannot be the same.");
     }
@@ -208,12 +211,16 @@ void Integrate::initialize(
   ensemble->total_steps = &this->total_steps;
   ensemble->thermo = &thermo;
   ensemble->fixed_group = fixed_group;
+  ensemble->fixed_grouping_method = fixed_grouping_method;
+  ensemble->move_grouping_method = move_grouping_method;
 }
 
 void Integrate::finalize()
 {
   fixed_group = -1; // no group has an index of -1
   move_group = -1;
+  fixed_grouping_method = 0;
+  move_grouping_method = 0;
   deform_x = 0;
   deform_y = 0;
   deform_z = 0;
@@ -990,64 +997,105 @@ void Integrate::parse_ensemble(
 
 void Integrate::parse_fix(const char** param, int num_param, std::vector<Group>& group)
 {
-  if (num_param != 2) {
-    PRINT_INPUT_ERROR("Keyword 'fix' should have 1 parameter.");
-  }
-
-  if (!is_valid_int(param[1], &fixed_group)) {
-    PRINT_INPUT_ERROR("Fixed group ID should be an integer.");
+  if (num_param != 2 && num_param != 3) {
+    PRINT_INPUT_ERROR("Keyword 'fix' should have 1 or 2 parameters.");
   }
 
   if (group.size() < 1) {
     PRINT_INPUT_ERROR("Cannot use 'fix' without grouping method.");
   }
 
+  if (num_param == 3) {
+    // fix grouping_method group_id
+    if (!is_valid_int(param[1], &fixed_grouping_method)) {
+      PRINT_INPUT_ERROR("Grouping method for 'fix' should be an integer.");
+    }
+    if (fixed_grouping_method < 0) {
+      PRINT_INPUT_ERROR("Grouping method for 'fix' should >= 0.");
+    }
+    if (fixed_grouping_method >= group.size()) {
+      PRINT_INPUT_ERROR("Grouping method for 'fix' should < number of grouping methods.");
+    }
+    if (!is_valid_int(param[2], &fixed_group)) {
+      PRINT_INPUT_ERROR("Fixed group ID should be an integer.");
+    }
+  } else {
+    // fix group_id (default grouping_method = 0)
+    fixed_grouping_method = 0;
+    if (!is_valid_int(param[1], &fixed_group)) {
+      PRINT_INPUT_ERROR("Fixed group ID should be an integer.");
+    }
+  }
+
   if (fixed_group < 0) {
     PRINT_INPUT_ERROR("Fixed group ID should >= 0.");
   }
 
-  if (fixed_group >= group[0].number) {
+  if (fixed_group >= group[fixed_grouping_method].number) {
     PRINT_INPUT_ERROR("Fixed group ID should < number of groups.");
   }
 
-  printf("Group %d in grouping method 0 will be fixed.\n", fixed_group);
+  printf(
+    "Group %d in grouping method %d will be fixed.\n", fixed_group, fixed_grouping_method);
 }
 
 void Integrate::parse_move(const char** param, int num_param, std::vector<Group>& group)
 {
-  if (num_param != 5) {
-    PRINT_INPUT_ERROR("Keyword 'move' should have 4 parameters.");
-  }
-
-  if (!is_valid_int(param[1], &move_group)) {
-    PRINT_INPUT_ERROR("Moving group ID should be an integer.");
+  if (num_param != 5 && num_param != 6) {
+    PRINT_INPUT_ERROR("Keyword 'move' should have 4 or 5 parameters.");
   }
 
   if (group.size() < 1) {
     PRINT_INPUT_ERROR("Cannot use 'move' without grouping method.");
   }
 
+  int vid; // index where vx starts
+  if (num_param == 6) {
+    // move grouping_method group_id vx vy vz
+    if (!is_valid_int(param[1], &move_grouping_method)) {
+      PRINT_INPUT_ERROR("Grouping method for 'move' should be an integer.");
+    }
+    if (move_grouping_method < 0) {
+      PRINT_INPUT_ERROR("Grouping method for 'move' should >= 0.");
+    }
+    if (move_grouping_method >= group.size()) {
+      PRINT_INPUT_ERROR("Grouping method for 'move' should < number of grouping methods.");
+    }
+    if (!is_valid_int(param[2], &move_group)) {
+      PRINT_INPUT_ERROR("Moving group ID should be an integer.");
+    }
+    vid = 3;
+  } else {
+    // move group_id vx vy vz (default grouping_method = 0)
+    move_grouping_method = 0;
+    if (!is_valid_int(param[1], &move_group)) {
+      PRINT_INPUT_ERROR("Moving group ID should be an integer.");
+    }
+    vid = 2;
+  }
+
   if (move_group < 0) {
     PRINT_INPUT_ERROR("Moving group ID should >= 0.");
   }
 
-  if (move_group >= group[0].number) {
+  if (move_group >= group[move_grouping_method].number) {
     PRINT_INPUT_ERROR("Moving group ID should < number of groups.");
   }
 
-  if (!is_valid_real(param[2], &move_velocity[0])) {
+  if (!is_valid_real(param[vid], &move_velocity[0])) {
     PRINT_INPUT_ERROR("Moving velocity in x direction should be a number.");
   }
-  if (!is_valid_real(param[3], &move_velocity[1])) {
+  if (!is_valid_real(param[vid + 1], &move_velocity[1])) {
     PRINT_INPUT_ERROR("Moving velocity in y direction should be a number.");
   }
-  if (!is_valid_real(param[4], &move_velocity[2])) {
+  if (!is_valid_real(param[vid + 2], &move_velocity[2])) {
     PRINT_INPUT_ERROR("Moving velocity in z direction should be a number.");
   }
 
   printf(
-    "Group %d in grouping method 0 will move with velocity vector (%g, %g, %g) A/fs.\n",
+    "Group %d in grouping method %d will move with velocity vector (%g, %g, %g) A/fs.\n",
     move_group,
+    move_grouping_method,
     move_velocity[0],
     move_velocity[1],
     move_velocity[2]);

src/integrate/integrate.cuh

@@ -74,6 +74,8 @@ public:
   int sink;
   int fixed_group = -1; // ID of the group in which the atoms will be fixed
   int move_group = -1;  // ID of the group in which the atoms will move with a constant velocity
+  int fixed_grouping_method = 0;
+  int move_grouping_method = 0;
   double move_velocity[3];
 
   double temperature;  // target temperature at a specific time

src/main_gpumd/cohesive.cu

@@ -359,7 +359,7 @@ void Cohesive::compute(
     Box new_box;
     deform_box(num_atoms, cpu_D[n], box, new_box, position_per_atom, old_box_inv);
 
-    Minimizer_SD minimizer(-1, num_atoms, 1000, 1.0e-5);
+    Minimizer_SD minimizer(-1, 0, num_atoms, 1000, 1.0e-5);
     minimizer.compute(
       force,
       new_box,