Refactoring of gridsys.c

c/bzgrid.c

@@ -176,10 +176,10 @@ static long get_ir_grid_map(long ir_mapping_table[],
                             const long PS[3],
                             const RotMats *rot_reciprocal);
 static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
-                                        const long (*grid_address)[3],
+                                        const long (*gr_grid_addresses)[3],
                                         const long Qinv[3][3]);
 static void get_bz_grid_addresses_type2(BZGrid *bzgrid,
-                                        const long (*grid_address)[3],
+                                        const long (*gr_grid_addresses)[3],
                                         const long Qinv[3][3]);
 static void set_bz_address(long address[3],
                            const long bz_index,
@@ -486,7 +486,7 @@ static long get_ir_grid_map(long ir_mapping_table[],
 }
 
 static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
-                                        const long (*grid_address)[3],
+                                        const long (*gr_grid_addresses)[3],
                                         const long Qinv[3][3])
 {
   double tolerance, min_distance;
@@ -513,8 +513,8 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
   bzgrid->gp_map[num_bzmesh] = 0;
   id_shift = 0;
   for (i = 0; i < total_num_gp; i++) {
-    min_distance = get_bz_distances(nint, distances, bzgrid, grid_address[i],
-                                    tolerance);
+    min_distance = get_bz_distances(nint, distances, bzgrid,
+                                    gr_grid_addresses[i], tolerance);
     count = 0;
     for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++) {
       if (distances[j] < min_distance + tolerance) {
@@ -527,7 +527,7 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
         count++;
         set_bz_address(bzgrid->addresses[gp],
                        j,
-                       grid_address[i],
+                       gr_grid_addresses[i],
                        bzgrid->D_diag,
                        nint,
                        Qinv);
@@ -548,7 +548,7 @@ static void get_bz_grid_addresses_type1(BZGrid *bzgrid,
 }
 
 static void get_bz_grid_addresses_type2(BZGrid *bzgrid,
-                                        const long (*grid_address)[3],
+                                        const long (*gr_grid_addresses)[3],
                                         const long Qinv[3][3])
 {
   double tolerance, min_distance;
@@ -563,13 +563,13 @@ static void get_bz_grid_addresses_type2(BZGrid *bzgrid,
 
   for (i = 0;
        i < bzgrid->D_diag[0] * bzgrid->D_diag[1] * bzgrid->D_diag[2]; i++) {
-    min_distance = get_bz_distances(nint, distances, bzgrid, grid_address[i],
-                                    tolerance);
+    min_distance = get_bz_distances(nint, distances, bzgrid,
+                                    gr_grid_addresses[i], tolerance);
     for (j = 0; j < BZG_NUM_BZ_SEARCH_SPACE; j++) {
       if (distances[j] < min_distance + tolerance) {
         set_bz_address(bzgrid->addresses[num_gp],
                        j,
-                       grid_address[i],
+                       gr_grid_addresses[i],
                        bzgrid->D_diag,
                        nint,
                        Qinv);

c/grgrid.c

@@ -64,7 +64,7 @@ static long rotate_grid_index(const long grid_index,
                               const long rotation[3][3],
                               const long D_diag[3],
                               const long PS[3]);
-static void get_ir_grid_map(long ir_grid_indices[],
+static void get_ir_grid_map(long *ir_grid_map,
                             const long (*rotations)[3][3],
                             const long num_rot,
                             const long D_diag[3],
@@ -237,13 +237,13 @@ long grg_rotate_grid_index(const long grid_index,
 /* -----------------------------*/
 /* Find irreducible grid points */
 /* -----------------------------*/
-void grg_get_ir_grid_map(long ir_grid_indices[],
+void grg_get_ir_grid_map(long *ir_grid_map,
                          const long (*rotations)[3][3],
                          const long num_rot,
                          const long D_diag[3],
                          const long PS[3])
 {
-  get_ir_grid_map(ir_grid_indices,
+  get_ir_grid_map(ir_grid_map,
                   rotations,
                   num_rot,
                   D_diag,
@@ -440,7 +440,10 @@ static long rotate_grid_index(const long grid_index,
   return get_double_grid_index(dadrs_rot, D_diag, PS);
 }
 
-static void get_ir_grid_map(long ir_grid_indices[],
+/* Find ir-grid points. */
+/* This algorithm relies on the ir-grid index is always smallest */
+/* number among symmetrically equivalent grid points. */
+static void get_ir_grid_map(long *ir_grid_map,
                             const long (*rotations)[3][3],
                             const long num_rot,
                             const long D_diag[3],
@@ -452,7 +455,7 @@ static void get_ir_grid_map(long ir_grid_indices[],
   num_gp = D_diag[0] * D_diag[1] * D_diag[2];
 
   for (gp = 0; gp < num_gp; gp++) {
-    ir_grid_indices[gp] = num_gp;
+    ir_grid_map[gp] = num_gp;
   }
 
   /* Do not simply multithreaded this for-loop. */
@@ -461,12 +464,12 @@ static void get_ir_grid_map(long ir_grid_indices[],
     for (i = 0; i < num_rot; i++) {
       r_gp = rotate_grid_index(gp, rotations[i], D_diag, PS);
       if (r_gp < gp) {
-        ir_grid_indices[gp] = ir_grid_indices[r_gp];
+        ir_grid_map[gp] = ir_grid_map[r_gp];
         break;
       }
     }
-    if (ir_grid_indices[gp] == num_gp) {
-      ir_grid_indices[gp] = gp;
+    if (ir_grid_map[gp] == num_gp) {
+      ir_grid_map[gp] = gp;
     }
   }
 

c/grgrid.h

@@ -67,7 +67,7 @@ long grg_rotate_grid_index(const long grid_index,
                            const long rotations[3][3],
                            const long D_diag[3],
                            const long PS[3]);
-void grg_get_ir_grid_map(long ir_grid_indices[],
+void grg_get_ir_grid_map(long *ir_grid_map,
                          const long (*rotations)[3][3],
                          const long num_rot,
                          const long D_diag[3],

c/gridsys.c

@@ -140,6 +140,8 @@ double gridsys_get_thm_integration_weight(const double omega,
 }
 
 
+/* Get one dataset of relative grid address used for tetrahedron */
+/* method. rec_lattice is used to choose the one. */
 /* rec_lattice : microzone basis vectors in column vectors */
 void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3],
                                            const double rec_lattice[3][3])
@@ -149,21 +151,81 @@ void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3
 
 
 /* The rotations are those after proper transformation in GRGrid. */
-void gridsys_get_ir_grid_map(long *ir_grid_indices,
+void gridsys_get_ir_grid_map(long *ir_grid_map,
                              const long (*rotations)[3][3],
                              const long num_rot,
                              const long D_diag[3],
                              const long PS[3])
 {
-  grg_get_ir_grid_map(ir_grid_indices, rotations, num_rot, D_diag, PS);
+  grg_get_ir_grid_map(ir_grid_map, rotations, num_rot, D_diag, PS);
 }
 
 
+/* Find shortest grid points from Gamma considering periodicity of */
+/* reciprocal lattice. See the details in docstring of BZGrid. */
+/* */
+/* Parameters */
+/* ---------- */
+/* type : Data structure type of bz_map. type=2 is always */
+/*        recommended.*/
+/* Returns */
+/* ------- */
+/* bz_grid_addresses : */
+/*     Grid point addresses of shortest grid points. */
+/*     len(bz_grid_addresses) <= product(D_diag + 1) */
+/*     for full gr_grid_addresses. */
+/* bz_map : */
+/*     List of accumulated numbers of BZ grid points from the first GR */
+/*     grid point to the last grid point. */
+/*     [0, 1, 3, 4, ...] means multiplicities of [1, 2, 1, ...] */
+/*     With type=2, len(bz_map)=product(D_diag) + 1. */
+/* bzg2grg : */
+/*     Mapping table of bz_grid_addresses to gr_grid_addresses. */
+/*     len(bzg2grg) == len(bz_grid_addresses) <= product(D_diag + 1). */
+/* (function return) size : */
+/*     Number of bz_grid_addresses stored. */
+long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
+                                   long *bz_map,
+                                   long *bzg2grg,
+                                   const long (*gr_grid_addresses)[3],
+                                   const long D_diag[3],
+                                   const long Q[3][3],
+                                   const long PS[3],
+                                   const double rec_lattice[3][3],
+                                   const long type)
+{
+  BZGrid *bzgrid;
+  long i, j, size;
 
+  if ((bzgrid = (BZGrid*) malloc(sizeof(BZGrid))) == NULL) {
+    warning_print("Memory could not be allocated.");
+    return 0;
+  }
 
+  bzgrid->addresses = bz_grid_addresses;
+  bzgrid->gp_map = bz_map;
+  bzgrid->bzg2grg = bzg2grg;
+  bzgrid->type = type;
+  for (i = 0; i < 3; i++) {
+    bzgrid->D_diag[i] = D_diag[i];
+    bzgrid->PS[i] = PS[i];
+    for (j = 0; j < 3; j++) {
+      bzgrid->Q[i][j] = Q[i][j];
+      bzgrid->reclat[i][j] = rec_lattice[i][j];
+    }
+  }
 
+  if (bzg_get_bz_grid_addresses(bzgrid, gr_grid_addresses)) {
+    size = bzgrid->size;
+  } else {
+    size = 0;
+  }
 
+  free(bzgrid);
+  bzgrid = NULL;
 
+  return size;
+}
 
 
 long gridsys_get_triplets_at_q(long *map_triplets,
@@ -226,9 +288,9 @@ long gridsys_get_BZ_triplets_at_q(long (*triplets)[3],
   return num_ir;
 }
 
-/* relative_grid_addresses are given as P multipled with those from dataset,
- * i.e.,
- *     np.dot(relative_grid_addresses, P.T) */
+/* relative_grid_addresses are given as P multipled with those from */
+/* dataset, i.e., */
+/*     np.dot(relative_grid_addresses, P.T) */
 long gridsys_get_integration_weight(double *iw,
                                     char *iw_zero,
                                     const double *frequency_points,
@@ -308,47 +370,3 @@ void gridsys_get_integration_weight_with_sigma(double *iw,
                                         num_band,
                                         tp_type);
 }
-
-
-long gridsys_get_bz_grid_address(long (*bz_grid_addresses)[3],
-                                 long *bz_map,
-                                 long *bzg2grg,
-                                 const long (*grid_address)[3],
-                                 const long D_diag[3],
-                                 const long Q[3][3],
-                                 const long PS[3],
-                                 const double rec_lattice[3][3],
-                                 const long type)
-{
-  BZGrid *bzgrid;
-  long i, j, size;
-
-  if ((bzgrid = (BZGrid*) malloc(sizeof(BZGrid))) == NULL) {
-    warning_print("Memory could not be allocated.");
-    return 0;
-  }
-
-  bzgrid->addresses = bz_grid_addresses;
-  bzgrid->gp_map = bz_map;
-  bzgrid->bzg2grg = bzg2grg;
-  bzgrid->type = type;
-  for (i = 0; i < 3; i++) {
-    bzgrid->D_diag[i] = D_diag[i];
-    bzgrid->PS[i] = PS[i];
-    for (j = 0; j < 3; j++) {
-      bzgrid->Q[i][j] = Q[i][j];
-      bzgrid->reclat[i][j] = rec_lattice[i][j];
-    }
-  }
-
-  if (bzg_get_bz_grid_addresses(bzgrid, grid_address)) {
-    size = bzgrid->size;
-  } else {
-    size = 0;
-  }
-
-  free(bzgrid);
-  bzgrid = NULL;
-
-  return size;
-}

c/gridsys.h

@@ -75,14 +75,20 @@ double gridsys_get_thm_integration_weight(const double omega,
                                           const char function);
 void gridsys_get_thm_relative_grid_address(long relative_grid_addresses[24][4][3],
                                            const double rec_lattice[3][3]);
-void gridsys_get_ir_grid_map(long *ir_grid_indices,
+void gridsys_get_ir_grid_map(long *ir_grid_map,
                              const long (*rotations)[3][3],
                              const long num_rot,
                              const long D_diag[3],
                              const long PS[3]);
-
-
-
+long gridsys_get_bz_grid_addresses(long (*bz_grid_addresses)[3],
+                                   long *bz_map,
+                                   long *bzg2grg,
+                                   const long (*gr_grid_addresses)[3],
+                                   const long D_diag[3],
+                                   const long Q[3][3],
+                                   const long PS[3],
+                                   const double rec_lattice[3][3],
+                                   const long type);
 long gridsys_get_triplets_at_q(long *map_triplets,
                                long *map_q,
                                const long grid_point,
@@ -129,15 +135,6 @@ void gridsys_get_integration_weight_with_sigma(double *iw,
                                                const double *frequencies,
                                                const long num_band,
                                                const long tp_type);
-long gridsys_get_bz_grid_address(long (*bz_grid_addresses)[3],
-                                 long *bz_map,
-                                 long *bzg2grg,
-                                 const long (*grid_address)[3],
-                                 const long D_diag[3],
-                                 const long Q[3][3],
-                                 const long PS[3],
-                                 const double rec_lattice[3][3],
-                                 const long type);
 
 
 #ifdef __cplusplus