Merge pull request #614 from Tingliangstu/master

allow compute SHC for all group IDs except group ID 0 in a given grouping method

doc/gpumd/input_parameters/compute_shc.rst

@@ -36,7 +36,8 @@ The angular frequency data will be :attr:`max_omega/num_omega, 2*max_omega/num_o
 
 This means that :math:`K(t)` will be calculated for atoms in group :attr:`group_id` of grouping method :attr:`grouping_method`.
 Usually, :attr:`group_id` should be :math:`\geq 0` and smaller than the number of groups in grouping method :attr:`grouping_method`.
-If :attr:`grouping_method` is 0 and :attr:`group_id` is -1, it means to calculate the :math:`K(t)` for every group.
+If :attr:`grouping_method` is assigned and :attr:`group_id` is -1, it means to calculate the :math:`K(t)` for every :attr:`group_id` except for :attr:`group_id` 0 in the assigned :attr:`grouping_method`.
+Since it is very time and memory consuming to calculate the all group :math:`K(t)` for a large system, so one can assign the part that don't want to calculate to :attr:`group_id` 0.
 Also, grouping method :attr:`grouping_method` must be defined in the :ref:`simulation model input file <model_xyz>`.
 If this option is missing, it means computing :math:`K(t)` for the whole system.
 
@@ -75,6 +76,22 @@ means that
 * you want to consider 500 frequency points
 * the maximum angular frequency is 200 THz
 
+Example 3
+^^^^^^^^^
+
+The command::
+
+  compute_shc 1 500 1 500 200.0 group 1 -1
+
+means that
+
+* you want to calculate :math:`K(t)` for all :attr:`group_id` except for :attr:`group_id` 0 defined in grouping method :attr:`1`
+* the sampling interval is 1 (sample the data at each time step)
+* the maximum number of correlation steps is 500
+* the transport direction is :math:`y`
+* you want to consider 500 frequency points
+* the maximum angular frequency is 200 THz
+
 Caveats
 -------
 This computation can be memory consuming.

doc/gpumd/output_files/shc_out.rst

@@ -32,3 +32,5 @@ In the next :attr:`num_omega` rows:
 :math:`J_q^{\rm in}(\omega) + J_q^{\rm out}(\omega) = J_q(\omega)` is exactly the left expression in Eq. (20) of [Fan2019]_.
 
 Only the potential part of the heat current has been included.
+
+If :attr: 'group_id' is -1, then the file follows the above rules and will contain :math:`K(t)` and :math:`J_q(\omega)` for each group id except for group id 0. And the contents of the :attr: 'group_id' are arranged from smallest to largest.

src/measure/shc.cu

@@ -41,7 +41,7 @@ void SHC::preprocess(const int N, const std::vector<Group>& group)
     group_size = N;
     group_num = 1;
   } else {
-    if (group_method == 0 && group_id == -1) {
+    if (group_id == -1) {
       group_size = N;
       group_num = group[group_method].number;
     } else {
@@ -185,8 +185,8 @@ void SHC::process(
     CHECK(cudaMemcpy(vy.data() + offset, vy_tmp, sizeof(double) * N, cudaMemcpyDeviceToDevice));
     CHECK(cudaMemcpy(vz.data() + offset, vz_tmp, sizeof(double) * N, cudaMemcpyDeviceToDevice));
   } else {
-    if (group_method == 0 && group_id == -1) {
-      for (int n = 0; n < group_num; ++n) {
+    if (group_id == -1) {
+      for (int n = 1; n < group_num; ++n) {
         int offset_s = Nc * group[group_method].cpu_size_sum[n] + 
                        correlation_step * group[group_method].cpu_size[n];
         gpu_copy_data<<<(group[group_method].cpu_size[n] - 1) / BLOCK_SIZE_SHC + 1, BLOCK_SIZE_SHC>>>(
@@ -231,8 +231,8 @@ void SHC::process(
   if (sample_step >= Nc - 1) {
     ++num_time_origins;
 
-    if (group_method == 0 && group_id == -1) {
-      for (int n = 0; n < group_num; ++n) {
+    if (group_id == -1) {
+      for (int n = 1; n < group_num; ++n) {
         int offset_s = Nc * group[group_method].cpu_size_sum[n];
         int offset_e = offset_s + correlation_step * group[group_method].cpu_size[n];
         gpu_find_k<<<Nc, BLOCK_SIZE_SHC>>>(
@@ -303,8 +303,8 @@ void SHC::average_k()
   ko_positive.copy_to_host(ko_positive_cpu.data());
 
   const double scalar = 1000.0 / TIME_UNIT_CONVERSION / num_time_origins;
-  if (group_method == 0 && group_id == -1) {
-    for (int n = 0; n < group_num; ++n) {
+  if (group_id == -1) {
+    for (int n = 1; n < group_num; ++n) {
       const int offset_k = (Nc * 2 - 1) * n;
       const int offset_n = Nc * n;
       for (int nc = 0; nc < Nc - 1; ++nc) {
@@ -331,8 +331,8 @@ void SHC::average_k()
 
 void SHC::find_shc(const double dt_in_ps, const double d_omega)
 {
-  if (group_method == 0 && group_id == -1) {
-    for (int n = 0; n < group_num; ++n) {
+  if (group_id == -1) {
+    for (int n = 1; n < group_num; ++n) {
       const int offset_k = (Nc * 2 - 1) * n;
       for (int nc = 0; nc < Nc * 2 - 1; ++nc) {
         const double hann_window = (cos(PI * (nc + 1 - Nc) / Nc) + 1.0) * 0.5;
@@ -341,7 +341,7 @@ void SHC::find_shc(const double dt_in_ps, const double d_omega)
       }
     }
 
-    for (int n = 0; n < group_num; ++n) {
+    for (int n = 1; n < group_num; ++n) {
       const int offset_s = num_omega * n;
       const int offset_k = (Nc * 2 - 1) * n;
       for (int nw = 0; nw < num_omega; ++nw) {
@@ -390,8 +390,8 @@ void SHC::postprocess(const double time_step)
 
   average_k();
   find_shc(dt_in_ps, d_omega);
-  if (group_method == 0 && group_id == -1) {
-    for (int n = 0; n < group_num; ++n) {
+  if (group_id == -1) {
+    for (int n = 1; n < group_num; ++n) {
       const int offset_k = (Nc * 2 - 1) * n;
       // ki and ko are in units of A*eV/ps
       for (int nc = 0; nc < Nc * 2 - 1; ++nc) {
@@ -487,4 +487,11 @@ void SHC::parse(const char** param, int num_param, const std::vector<Group>& gro
       PRINT_INPUT_ERROR("Unrecognized argument in compute_shc.\n");
     }
   }
+
+  if (group_id == -1) {
+    printf("    compute SHC for all group IDs except for group ID 0 in grouping method %d.\n", group_method);
+  }
+  if (group_id < -1) {
+    PRINT_INPUT_ERROR("group ID should >= -1 for computing SHC.");
+  }
 }