int -> long in _phono3py.c under way

c/_phono3py.c

@@ -91,12 +91,12 @@ static void pinv_from_eigensolution(double *data,
                                     const double *eigvals,
                                     const long size,
                                     const double cutoff,
-                                    const int pinv_method);
+                                    const long pinv_method);
 static void show_colmat_info(const PyArrayObject *collision_matrix_py,
                              const long i_sigma,
                              const long i_temp,
                              const long adrs_shift);
-static Iarray* convert_to_iarray(const PyArrayObject* npyary);
+static Larray* convert_to_larray(const PyArrayObject* npyary);
 static Darray* convert_to_darray(const PyArrayObject* npyary);
 
 
@@ -306,7 +306,7 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
   PyArrayObject *py_s2p_map;
   PyArrayObject *py_band_indices;
   double cutoff_frequency;
-  int symmetrize_fc3_q;
+  long symmetrize_fc3_q;
 
   Darray *fc3_normal_squared;
   Darray *freqs;
@@ -314,20 +314,20 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
   long (*triplets)[3];
   long num_triplets;
   char* g_zero;
-  int *grid_address;
+  long *grid_address;
-  int *mesh;
+  long *mesh;
   double *fc3;
   double *svecs;
-  int *multi;
+  long *multi;
   double *masses;
-  int *p2s;
+  long *p2s;
-  int *s2p;
+  long *s2p;
-  int *band_indices;
+  long *band_indices;
-  int svecs_dims[3];
+  long svecs_dims[3];
-  int i;
+  long i;
-  int is_compact_fc3;
+  long is_compact_fc3;
 
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOid",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOld",
                         &py_fc3_normal_squared,
                         &py_g_zero,
                         &py_frequencies,
@@ -356,8 +356,8 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   num_triplets = (long)PyArray_DIMS(py_triplets)[0];
   g_zero = (char*)PyArray_DATA(py_g_zero);
-  grid_address = (int*)PyArray_DATA(py_grid_address);
+  grid_address = (long*)PyArray_DATA(py_grid_address);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
   fc3 = (double*)PyArray_DATA(py_fc3);
   if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
     is_compact_fc3 = 0;
@@ -368,11 +368,11 @@ static PyObject * py_get_interaction(PyObject *self, PyObject *args)
   for (i = 0; i < 3; i++) {
     svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
   }
-  multi = (int*)PyArray_DATA(py_multiplicities);
+  multi = (long*)PyArray_DATA(py_multiplicities);
   masses = (double*)PyArray_DATA(py_masses);
-  p2s = (int*)PyArray_DATA(py_p2s_map);
+  p2s = (long*)PyArray_DATA(py_p2s_map);
-  s2p = (int*)PyArray_DATA(py_s2p_map);
+  s2p = (long*)PyArray_DATA(py_s2p_map);
-  band_indices = (int*)PyArray_DATA(py_band_indices);
+  band_indices = (long*)PyArray_DATA(py_band_indices);
 
   ph3py_get_interaction(fc3_normal_squared,
                         g_zero,
@@ -422,32 +422,32 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
   PyArrayObject *py_band_indices;
   PyArrayObject *py_temperatures;
   double cutoff_frequency;
-  int is_NU;
+  long is_NU;
-  int symmetrize_fc3_q;
+  long symmetrize_fc3_q;
 
   double *gamma;
-  int (*relative_grid_address)[4][3];
+  long (*relative_grid_address)[4][3];
   double *frequencies;
   lapack_complex_double *eigenvectors;
   long (*triplets)[3];
   long num_triplets;
   long *triplet_weights;
-  int *grid_address;
+  long *grid_address;
   long *bz_map;
-  int *mesh;
+  long *mesh;
   double *fc3;
   double *svecs;
-  int *multi;
+  long *multi;
   double *masses;
-  int *p2s;
+  long *p2s;
-  int *s2p;
+  long *s2p;
-  Iarray *band_indices;
+  Larray *band_indices;
   Darray *temperatures;
-  int svecs_dims[3];
+  long svecs_dims[3];
-  int i;
+  long i;
-  int is_compact_fc3;
+  long is_compact_fc3;
 
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOiid",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOOOOOOOOlld",
                         &py_gamma,
                         &py_relative_grid_address,
                         &py_frequencies,
@@ -472,15 +472,15 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
   }
 
   gamma = (double*)PyArray_DATA(py_gamma);
-  relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
+  relative_grid_address = (long(*)[4][3])PyArray_DATA(py_relative_grid_address);
   frequencies = (double*)PyArray_DATA(py_frequencies);
   eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   num_triplets = (long)PyArray_DIMS(py_triplets)[0];
   triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
-  grid_address = (int*)PyArray_DATA(py_grid_address);
+  grid_address = (long*)PyArray_DATA(py_grid_address);
   bz_map = (long*)PyArray_DATA(py_bz_map);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
   fc3 = (double*)PyArray_DATA(py_fc3);
   if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
     is_compact_fc3 = 0;
@@ -491,11 +491,11 @@ static PyObject * py_get_pp_collision(PyObject *self, PyObject *args)
   for (i = 0; i < 3; i++) {
     svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
   }
-  multi = (int*)PyArray_DATA(py_multiplicities);
+  multi = (long*)PyArray_DATA(py_multiplicities);
   masses = (double*)PyArray_DATA(py_masses);
-  p2s = (int*)PyArray_DATA(py_p2s_map);
+  p2s = (long*)PyArray_DATA(py_p2s_map);
-  s2p = (int*)PyArray_DATA(py_s2p_map);
+  s2p = (long*)PyArray_DATA(py_s2p_map);
-  band_indices = convert_to_iarray(py_band_indices);
+  band_indices = convert_to_larray(py_band_indices);
   temperatures = convert_to_darray(py_temperatures);
 
   ph3py_get_pp_collision(gamma,
@@ -547,8 +547,8 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
   PyArrayObject *py_s2p_map;
   PyArrayObject *py_band_indices;
   PyArrayObject *py_temperatures;
-  int is_NU;
+  long is_NU;
-  int symmetrize_fc3_q;
+  long symmetrize_fc3_q;
   double sigma;
   double sigma_cutoff;
   double cutoff_frequency;
@@ -559,21 +559,21 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
   long (*triplets)[3];
   long num_triplets;
   long *triplet_weights;
-  int *grid_address;
+  long *grid_address;
-  int *mesh;
+  long *mesh;
   double *fc3;
   double *svecs;
-  int *multi;
+  long *multi;
   double *masses;
-  int *p2s;
+  long *p2s;
-  int *s2p;
+  long *s2p;
-  Iarray *band_indices;
+  Larray *band_indices;
   Darray *temperatures;
-  int svecs_dims[3];
+  long svecs_dims[3];
-  int i;
+  long i;
-  int is_compact_fc3;
+  long is_compact_fc3;
 
-  if (!PyArg_ParseTuple(args, "OddOOOOOOOOOOOOOOiid",
+  if (!PyArg_ParseTuple(args, "OddOOOOOOOOOOOOOOlld",
                         &py_gamma,
                         &sigma,
                         &sigma_cutoff,
@@ -603,8 +603,8 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   num_triplets = (long)PyArray_DIMS(py_triplets)[0];
   triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
-  grid_address = (int*)PyArray_DATA(py_grid_address);
+  grid_address = (long*)PyArray_DATA(py_grid_address);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
   fc3 = (double*)PyArray_DATA(py_fc3);
   if (PyArray_DIMS(py_fc3)[0] == PyArray_DIMS(py_fc3)[1]) {
     is_compact_fc3 = 0;
@@ -615,11 +615,11 @@ static PyObject * py_get_pp_collision_with_sigma(PyObject *self, PyObject *args)
   for (i = 0; i < 3; i++) {
     svecs_dims[i] = PyArray_DIMS(py_shortest_vectors)[i];
   }
-  multi = (int*)PyArray_DATA(py_multiplicities);
+  multi = (long*)PyArray_DATA(py_multiplicities);
   masses = (double*)PyArray_DATA(py_masses);
-  p2s = (int*)PyArray_DATA(py_p2s_map);
+  p2s = (long*)PyArray_DATA(py_p2s_map);
-  s2p = (int*)PyArray_DATA(py_s2p_map);
+  s2p = (long*)PyArray_DATA(py_s2p_map);
-  band_indices = convert_to_iarray(py_band_indices);
+  band_indices = convert_to_larray(py_band_indices);
   temperatures = convert_to_darray(py_temperatures);
 
   ph3py_get_pp_collision_with_sigma(gamma,
@@ -664,7 +664,7 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
   PyArrayObject *py_g;
   PyArrayObject *py_g_zero;
   double cutoff_frequency, temperature;
-  int frequency_point_index;
+  long frequency_point_index;
 
   Darray *fc3_normal_squared;
   double *gamma;
@@ -675,7 +675,7 @@ static PyObject * py_get_imag_self_energy_with_g(PyObject *self, PyObject *args)
   long *triplet_weights;
   long num_frequency_points;
 
-  if (!PyArg_ParseTuple(args, "OOOOOdOOdi",
+  if (!PyArg_ParseTuple(args, "OOOOOdOOdl",
                         &py_gamma,
                         &py_fc3_normal_squared,
                         &py_triplets,
@@ -740,7 +740,7 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
   double *frequencies;
   long (*triplets)[3];
   long *triplet_weights;
-  int *grid_address;
+  long *grid_address;
 
   if (!PyArg_ParseTuple(args, "OOOOOOOOdOOd",
                         &py_gamma_detail,
@@ -767,7 +767,7 @@ py_get_detailed_imag_self_energy_with_g(PyObject *self, PyObject *args)
   frequencies = (double*)PyArray_DATA(py_frequencies);
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
-  grid_address = (int*)PyArray_DATA(py_grid_address);
+  grid_address = (long*)PyArray_DATA(py_grid_address);
 
   ph3py_get_detailed_imag_self_energy_at_bands_with_g(gamma_detail,
                                                       gamma_N,
@@ -802,7 +802,7 @@ static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
   Darray *fc3_normal_squared;
   double *shift;
   double *frequencies;
-  int *band_indices;
+  long *band_indices;
   long (*triplets)[3];
   long *triplet_weights;
 
@@ -824,7 +824,7 @@ static PyObject * py_get_real_self_energy_at_bands(PyObject *self,
   fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
   shift = (double*)PyArray_DATA(py_shift);
   frequencies = (double*)PyArray_DATA(py_frequencies);
-  band_indices = (int*)PyArray_DATA(py_band_indices);
+  band_indices = (long*)PyArray_DATA(py_band_indices);
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
 
@@ -860,7 +860,7 @@ static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
   Darray *fc3_normal_squared;
   double *shift;
   double *frequencies;
-  int *band_indices;
+  long *band_indices;
   long (*triplets)[3];
   long *triplet_weights;
 
@@ -883,7 +883,7 @@ static PyObject * py_get_real_self_energy_at_frequency_point(PyObject *self,
   fc3_normal_squared = convert_to_darray(py_fc3_normal_squared);
   shift = (double*)PyArray_DATA(py_shift);
   frequencies = (double*)PyArray_DATA(py_frequencies);
-  band_indices = (int*)PyArray_DATA(py_band_indices);
+  band_indices = (long*)PyArray_DATA(py_band_indices);
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   triplet_weights = (long*)PyArray_DATA(py_triplet_weights);
 
@@ -1058,10 +1058,10 @@ static PyObject * py_symmetrize_collision_matrix(PyObject *self, PyObject *args)
   }
 
   collision_matrix = (double*)PyArray_DATA(py_collision_matrix);
-  num_sigma = PyArray_DIMS(py_collision_matrix)[0];
+  num_sigma = (long)PyArray_DIMS(py_collision_matrix)[0];
-  num_temp = PyArray_DIMS(py_collision_matrix)[1];
+  num_temp = (long)PyArray_DIMS(py_collision_matrix)[1];
-  num_grid_points = PyArray_DIMS(py_collision_matrix)[2];
+  num_grid_points = (long)PyArray_DIMS(py_collision_matrix)[2];
-  num_band = PyArray_DIMS(py_collision_matrix)[3];
+  num_band = (long)PyArray_DIMS(py_collision_matrix)[3];
 
   if (PyArray_NDIM(py_collision_matrix) == 8) {
     num_column = num_grid_points * num_band * 3;
@@ -1133,7 +1133,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
   double *gamma;
   double *frequencies;
   lapack_complex_double *eigenvectors;
-  int *band_indices;
+  long *band_indices;
   double *mass_variances;
   long num_band, num_band0;
 
@@ -1154,7 +1154,7 @@ static PyObject * py_get_isotope_strength(PyObject *self, PyObject *args)
   gamma = (double*)PyArray_DATA(py_gamma);
   frequencies = (double*)PyArray_DATA(py_frequencies);
   eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
-  band_indices = (int*)PyArray_DATA(py_band_indices);
+  band_indices = (long*)PyArray_DATA(py_band_indices);
   mass_variances = (double*)PyArray_DATA(py_mass_variances);
   num_band = (long)PyArray_DIMS(py_frequencies)[1];
   num_band0 = (long)PyArray_DIMS(py_band_indices)[0];
@@ -1192,7 +1192,7 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
   long *ir_grid_points;
   long *weights;
   lapack_complex_double *eigenvectors;
-  int *band_indices;
+  long *band_indices;
   double *mass_variances;
   long num_band, num_band0, num_ir_grid_points;
   double *integration_weights;
@@ -1217,7 +1217,7 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
   ir_grid_points = (long*)PyArray_DATA(py_ir_grid_points);
   weights = (long*)PyArray_DATA(py_weights);
   eigenvectors = (lapack_complex_double*)PyArray_DATA(py_eigenvectors);
-  band_indices = (int*)PyArray_DATA(py_band_indices);
+  band_indices = (long*)PyArray_DATA(py_band_indices);
   mass_variances = (double*)PyArray_DATA(py_mass_variances);
   num_band = (long)PyArray_DIMS(py_frequencies)[1];
   num_band0 = (long)PyArray_DIMS(py_band_indices)[0];
@@ -1244,17 +1244,17 @@ static PyObject * py_get_thm_isotope_strength(PyObject *self, PyObject *args)
 static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
 {
   PyArrayObject *force_constants_third;
-  int target;
+  long target;
-  int source;
+  long source;
   PyArrayObject *rotation_cart_inv;
   PyArrayObject *atom_mapping_py;
 
   double *fc3;
   double *rot_cart_inv;
-  int *atom_mapping;
+  long *atom_mapping;
   long num_atom;
 
-  if (!PyArg_ParseTuple(args, "OiiOO",
+  if (!PyArg_ParseTuple(args, "OllOO",
                         &force_constants_third,
                         &target,
                         &source,
@@ -1265,7 +1265,7 @@ static PyObject * py_distribute_fc3(PyObject *self, PyObject *args)
 
   fc3 = (double*)PyArray_DATA(force_constants_third);
   rot_cart_inv = (double*)PyArray_DATA(rotation_cart_inv);
-  atom_mapping = (int*)PyArray_DATA(atom_mapping_py);
+  atom_mapping = (long*)PyArray_DATA(atom_mapping_py);
   num_atom = (long)PyArray_DIMS(atom_mapping_py)[0];
 
   ph3py_distribute_fc3(fc3,
@@ -1290,7 +1290,7 @@ static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
   double (*delta_fc2s)[3][3];
   double *inv_U;
   double (*site_sym_cart)[3][3];
-  int *rot_map_syms;
+  long *rot_map_syms;
   long num_atom, num_disp, num_site_sym;
 
   if (!PyArg_ParseTuple(args, "OOOOO",
@@ -1311,7 +1311,7 @@ static PyObject * py_rotate_delta_fc2s(PyObject *self, PyObject *args)
   /* (n_sym, 3, 3) */
   site_sym_cart = (double(*)[3][3])PyArray_DATA(py_site_sym_cart);
   /* (n_sym, natom) */
-  rot_map_syms = (int*)PyArray_DATA(py_rot_map_syms);
+  rot_map_syms = (long*)PyArray_DATA(py_rot_map_syms);
 
   num_atom = (long)PyArray_DIMS(py_fc3)[0];
   num_disp = (long)PyArray_DIMS(py_delta_fc2s)[0];
@@ -1359,10 +1359,10 @@ py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
   PyArrayObject* py_nsym_list;
 
   double *fc3;
-  int *s2pp;
+  long *s2pp;
-  int *p2s;
+  long *p2s;
-  int *nsym_list;
+  long *nsym_list;
-  int *perms;
+  long *perms;
   long n_patom, n_satom;
 
   if (!PyArg_ParseTuple(args, "OOOOO",
@@ -1375,10 +1375,10 @@ py_set_permutation_symmetry_compact_fc3(PyObject *self, PyObject *args)
   }
 
   fc3 = (double*)PyArray_DATA(py_fc3);
-  perms = (int*)PyArray_DATA(py_permutations);
+  perms = (long*)PyArray_DATA(py_permutations);
-  s2pp = (int*)PyArray_DATA(py_s2pp_map);
+  s2pp = (long*)PyArray_DATA(py_s2pp_map);
-  p2s = (int*)PyArray_DATA(py_p2s_map);
+  p2s = (long*)PyArray_DATA(py_p2s_map);
-  nsym_list = (int*)PyArray_DATA(py_nsym_list);
+  nsym_list = (long*)PyArray_DATA(py_nsym_list);
   n_patom = (long)PyArray_DIMS(py_fc3)[0];
   n_satom = (long)PyArray_DIMS(py_fc3)[1];
 
@@ -1400,16 +1400,16 @@ static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
   PyArrayObject* py_s2pp_map;
   PyArrayObject* py_p2s_map;
   PyArrayObject* py_nsym_list;
-  int t_type;
+  long t_type;
 
   double *fc3;
-  int *s2pp;
+  long *s2pp;
-  int *p2s;
+  long *p2s;
-  int *nsym_list;
+  long *nsym_list;
-  int *perms;
+  long *perms;
   long n_patom, n_satom;
 
-  if (!PyArg_ParseTuple(args, "OOOOOi",
+  if (!PyArg_ParseTuple(args, "OOOOOl",
                         &py_fc3,
                         &py_permutations,
                         &py_s2pp_map,
@@ -1420,10 +1420,10 @@ static PyObject * py_transpose_compact_fc3(PyObject *self, PyObject *args)
   }
 
   fc3 = (double*)PyArray_DATA(py_fc3);
-  perms = (int*)PyArray_DATA(py_permutations);
+  perms = (long*)PyArray_DATA(py_permutations);
-  s2pp = (int*)PyArray_DATA(py_s2pp_map);
+  s2pp = (long*)PyArray_DATA(py_s2pp_map);
-  p2s = (int*)PyArray_DATA(py_p2s_map);
+  p2s = (long*)PyArray_DATA(py_p2s_map);
-  nsym_list = (int*)PyArray_DATA(py_nsym_list);
+  nsym_list = (long*)PyArray_DATA(py_nsym_list);
   n_patom = (long)PyArray_DIMS(py_fc3)[0];
   n_satom = (long)PyArray_DIMS(py_fc3)[1];
 
@@ -1451,9 +1451,9 @@ static PyObject * py_get_neighboring_grid_points(PyObject *self, PyObject *args)
   long *relative_grid_points;
   long *grid_points;
   long num_grid_points, num_relative_grid_address;
-  int (*relative_grid_address)[3];
+  long (*relative_grid_address)[3];
-  int *mesh;
+  long *mesh;
-  int (*bz_grid_address)[3];
+  long (*bz_grid_address)[3];
   long *bz_map;
 
   if (!PyArg_ParseTuple(args, "OOOOOO",
@@ -1469,10 +1469,10 @@ static PyObject * py_get_neighboring_grid_points(PyObject *self, PyObject *args)
   relative_grid_points = (long*)PyArray_DATA(py_relative_grid_points);
   grid_points = (long*)PyArray_DATA(py_grid_points);
   num_grid_points = (long)PyArray_DIMS(py_grid_points)[0];
-  relative_grid_address = (int(*)[3])PyArray_DATA(py_relative_grid_address);
+  relative_grid_address = (long(*)[3])PyArray_DATA(py_relative_grid_address);
   num_relative_grid_address = (long)PyArray_DIMS(py_relative_grid_address)[0];
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
-  bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
+  bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address);
   bz_map = (long*)PyArray_DATA(py_bz_map);
 
   ph3py_get_neighboring_gird_points(relative_grid_points,
@@ -1501,10 +1501,10 @@ static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
   double *iw;
   double *frequency_points;
   long num_band0, num_band, num_gp;
-  int (*relative_grid_address)[4][3];
+  long (*relative_grid_address)[4][3];
-  int *mesh;
+  long *mesh;
   long *grid_points;
-  int (*bz_grid_address)[3];
+  long (*bz_grid_address)[3];
   long *bz_map;
   double *frequencies;
 
@@ -1523,11 +1523,11 @@ static PyObject * py_set_integration_weights(PyObject *self, PyObject *args)
   iw = (double*)PyArray_DATA(py_iw);
   frequency_points = (double*)PyArray_DATA(py_frequency_points);
   num_band0 = (long)PyArray_DIMS(py_frequency_points)[0];
-  relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
+  relative_grid_address = (long(*)[4][3])PyArray_DATA(py_relative_grid_address);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
   grid_points = (long*)PyArray_DATA(py_grid_points);
   num_gp = (long)PyArray_DIMS(py_grid_points)[0];
-  bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
+  bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address);
   bz_map = (long*)PyArray_DATA(py_bz_map);
   frequencies = (double*)PyArray_DATA(py_frequencies);
   num_band = (long)PyArray_DIMS(py_frequencies)[1];
@@ -1556,18 +1556,18 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
   PyArrayObject *py_mesh;
   PyArrayObject *py_rotations;
   long fixed_grid_number;
-  int is_time_reversal;
+  long is_time_reversal;
-  int swappable;
+  long swappable;
 
-  int (*grid_address)[3];
+  long (*grid_address)[3];
   long *map_triplets;
   long *map_q;
-  int *mesh;
+  long *mesh;
-  int (*rot)[3][3];
+  long (*rot)[3][3];
   long num_rot;
   long num_ir;
 
-  if (!PyArg_ParseTuple(args, "OOOlOiOi",
+  if (!PyArg_ParseTuple(args, "OOOlOlOl",
                         &py_map_triplets,
                         &py_map_q,
                         &py_grid_address,
@@ -1579,11 +1579,11 @@ py_tpl_get_triplets_reciprocal_mesh_at_q(PyObject *self, PyObject *args)
     return NULL;
   }
 
-  grid_address = (int(*)[3])PyArray_DATA(py_grid_address);
+  grid_address = (long(*)[3])PyArray_DATA(py_grid_address);
   map_triplets = (long*)PyArray_DATA(py_map_triplets);
   map_q = (long*)PyArray_DATA(py_map_q);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
-  rot = (int(*)[3][3])PyArray_DATA(py_rotations);
+  rot = (long(*)[3][3])PyArray_DATA(py_rotations);
   num_rot = (long)PyArray_DIMS(py_rotations)[0];
   num_ir = ph3py_get_triplets_reciprocal_mesh_at_q(map_triplets,
                                                    map_q,
@@ -1608,11 +1608,11 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
   long grid_point;
 
   long (*triplets)[3];
-  int (*bz_grid_address)[3];
+  long (*bz_grid_address)[3];
   long *bz_map;
   long *map_triplets;
   long num_map_triplets;
-  int *mesh;
+  long *mesh;
   long num_ir;
 
   if (!PyArg_ParseTuple(args, "OlOOOO",
@@ -1626,11 +1626,11 @@ static PyObject * py_tpl_get_BZ_triplets_at_q(PyObject *self, PyObject *args)
   }
 
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
-  bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
+  bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address);
   bz_map = (long*)PyArray_DATA(py_bz_map);
   map_triplets = (long*)PyArray_DATA(py_map_triplets);
   num_map_triplets = (long)PyArray_DIMS(py_map_triplets)[0];
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
 
   num_ir = ph3py_get_BZ_triplets_at_q(triplets,
                                       grid_point,
@@ -1656,20 +1656,20 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
   PyArrayObject *py_frequencies2;
   PyArrayObject *py_bz_grid_address;
   PyArrayObject *py_bz_map;
-  int tp_type;
+  long tp_type;
 
   double *iw;
   char *iw_zero;
   double *frequency_points;
-  int (*relative_grid_address)[4][3];
+  long (*relative_grid_address)[4][3];
-  int *mesh;
+  long *mesh;
   long (*triplets)[3];
-  int (*bz_grid_address)[3];
+  long (*bz_grid_address)[3];
   long *bz_map;
   double *frequencies1, *frequencies2;
   long num_band0, num_band1, num_band2, num_triplets;
 
-  if (!PyArg_ParseTuple(args, "OOOOOOOOOOi",
+  if (!PyArg_ParseTuple(args, "OOOOOOOOOOl",
                         &py_iw,
                         &py_iw_zero,
                         &py_frequency_points,
@@ -1688,11 +1688,11 @@ py_set_triplets_integration_weights(PyObject *self, PyObject *args)
   iw_zero = (char*)PyArray_DATA(py_iw_zero);
   frequency_points = (double*)PyArray_DATA(py_frequency_points);
   num_band0 = (long)PyArray_DIMS(py_frequency_points)[0];
-  relative_grid_address = (int(*)[4][3])PyArray_DATA(py_relative_grid_address);
+  relative_grid_address = (long(*)[4][3])PyArray_DATA(py_relative_grid_address);
-  mesh = (int*)PyArray_DATA(py_mesh);
+  mesh = (long*)PyArray_DATA(py_mesh);
   triplets = (long(*)[3])PyArray_DATA(py_triplets);
   num_triplets = (long)PyArray_DIMS(py_triplets)[0];
-  bz_grid_address = (int(*)[3])PyArray_DATA(py_bz_grid_address);
+  bz_grid_address = (long(*)[3])PyArray_DATA(py_bz_grid_address);
   bz_map = (long*)PyArray_DATA(py_bz_map);
   frequencies1 = (double*)PyArray_DATA(py_frequencies1);
   frequencies2 = (double*)PyArray_DATA(py_frequencies2);
@@ -1779,15 +1779,15 @@ py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
   PyArrayObject *py_collision_matrix;
   PyArrayObject *py_eigenvalues;
   double cutoff;
-  int i_sigma, i_temp, is_pinv, solver;
+  long i_sigma, i_temp, is_pinv, solver;
 
   double *collision_matrix;
   double *eigvals;
   long num_temp, num_grid_point, num_band;
   long num_column, adrs_shift;
-  int info;
+  long info;
 
-  if (!PyArg_ParseTuple(args, "OOiidii",
+  if (!PyArg_ParseTuple(args, "OOlldll",
                         &py_collision_matrix,
                         &py_eigenvalues,
                         &i_sigma,
@@ -1826,7 +1826,7 @@ py_diagonalize_collision_matrix(PyObject *self, PyObject *args)
                             eigvals, num_column, cutoff, 0);
   }
 
-  return PyLong_FromLong((long) info);
+  return PyLong_FromLong(info);
 }
 
 static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
@@ -1834,14 +1834,14 @@ static PyObject * py_pinv_from_eigensolution(PyObject *self, PyObject *args)
   PyArrayObject *py_collision_matrix;
   PyArrayObject *py_eigenvalues;
   double cutoff;
-  int i_sigma, i_temp, pinv_method;
+  long i_sigma, i_temp, pinv_method;
 
   double *collision_matrix;
   double *eigvals;
   long num_temp, num_grid_point, num_band;
   long num_column, adrs_shift;
 
-  if (!PyArg_ParseTuple(args, "OOiidi",
+  if (!PyArg_ParseTuple(args, "OOlldl",
                         &py_collision_matrix,
                         &py_eigenvalues,
                         &i_sigma,
@@ -1891,7 +1891,7 @@ static void pinv_from_eigensolution(double *data,
                                     const double *eigvals,
                                     const long size,
                                     const double cutoff,
-                                    const int pinv_method)
+                                    const long pinv_method)
 {
   long i, ib, j, k, max_l, i_s, j_s;
   double *tmp_data;
@@ -1976,7 +1976,7 @@ static void show_colmat_info(const PyArrayObject *py_collision_matrix,
                              const long i_temp,
                              const long adrs_shift)
 {
-  int i;
+  long i;
 
   printf(" Array_shape:(");
   for (i = 0; i < PyArray_NDIM(py_collision_matrix); i++) {
@@ -1991,16 +1991,16 @@ static void show_colmat_info(const PyArrayObject *py_collision_matrix,
 }
 
 
-static Iarray* convert_to_iarray(const PyArrayObject* npyary)
+static Larray* convert_to_larray(const PyArrayObject* npyary)
 {
-  int i;
+  long i;
-  Iarray *ary;
+  Larray *ary;
 
-  ary = (Iarray*) malloc(sizeof(Iarray));
+  ary = (Larray*) malloc(sizeof(Larray));
   for (i = 0; i < PyArray_NDIM(npyary); i++) {
     ary->dims[i] = PyArray_DIMS(npyary)[i];
   }
-  ary->data = (int*)PyArray_DATA(npyary);
+  ary->data = (long*)PyArray_DATA(npyary);
   return ary;
 }
 

c/fc3.c

@@ -41,7 +41,7 @@ static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
                               const long j_atom,
                               PHPYCONST double (*delta_fc2s)[3][3],
                               PHPYCONST double (*site_sym_cart)[3][3],
-                              const int *rot_map_sym,
+                              const long *rot_map_sym,
                               const long num_atom,
                               const long num_site_sym,
                               const long num_disp);
@@ -53,46 +53,46 @@ static void tensor3_rotation(double *rot_tensor,
                              const double *rot_cartesian);
 static double tensor3_rotation_elem(const double *tensor,
                                     const double *r,
-                                    const int pos);
+                                    const long pos);
 static void copy_permutation_symmetry_fc3_elem(double *fc3,
                                                const double fc3_elem[27],
-                                               const int a,
+                                               const long a,
-                                               const int b,
+                                               const long b,
-                                               const int c,
+                                               const long c,
                                                const long num_atom);
 static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
                                               const double *fc3,
-                                              const int a,
+                                              const long a,
-                                              const int b,
+                                              const long b,
-                                              const int c,
+                                              const long c,
                                               const long num_atom);
 static void set_permutation_symmetry_compact_fc3(double * fc3,
-                                                 const int p2s[],
+                                                 const long p2s[],
-                                                 const int s2pp[],
+                                                 const long s2pp[],
-                                                 const int nsym_list[],
+                                                 const long nsym_list[],
-                                                 const int perms[],
+                                                 const long perms[],
                                                  const long n_satom,
                                                  const long n_patom);
 static void transpose_compact_fc3_type01(double * fc3,
-                                         const int p2s[],
+                                         const long p2s[],
-                                         const int s2pp[],
+                                         const long s2pp[],
-                                         const int nsym_list[],
+                                         const long nsym_list[],
-                                         const int perms[],
+                                         const long perms[],
                                          const long n_satom,
                                          const long n_patom,
-                                         const int t_type);
+                                         const long t_type);
 static void transpose_compact_fc3_type2(double * fc3,
-                                        const int p2s[],
+                                        const long p2s[],
-                                        const int s2pp[],
+                                        const long s2pp[],
-                                        const int nsym_list[],
+                                        const long nsym_list[],
-                                        const int perms[],
+                                        const long perms[],
                                         const long n_satom,
                                         const long n_patom);
 
 void fc3_distribute_fc3(double *fc3,
-                        const int target,
+                        const long target,
-                        const int source,
+                        const long source,
-                        const int *atom_mapping,
+                        const long *atom_mapping,
                         const long num_atom,
                         const double *rot_cart)
 {
@@ -113,7 +113,7 @@ void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
                           PHPYCONST double (*delta_fc2s)[3][3],
                           const double *inv_U,
                           PHPYCONST double (*site_sym_cart)[3][3],
-                          const int *rot_map_syms,
+                          const long *rot_map_syms,
                           const long num_atom,
                           const long num_site_sym,
                           const long num_disp)
@@ -171,10 +171,10 @@ void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
 }
 
 void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
-                                              const int p2s[],
+                                              const long p2s[],
-                                              const int s2pp[],
+                                              const long s2pp[],
-                                              const int nsym_list[],
+                                              const long nsym_list[],
-                                              const int perms[],
+                                              const long perms[],
                                               const long n_satom,
                                               const long n_patom)
 {
@@ -188,13 +188,13 @@ void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
 }
 
 void fc3_transpose_compact_fc3(double * fc3,
-                               const int p2s[],
+                               const long p2s[],
-                               const int s2pp[],
+                               const long s2pp[],
-                               const int nsym_list[],
+                               const long nsym_list[],
-                               const int perms[],
+                               const long perms[],
                                const long n_satom,
                                const long n_patom,
-                               const int t_type)
+                               const long t_type)
 {
   /* Three types of index permutations                       */
   /*     t_type=0: dim[0] <-> dim[1]                         */
@@ -227,7 +227,7 @@ static void rotate_delta_fc2s(double (*rot_delta_fc2s)[3][3],
                               const long j_atom,
                               PHPYCONST double (*delta_fc2s)[3][3],
                               PHPYCONST double (*site_sym_cart)[3][3],
-                              const int *rot_map_sym,
+                              const long *rot_map_sym,
                               const long num_atom,
                               const long num_site_sym,
                               const long num_disp)
@@ -281,7 +281,7 @@ static void tensor3_rotation(double *rot_tensor,
 
 static double tensor3_rotation_elem(const double *tensor,
                                     const double *r,
-                                    const int pos)
+                                    const long pos)
 {
   long i, j, k, l, m, n;
   double sum;
@@ -304,9 +304,9 @@ static double tensor3_rotation_elem(const double *tensor,
 
 static void copy_permutation_symmetry_fc3_elem(double *fc3,
                                                const double fc3_elem[27],
-                                               const int a,
+                                               const long a,
-                                               const int b,
+                                               const long b,
-                                               const int c,
+                                               const long c,
                                                const long num_atom)
 {
   long i, j, k;
@@ -345,9 +345,9 @@ static void copy_permutation_symmetry_fc3_elem(double *fc3,
 
 static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
                                               const double *fc3,
-                                              const int a,
+                                              const long a,
-                                              const int b,
+                                              const long b,
-                                              const int c,
+                                              const long c,
                                               const long num_atom)
 {
   long i, j, k;
@@ -380,10 +380,10 @@ static void set_permutation_symmetry_fc3_elem(double *fc3_elem,
 }
 
 static void set_permutation_symmetry_compact_fc3(double * fc3,
-                                                 const int p2s[],
+                                                 const long p2s[],
-                                                 const int s2pp[],
+                                                 const long s2pp[],
-                                                 const int nsym_list[],
+                                                 const long nsym_list[],
-                                                 const int perms[],
+                                                 const long perms[],
                                                  const long n_satom,
                                                  const long n_patom)
 {
@@ -392,7 +392,7 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
   /*     i * n_satom * n_satom * 27 + j * n_satom * 27 +     */
   /*     k * 27 + l * 9 + m * 3 + n                          */
   long i, j, k, l, m, n, i_p, j_p, k_p;
-  int done_any;
+  long done_any;
   long i_trans_j, k_trans_j, i_trans_k, j_trans_k;
   long adrs[6];
   double fc3_elem[3][3][3];
@@ -475,13 +475,13 @@ static void set_permutation_symmetry_compact_fc3(double * fc3,
 }
 
 static void transpose_compact_fc3_type01(double * fc3,
-                                         const int p2s[],
+                                         const long p2s[],
-                                         const int s2pp[],
+                                         const long s2pp[],
-                                         const int nsym_list[],
+                                         const long nsym_list[],
-                                         const int perms[],
+                                         const long perms[],
                                          const long n_satom,
                                          const long n_patom,
-                                         const int t_type)
+                                         const long t_type)
 {
   /* Three types of index permutations                       */
   /*     t_type=0: dim[0] <-> dim[1]                         */
@@ -581,10 +581,10 @@ static void transpose_compact_fc3_type01(double * fc3,
 }
 
 static void transpose_compact_fc3_type2(double * fc3,
-                                        const int p2s[],
+                                        const long p2s[],
-                                        const int s2pp[],
+                                        const long s2pp[],
-                                        const int nsym_list[],
+                                        const long nsym_list[],
-                                        const int perms[],
+                                        const long perms[],
                                         const long n_satom,
                                         const long n_patom)
 {

c/fc3.h

@@ -36,34 +36,34 @@
 #define __fc3_H__
 
 void fc3_distribute_fc3(double *fc3,
-                        const int target,
+                        const long target,
-                        const int source,
+                        const long source,
-                        const int *atom_mapping,
+                        const long *atom_mapping,
                         const long num_atom,
                         const double *rot_cart);
 void fc3_rotate_delta_fc2(double (*fc3)[3][3][3],
                           PHPYCONST double (*delta_fc2s)[3][3],
                           const double *inv_U,
                           PHPYCONST double (*site_sym_cart)[3][3],
-                          const int *rot_map_syms,
+                          const long *rot_map_syms,
                           const long num_atom,
                           const long num_site_sym,
                           const long num_disp);
 void fc3_set_permutation_symmetry_fc3(double *fc3, const long num_atom);
 void fc3_set_permutation_symmetry_compact_fc3(double * fc3,
-                                              const int p2s[],
+                                              const long p2s[],
-                                              const int s2pp[],
+                                              const long s2pp[],
-                                              const int nsym_list[],
+                                              const long nsym_list[],
-                                              const int perms[],
+                                              const long perms[],
                                               const long n_satom,
                                               const long n_patom);
 void fc3_transpose_compact_fc3(double * fc3,
-                               const int p2s[],
+                               const long p2s[],
-                               const int s2pp[],
+                               const long s2pp[],
-                               const int nsym_list[],
+                               const long nsym_list[],
-                               const int perms[],
+                               const long perms[],
                                const long n_satom,
                                const long n_patom,
-                               const int t_type);
+                               const long t_type);
 
 #endif

c/imag_self_energy_with_g.c

@@ -41,10 +41,10 @@
 #include "imag_self_energy_with_g.h"
 #include "triplet.h"
 
-static int ise_set_g_pos_frequency_point(int (*g_pos)[4],
+static long ise_set_g_pos_frequency_point(long (*g_pos)[4],
-                                         const long num_band0,
+                                          const long num_band0,
-                                         const long num_band,
+                                          const long num_band,
-                                         const char *g_zero);
+                                          const char *g_zero);
 static void
 detailed_imag_self_energy_at_triplet(double *detailed_imag_self_energy,
                                      double *imag_self_energy,
@@ -93,14 +93,14 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
                                               const char *g_zero,
                                               const double temperature,
                                               const double cutoff_frequency,
-                                              const int num_frequency_points,
+                                              const long num_frequency_points,
-                                              const int frequency_point_index)
+                                              const long frequency_point_index)
 {
   long i, j, num_triplets, num_band0, num_band, num_band_prod;
   long num_g_pos, g_index_dims, g_index_shift;
-  int (*g_pos)[4];
+  long (*g_pos)[4];
   double *ise;
-  int at_a_frequency_point;
+  long at_a_frequency_point;
 
   g_pos = NULL;
   ise = NULL;
@@ -125,7 +125,7 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
 
 #pragma omp parallel for private(num_g_pos, j, g_pos)
   for (i = 0; i < num_triplets; i++) {
-    g_pos = (int(*)[4])malloc(sizeof(int[4]) * num_band_prod);
+    g_pos = (long(*)[4])malloc(sizeof(long[4]) * num_band_prod);
     /* ise_set_g_pos only works for the case of frquency points at */
     /* bands. For frequency sampling mode, g_zero is assumed all */
     /* with the array shape of (num_triplets, num_band0, num_band, */
@@ -187,7 +187,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
  const double *frequencies,
  const long (*triplets)[3],
  const long *triplet_weights,
- const int *grid_address,
+ const long *grid_address,
  const double *g,
  const char *g_zero,
  const double temperature,
@@ -195,7 +195,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
 {
   double *ise;
   long i, j, num_triplets, num_band0, num_band, num_band_prod;
-  int *is_N;
+  long *is_N;
   double ise_tmp, N, U;
 
   ise = NULL;
@@ -227,7 +227,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
        cutoff_frequency);
   }
 
-  is_N = (int*)malloc(sizeof(int) * num_triplets);
+  is_N = (long*)malloc(sizeof(long) * num_triplets);
   for (i = 0; i < num_triplets; i++) {
     is_N[i] = tpl_is_N(triplets[i], grid_address);
   }
@@ -263,17 +263,17 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
                                      const long triplet_weight,
                                      const double *g1,
                                      const double *g2_3,
-                                     PHPYCONST int (*g_pos)[4],
+                                     PHPYCONST long (*g_pos)[4],
                                      const long num_g_pos,
                                      const double *temperatures,
                                      const long num_temps,
                                      const double cutoff_frequency,
-                                     const int openmp_at_bands,
+                                     const long openmp_at_bands,
-                                     const int at_a_frequency_point)
+                                     const long at_a_frequency_point)
 {
   long i, j;
   double *n1, *n2;
-  int g_pos_3;
+  long g_pos_3;
 
   n1 = (double*)malloc(sizeof(double) * num_temps * num_band);
   n2 = (double*)malloc(sizeof(double) * num_temps * num_band);
@@ -328,10 +328,10 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
   n2 = NULL;
 }
 
-int ise_set_g_pos(int (*g_pos)[4],
+long ise_set_g_pos(long (*g_pos)[4],
-                  const long num_band0,
+                   const long num_band0,
-                  const long num_band,
+                   const long num_band,
-                  const char *g_zero)
+                   const char *g_zero)
 {
   long num_g_pos, j, k, l, jkl;
 
@@ -354,10 +354,10 @@ int ise_set_g_pos(int (*g_pos)[4],
   return num_g_pos;
 }
 
-static int ise_set_g_pos_frequency_point(int (*g_pos)[4],
+static long ise_set_g_pos_frequency_point(long (*g_pos)[4],
-                                         const long num_band0,
+                                          const long num_band0,
-                                         const long num_band,
+                                          const long num_band,
-                                         const char *g_zero)
+                                          const char *g_zero)
 {
   long num_g_pos, j, k, l, kl, jkl;
 

c/imag_self_energy_with_g.h

@@ -47,8 +47,8 @@ void ise_get_imag_self_energy_at_bands_with_g(double *imag_self_energy,
                                               const char *g_zero,
                                               const double temperature,
                                               const double cutoff_frequency,
-                                              const int num_frequency_points,
+                                              const long num_frequency_points,
-                                              const int frequency_point_index);
+                                              const long frequency_point_index);
 void ise_get_detailed_imag_self_energy_at_bands_with_g
 (double *detailed_imag_self_energy,
  double *imag_self_energy_N,
@@ -57,7 +57,7 @@ void ise_get_detailed_imag_self_energy_at_bands_with_g
  const double *frequencies,
  const long (*triplets)[3],
  const long *triplet_weights,
- const int *grid_address,
+ const long *grid_address,
  const double *g,
  const char *g_zero,
  const double temperature,
@@ -71,16 +71,16 @@ void ise_imag_self_energy_at_triplet(double *imag_self_energy,
                                      const long triplet_weight,
                                      const double *g1,
                                      const double *g2_3,
-                                     PHPYCONST int (*g_pos)[4],
+                                     PHPYCONST long (*g_pos)[4],
                                      const long num_g_pos,
                                      const double *temperatures,
                                      const long num_temps,
                                      const double cutoff_frequency,
-                                     const int openmp_at_bands,
+                                     const long openmp_at_bands,
-                                     const int at_a_frequency_point);
+                                     const long at_a_frequency_point);
-int ise_set_g_pos(int (*g_pos)[4],
+long ise_set_g_pos(long (*g_pos)[4],
-                  const long num_band0,
+                   const long num_band0,
-                  const long num_band,
+                   const long num_band,
-                  const char *g_zero);
+                   const char *g_zero);
 
 #endif

c/interaction.c

@@ -42,14 +42,14 @@
 #include "reciprocal_to_normal.h"
 #include "lapack_wrapper.h"
 
-static const int index_exchange[6][3] = {{0, 1, 2},
+static const long index_exchange[6][3] = {{0, 1, 2},
-                                         {2, 0, 1},
+                                          {2, 0, 1},
-                                         {1, 2, 0},
+                                          {1, 2, 0},
-                                         {2, 1, 0},
+                                          {2, 1, 0},
-                                         {0, 2, 1},
+                                          {0, 2, 1},
-                                         {1, 0, 2}};
+                                          {1, 0, 2}};
 static void real_to_normal(double *fc3_normal_squared,
-                           PHPYCONST int (*g_pos)[4],
+                           PHPYCONST long (*g_pos)[4],
                            const long num_g_pos,
                            const double *freqs0,
                            const double *freqs1,
@@ -58,42 +58,42 @@ static void real_to_normal(double *fc3_normal_squared,
                            const lapack_complex_double *eigvecs1,
                            const lapack_complex_double *eigvecs2,
                            const double *fc3,
-                           const int is_compact_fc3,
+                           const long is_compact_fc3,
                            const double q[9], /* q0, q1, q2 */
                            const double *shortest_vectors,
-                           const int svecs_dims[3],
+                           const long svecs_dims[3],
-                           const int *multiplicity,
+                           const long *multiplicity,
                            const double *masses,
-                           const int *p2s_map,
+                           const long *p2s_map,
-                           const int *s2p_map,
+                           const long *s2p_map,
-                           const int *band_indices,
+                           const long *band_indices,
                            const long num_band0,
                            const long num_band,
                            const double cutoff_frequency,
                            const long triplet_index,
                            const long num_triplets,
-                           const int openmp_at_bands);
+                           const long openmp_at_bands);
 static void real_to_normal_sym_q(double *fc3_normal_squared,
-                                 PHPYCONST int (*g_pos)[4],
+                                 PHPYCONST long (*g_pos)[4],
                                  const long num_g_pos,
                                  PHPYCONST double *freqs[3],
                                  PHPYCONST lapack_complex_double *eigvecs[3],
                                  const double *fc3,
-                                 const int is_compact_fc3,
+                                 const long is_compact_fc3,
                                  const double q[9], /* q0, q1, q2 */
                                  const double *shortest_vectors,
-                                 const int svecs_dims[3],
+                                 const long svecs_dims[3],
-                                 const int *multiplicity,
+                                 const long *multiplicity,
                                  const double *masses,
-                                 const int *p2s_map,
+                                 const long *p2s_map,
-                                 const int *s2p_map,
+                                 const long *s2p_map,
-                                 const int *band_indices,
+                                 const long *band_indices,
                                  const long num_band0,
                                  const long num_band,
                                  const double cutoff_frequency,
                                  const long triplet_index,
                                  const long num_triplets,
-                                 const int openmp_at_bands);
+                                 const long openmp_at_bands);
 
 /* fc3_normal_squared[num_triplets, num_band0, num_band, num_band] */
 void itr_get_interaction(Darray *fc3_normal_squared,
@@ -102,22 +102,22 @@ void itr_get_interaction(Darray *fc3_normal_squared,
                          const lapack_complex_double *eigenvectors,
                          const long (*triplets)[3],
                          const long num_triplets,
-                         const int *grid_address,
+                         const long *grid_address,
-                         const int *mesh,
+                         const long *mesh,
                          const double *fc3,
-                         const int is_compact_fc3,
+                         const long is_compact_fc3,
                          const double *shortest_vectors,
-                         const int svecs_dims[3],
+                         const long svecs_dims[3],
-                         const int *multiplicity,
+                         const long *multiplicity,
                          const double *masses,
-                         const int *p2s_map,
+                         const long *p2s_map,
-                         const int *s2p_map,
+                         const long *s2p_map,
-                         const int *band_indices,
+                         const long *band_indices,
-                         const int symmetrize_fc3_q,
+                         const long symmetrize_fc3_q,
                          const double cutoff_frequency)
 {
-  int openmp_per_triplets;
+  long openmp_per_triplets;
-  int (*g_pos)[4];
+  long (*g_pos)[4];
   long i;
   long num_band, num_band0, num_band_prod, num_g_pos;
 
@@ -135,7 +135,7 @@ void itr_get_interaction(Darray *fc3_normal_squared,
 
 #pragma omp parallel for schedule(guided) private(num_g_pos, g_pos) if (openmp_per_triplets)
   for (i = 0; i < num_triplets; i++) {
-    g_pos = (int(*)[4])malloc(sizeof(int[4]) * num_band_prod);
+    g_pos = (long(*)[4])malloc(sizeof(long[4]) * num_band_prod);
     num_g_pos = ise_set_g_pos(g_pos,
                               num_band0,
                               num_band,
@@ -175,27 +175,27 @@ void itr_get_interaction(Darray *fc3_normal_squared,
 void itr_get_interaction_at_triplet(double *fc3_normal_squared,
                                     const long num_band0,
                                     const long num_band,
-                                    PHPYCONST int (*g_pos)[4],
+                                    PHPYCONST long (*g_pos)[4],
                                     const long num_g_pos,
                                     const double *frequencies,
                                     const lapack_complex_double *eigenvectors,
                                     const long triplet[3],
-                                    const int *grid_address,
+                                    const long *grid_address,
-                                    const int *mesh,
+                                    const long *mesh,
                                     const double *fc3,
-                                    const int is_compact_fc3,
+                                    const long is_compact_fc3,
                                     const double *shortest_vectors,
-                                    const int svecs_dims[3],
+                                    const long svecs_dims[3],
-                                    const int *multiplicity,
+                                    const long *multiplicity,
                                     const double *masses,
-                                    const int *p2s_map,
+                                    const long *p2s_map,
-                                    const int *s2p_map,
+                                    const long *s2p_map,
-                                    const int *band_indices,
+                                    const long *band_indices,
-                                    const int symmetrize_fc3_q,
+                                    const long symmetrize_fc3_q,
                                     const double cutoff_frequency,
                                     const long triplet_index, /* only for print */
                                     const long num_triplets, /* only for print */
-                                    const int openmp_at_bands)
+                                    const long openmp_at_bands)
 {
   long j, k;
   double *freqs[3];
@@ -277,7 +277,7 @@ void itr_get_interaction_at_triplet(double *fc3_normal_squared,
 }
 
 static void real_to_normal(double *fc3_normal_squared,
-                           PHPYCONST int (*g_pos)[4],
+                           PHPYCONST long (*g_pos)[4],
                            const long num_g_pos,
                            const double *freqs0,
                            const double *freqs1,
@@ -286,21 +286,21 @@ static void real_to_normal(double *fc3_normal_squared,
                            const lapack_complex_double *eigvecs1,
                            const lapack_complex_double *eigvecs2,
                            const double *fc3,
-                           const int is_compact_fc3,
+                           const long is_compact_fc3,
                            const double q[9], /* q0, q1, q2 */
                            const double *shortest_vectors,
-                           const int svecs_dims[3],
+                           const long svecs_dims[3],
-                           const int *multiplicity,
+                           const long *multiplicity,
                            const double *masses,
-                           const int *p2s_map,
+                           const long *p2s_map,
-                           const int *s2p_map,
+                           const long *s2p_map,
-                           const int *band_indices,
+                           const long *band_indices,
                            const long num_band0,
                            const long num_band,
                            const double cutoff_frequency,
                            const long triplet_index,
                            const long num_triplets,
-                           const int openmp_at_bands)
+                           const long openmp_at_bands)
 {
   long num_patom;
   lapack_complex_double *fc3_reciprocal;
@@ -349,26 +349,26 @@ static void real_to_normal(double *fc3_normal_squared,
 }
 
 static void real_to_normal_sym_q(double *fc3_normal_squared,
-                                 PHPYCONST int (*g_pos)[4],
+                                 PHPYCONST long (*g_pos)[4],
                                  const long num_g_pos,
                                  PHPYCONST double *freqs[3],
                                  PHPYCONST lapack_complex_double *eigvecs[3],
                                  const double *fc3,
-                                 const int is_compact_fc3,
+                                 const long is_compact_fc3,
                                  const double q[9], /* q0, q1, q2 */
                                  const double *shortest_vectors,
-                                 const int svecs_dims[3],
+                                 const long svecs_dims[3],
-                                 const int *multiplicity,
+                                 const long *multiplicity,
                                  const double *masses,
-                                 const int *p2s_map,
+                                 const long *p2s_map,
-                                 const int *s2p_map,
+                                 const long *s2p_map,
-                                 const int *band_indices,
+                                 const long *band_indices,
                                  const long num_band0,
                                  const long num_band,
                                  const double cutoff_frequency,
                                  const long triplet_index,
                                  const long num_triplets,
-                                 const int openmp_at_bands)
+                                 const long openmp_at_bands)
 {
   long i, j, k, l;
   long band_ex[3];

c/interaction.h

@@ -44,42 +44,42 @@ void itr_get_interaction(Darray *fc3_normal_squared,
                          const lapack_complex_double *eigenvectors,
                          const long (*triplets)[3],
                          const long num_triplets,
-                         const int *grid_address,
+                         const long *grid_address,
-                         const int *mesh,
+                         const long *mesh,
                          const double *fc3,
-                         const int is_compact_fc3,
+                         const long is_compact_fc3,
                          const double *shortest_vectors,
-                         const int svecs_dims[3],
+                         const long svecs_dims[3],
-                         const int *multiplicity,
+                         const long *multiplicity,
                          const double *masses,
-                         const int *p2s_map,
+                         const long *p2s_map,
-                         const int *s2p_map,
+                         const long *s2p_map,
-                         const int *band_indices,
+                         const long *band_indices,
-                         const int symmetrize_fc3_q,
+                         const long symmetrize_fc3_q,
                          const double cutoff_frequency);
 void itr_get_interaction_at_triplet(double *fc3_normal_squared,
                                     const long num_band0,
                                     const long num_band,
-                                    PHPYCONST int (*g_pos)[4],
+                                    PHPYCONST long (*g_pos)[4],
                                     const long num_g_pos,
                                     const double *frequencies,
                                     const lapack_complex_double *eigenvectors,
                                     const long triplet[3],
-                                    const int *grid_address,
+                                    const long *grid_address,
-                                    const int *mesh,
+                                    const long *mesh,
                                     const double *fc3,
-                                    const int is_compact_fc3,
+                                    const long is_compact_fc3,
                                     const double *shortest_vectors,
-                                    const int svecs_dims[3],
+                                    const long svecs_dims[3],
-                                    const int *multiplicity,
+                                    const long *multiplicity,
                                     const double *masses,
-                                    const int *p2s_map,
+                                    const long *p2s_map,
-                                    const int *s2p_map,
+                                    const long *s2p_map,
-                                    const int *band_indices,
+                                    const long *band_indices,
-                                    const int symmetrize_fc3_q,
+                                    const long symmetrize_fc3_q,
                                     const double cutoff_frequency,
                                     const long triplet_index, /* only for print */
                                     const long num_triplets, /* only for print */
-                                    const int openmp_at_bands);
+                                    const long openmp_at_bands);
 
 #endif

c/isotope.c

@@ -45,7 +45,7 @@ iso_get_isotope_scattering_strength(double *gamma,
                                     const double *frequencies,
                                     const lapack_complex_double *eigenvectors,
                                     const long num_grid_points,
-                                    const int *band_indices,
+                                    const long *band_indices,
                                     const long num_band,
                                     const long num_band0,
                                     const double sigma,
@@ -135,7 +135,7 @@ void iso_get_thm_isotope_scattering_strength
  const double *frequencies,
  const lapack_complex_double *eigenvectors,
  const long num_grid_points,
- const int *band_indices,
+ const long *band_indices,
  const long num_band,
  const long num_band0,
  const double *integration_weights,

c/isotope.h

@@ -44,7 +44,7 @@ iso_get_isotope_scattering_strength(double *gamma,
                                     const double *frequencies,
                                     const lapack_complex_double *eigenvectors,
                                     const long num_grid_points,
-                                    const int *band_indices,
+                                    const long *band_indices,
                                     const long num_band,
                                     const long num_band0,
                                     const double sigma,
@@ -58,7 +58,7 @@ void iso_get_thm_isotope_scattering_strength
  const double *frequencies,
  const lapack_complex_double *eigenvectors,
  const long num_grid_points,
- const int *band_indices,
+ const long *band_indices,
  const long num_band,
  const long num_band0,
  const double *integration_weights,

c/kpoint.c

@@ -37,7 +37,7 @@
 #include <stddef.h>
 #include "mathfunc.h"
 #include "kpoint.h"
-#include "kgrid.h"
+#include "rgrid.h"
 
 #ifdef KPTWARNING
 #include <stdio.h>
@@ -46,39 +46,51 @@
 #define warning_print(...)
 #endif
 
-static MatINT *get_point_group_reciprocal(const MatINT * rotations,
+static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
-                                          const int is_time_reversal);
+                                           const long is_time_reversal);
-static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
+static MatLONG *get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
-                                                 const double symprec,
+                                                  const double symprec,
-                                                 const long num_q,
+                                                  const long num_q,
-                                                 SPGCONST double qpoints[][3]);
+                                                  KPTCONST double qpoints[][3]);
-static long get_dense_ir_reciprocal_mesh(int grid_address[][3],
+static long get_dense_ir_reciprocal_mesh(long grid_address[][3],
                                          long ir_mapping_table[],
-                                         const int mesh[3],
+                                         const long mesh[3],
-                                         const int is_shift[3],
+                                         const long is_shift[3],
-                                         const MatINT *rot_reciprocal);
+                                         const MatLONG *rot_reciprocal);
-static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
+static long get_dense_ir_reciprocal_mesh_normal(long grid_address[][3],
                                                 long ir_mapping_table[],
-                                                const int mesh[3],
+                                                const long mesh[3],
-                                                const int is_shift[3],
+                                                const long is_shift[3],
-                                                const MatINT *rot_reciprocal);
+                                                const MatLONG *rot_reciprocal);
-static long get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
+static long get_dense_ir_reciprocal_mesh_distortion(long grid_address[][3],
                                                     long ir_mapping_table[],
-                                                    const int mesh[3],
+                                                    const long mesh[3],
-                                                    const int is_shift[3],
+                                                    const long is_shift[3],
-                                                    const MatINT *rot_reciprocal);
+                                                    const MatLONG *rot_reciprocal);
-static long get_dense_num_ir(long ir_mapping_table[], const int mesh[3]);
+static long get_dense_num_ir(long ir_mapping_table[], const long mesh[3]);
-static int check_mesh_symmetry(const int mesh[3],
+static long check_mesh_symmetry(const long mesh[3],
-                               const int is_shift[3],
+                                const long is_shift[3],
-                               const MatINT *rot_reciprocal);
+                                const MatLONG *rot_reciprocal);
+static long Nint(const double a);
+static double Dabs(const double a);
+static void transpose_matrix_l3(long a[3][3], KPTCONST long b[3][3]);
+static void multiply_matrix_l3(long m[3][3],
+                               KPTCONST long a[3][3], KPTCONST long b[3][3]);
+static long check_identity_matrix_l3(KPTCONST long a[3][3],
+                                     KPTCONST long b[3][3]);
+static void multiply_matrix_vector_ld3(double v[3],
+                                       KPTCONST long a[3][3],
+                                       const double b[3]);
+static void multiply_matrix_vector_l3(long v[3],
+                                      KPTCONST long a[3][3],
+                                      const long b[3]);
 
 
-
+long kpt_get_dense_irreducible_reciprocal_mesh(long grid_address[][3],
-long kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
                                                long ir_mapping_table[],
-                                               const int mesh[3],
+                                               const long mesh[3],
-                                               const int is_shift[3],
+                                               const long is_shift[3],
-                                               const MatINT *rot_reciprocal)
+                                               const MatLONG *rot_reciprocal)
 {
   long num_ir;
 
@@ -91,16 +103,16 @@ long kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
   return num_ir;
 }
 
-MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
+MatLONG *kpt_get_point_group_reciprocal(const MatLONG * rotations,
-                                       const int is_time_reversal)
+                                        const long is_time_reversal)
 {
   return get_point_group_reciprocal(rotations, is_time_reversal);
 }
 
-MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
+MatLONG *kpt_get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
-                                              const double symprec,
+                                               const double symprec,
-                                              const long num_q,
+                                               const long num_q,
-                                              SPGCONST double qpoints[][3])
+                                               KPTCONST double qpoints[][3])
 {
   return get_point_group_reciprocal_with_q(rot_reciprocal,
                                            symprec,
@@ -108,14 +120,62 @@ MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
                                            qpoints);
 }
 
-/* Return NULL if failed */
+void kpt_copy_matrix_l3(long a[3][3], KPTCONST long b[3][3])
-static MatINT *get_point_group_reciprocal(const MatINT * rotations,
-                                          const int is_time_reversal)
 {
-  int i, j, num_rot;
+  a[0][0] = b[0][0];
-  MatINT *rot_reciprocal, *rot_return;
+  a[0][1] = b[0][1];
-  int *unique_rot;
+  a[0][2] = b[0][2];
-  SPGCONST int inversion[3][3] = {
+  a[1][0] = b[1][0];
+  a[1][1] = b[1][1];
+  a[1][2] = b[1][2];
+  a[2][0] = b[2][0];
+  a[2][1] = b[2][1];
+  a[2][2] = b[2][2];
+}
+
+MatLONG * kpt_alloc_MatLONG(const long size)
+{
+  MatLONG *matlong;
+
+  matlong = NULL;
+
+  if ((matlong = (MatLONG*) malloc(sizeof(MatLONG))) == NULL) {
+    warning_print("spglib: Memory could not be allocated.");
+    return NULL;
+  }
+
+  matlong->size = size;
+  if (size > 0) {
+    if ((matlong->mat = (long (*)[3][3]) malloc(sizeof(long[3][3]) * size))
+        == NULL) {
+      warning_print("spglib: Memory could not be allocated ");
+      warning_print("(MatLONG, line %d, %s).\n", __LINE__, __FILE__);
+      free(matlong);
+      matlong = NULL;
+      return NULL;
+    }
+  }
+  return matlong;
+}
+
+void kpt_free_MatLONG(MatLONG * matlong)
+{
+  if (matlong->size > 0) {
+    free(matlong->mat);
+    matlong->mat = NULL;
+  }
+  free(matlong);
+}
+
+
+/* Return NULL if failed */
+static MatLONG *get_point_group_reciprocal(const MatLONG * rotations,
+                                           const long is_time_reversal)
+{
+  long i, j, num_rot;
+  MatLONG *rot_reciprocal, *rot_return;
+  long *unique_rot;
+  KPTCONST long inversion[3][3] = {
     {-1, 0, 0 },
     { 0,-1, 0 },
     { 0, 0,-1 }
@@ -126,18 +186,18 @@ static MatINT *get_point_group_reciprocal(const MatINT * rotations,
   unique_rot = NULL;
 
   if (is_time_reversal) {
-    if ((rot_reciprocal = mat_alloc_MatINT(rotations->size * 2)) == NULL) {
+    if ((rot_reciprocal = kpt_alloc_MatLONG(rotations->size * 2)) == NULL) {
       return NULL;
     }
   } else {
-    if ((rot_reciprocal = mat_alloc_MatINT(rotations->size)) == NULL) {
+    if ((rot_reciprocal = kpt_alloc_MatLONG(rotations->size)) == NULL) {
       return NULL;
     }
   }
 
-  if ((unique_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) {
+  if ((unique_rot = (long*)malloc(sizeof(long) * rot_reciprocal->size)) == NULL) {
     warning_print("spglib: Memory of unique_rot could not be allocated.");
-    mat_free_MatINT(rot_reciprocal);
+    kpt_free_MatLONG(rot_reciprocal);
     rot_reciprocal = NULL;
     return NULL;
   }
@@ -147,20 +207,20 @@ static MatINT *get_point_group_reciprocal(const MatINT * rotations,
   }
 
   for (i = 0; i < rotations->size; i++) {
-    mat_transpose_matrix_i3(rot_reciprocal->mat[i], rotations->mat[i]);
+    transpose_matrix_l3(rot_reciprocal->mat[i], rotations->mat[i]);
 
     if (is_time_reversal) {
-      mat_multiply_matrix_i3(rot_reciprocal->mat[rotations->size+i],
+      multiply_matrix_l3(rot_reciprocal->mat[rotations->size+i],
-                             inversion,
+                         inversion,
-                             rot_reciprocal->mat[i]);
+                         rot_reciprocal->mat[i]);
     }
   }
 
   num_rot = 0;
   for (i = 0; i < rot_reciprocal->size; i++) {
     for (j = 0; j < num_rot; j++) {
-      if (mat_check_identity_matrix_i3(rot_reciprocal->mat[unique_rot[j]],
+      if (check_identity_matrix_l3(rot_reciprocal->mat[unique_rot[j]],
-                                       rot_reciprocal->mat[i])) {
+                                   rot_reciprocal->mat[i])) {
         goto escape;
       }
     }
@@ -170,37 +230,37 @@ static MatINT *get_point_group_reciprocal(const MatINT * rotations,
     ;
   }
 
-  if ((rot_return = mat_alloc_MatINT(num_rot)) != NULL) {
+  if ((rot_return = kpt_alloc_MatLONG(num_rot)) != NULL) {
     for (i = 0; i < num_rot; i++) {
-      mat_copy_matrix_i3(rot_return->mat[i], rot_reciprocal->mat[unique_rot[i]]);
+      kpt_copy_matrix_l3(rot_return->mat[i], rot_reciprocal->mat[unique_rot[i]]);
     }
   }
 
   free(unique_rot);
   unique_rot = NULL;
-  mat_free_MatINT(rot_reciprocal);
+  kpt_free_MatLONG(rot_reciprocal);
   rot_reciprocal = NULL;
 
   return rot_return;
 }
 
 /* Return NULL if failed */
-static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
+static MatLONG *get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
-                                                 const double symprec,
+                                                  const double symprec,
-                                                 const long num_q,
+                                                  const long num_q,
-                                                 SPGCONST double qpoints[][3])
+                                                  KPTCONST double qpoints[][3])
 {
-  int i, j, k, l, is_all_ok, num_rot;
+  long i, j, k, l, is_all_ok, num_rot;
-  int *ir_rot;
+  long *ir_rot;
   double q_rot[3], diff[3];
-  MatINT * rot_reciprocal_q;
+  MatLONG * rot_reciprocal_q;
 
   ir_rot = NULL;
   rot_reciprocal_q = NULL;
   is_all_ok = 0;
   num_rot = 0;
 
-  if ((ir_rot = (int*)malloc(sizeof(int) * rot_reciprocal->size)) == NULL) {
+  if ((ir_rot = (long*)malloc(sizeof(long) * rot_reciprocal->size)) == NULL) {
     warning_print("spglib: Memory of ir_rot could not be allocated.");
     return NULL;
   }
@@ -211,19 +271,19 @@ static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
   for (i = 0; i < rot_reciprocal->size; i++) {
     for (j = 0; j < num_q; j++) {
       is_all_ok = 0;
-      mat_multiply_matrix_vector_id3(q_rot,
+      multiply_matrix_vector_ld3(q_rot,
-                                     rot_reciprocal->mat[i],
+                                 rot_reciprocal->mat[i],
-                                     qpoints[j]);
+                                 qpoints[j]);
 
       for (k = 0; k < num_q; k++) {
         for (l = 0; l < 3; l++) {
           diff[l] = q_rot[l] - qpoints[k][l];
-          diff[l] -= mat_Nint(diff[l]);
+          diff[l] -= Nint(diff[l]);
         }
 
-        if (mat_Dabs(diff[0]) < symprec &&
+        if (Dabs(diff[0]) < symprec &&
-            mat_Dabs(diff[1]) < symprec &&
+            Dabs(diff[1]) < symprec &&
-            mat_Dabs(diff[2]) < symprec) {
+            Dabs(diff[2]) < symprec) {
           is_all_ok = 1;
           break;
         }
@@ -240,9 +300,9 @@ static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
     }
   }
 
-  if ((rot_reciprocal_q = mat_alloc_MatINT(num_rot)) != NULL) {
+  if ((rot_reciprocal_q = kpt_alloc_MatLONG(num_rot)) != NULL) {
     for (i = 0; i < num_rot; i++) {
-      mat_copy_matrix_i3(rot_reciprocal_q->mat[i],
+      kpt_copy_matrix_l3(rot_reciprocal_q->mat[i],
                          rot_reciprocal->mat[ir_rot[i]]);
     }
   }
@@ -253,11 +313,11 @@ static MatINT *get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
   return rot_reciprocal_q;
 }
 
-static long get_dense_ir_reciprocal_mesh(int grid_address[][3],
+static long get_dense_ir_reciprocal_mesh(long grid_address[][3],
                                          long ir_mapping_table[],
-                                         const int mesh[3],
+                                         const long mesh[3],
-                                         const int is_shift[3],
+                                         const long is_shift[3],
-                                         const MatINT *rot_reciprocal)
+                                         const MatLONG *rot_reciprocal)
 {
   if (check_mesh_symmetry(mesh, is_shift, rot_reciprocal)) {
     return get_dense_ir_reciprocal_mesh_normal(grid_address,
@@ -274,11 +334,11 @@ static long get_dense_ir_reciprocal_mesh(int grid_address[][3],
   }
 }
 
-static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
+static long get_dense_ir_reciprocal_mesh_normal(long grid_address[][3],
                                                 long ir_mapping_table[],
-                                                const int mesh[3],
+                                                const long mesh[3],
-                                                const int is_shift[3],
+                                                const long is_shift[3],
-                                                const MatINT *rot_reciprocal)
+                                                const MatLONG *rot_reciprocal)
 {
   /* In the following loop, mesh is doubled. */
   /* Even and odd mesh numbers correspond to */
@@ -288,23 +348,23 @@ static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
   /* ir_mapping_table: the mapping from each point to ir-point. */
 
   long i, grid_point_rot;
-  int j;
+  long j;
-  int address_double[3], address_double_rot[3];
+  long address_double[3], address_double_rot[3];
 
-  kgd_get_all_grid_addresses(grid_address, mesh);
+  rgd_get_all_grid_addresses(grid_address, mesh);
 
 #pragma omp parallel for private(j, grid_point_rot, address_double, address_double_rot)
   for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
-    kgd_get_grid_address_double_mesh(address_double,
+    rgd_get_double_grid_address(address_double,
-                                     grid_address[i],
+                                grid_address[i],
-                                     mesh,
+                                mesh,
-                                     is_shift);
+                                is_shift);
     ir_mapping_table[i] = i;
     for (j = 0; j < rot_reciprocal->size; j++) {
-      mat_multiply_matrix_vector_i3(address_double_rot,
+      multiply_matrix_vector_l3(address_double_rot,
-                                    rot_reciprocal->mat[j],
+                                rot_reciprocal->mat[j],
-                                    address_double);
+                                address_double);
-      grid_point_rot = kgd_get_dense_grid_point_double_mesh(address_double_rot, mesh);
+      grid_point_rot = rgd_get_double_grid_index(address_double_rot, mesh);
       if (grid_point_rot < ir_mapping_table[i]) {
 #ifdef _OPENMP
         ir_mapping_table[i] = grid_point_rot;
@@ -320,21 +380,21 @@ static long get_dense_ir_reciprocal_mesh_normal(int grid_address[][3],
 }
 
 static long
-get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
+get_dense_ir_reciprocal_mesh_distortion(long grid_address[][3],
                                         long ir_mapping_table[],
-                                        const int mesh[3],
+                                        const long mesh[3],
-                                        const int is_shift[3],
+                                        const long is_shift[3],
-                                        const MatINT *rot_reciprocal)
+                                        const MatLONG *rot_reciprocal)
 {
   long i, grid_point_rot;
-  int j, k, indivisible;
+  long j, k, indivisible;
-  int address_double[3], address_double_rot[3];
+  long address_double[3], address_double_rot[3];
   long long_address_double[3], long_address_double_rot[3], divisor[3];
 
   /* divisor, long_address_double, and long_address_double_rot have */
   /* long integer type to treat dense mesh. */
 
-  kgd_get_all_grid_addresses(grid_address, mesh);
+  rgd_get_all_grid_addresses(grid_address, mesh);
 
   for (j = 0; j < 3; j++) {
     divisor[j] = mesh[(j + 1) % 3] * mesh[(j + 2) % 3];
@@ -342,10 +402,10 @@ get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
 
 #pragma omp parallel for private(j, k, grid_point_rot, address_double, address_double_rot, long_address_double, long_address_double_rot)
   for (i = 0; i < mesh[0] * mesh[1] * (long)(mesh[2]); i++) {
-    kgd_get_grid_address_double_mesh(address_double,
+    rgd_get_double_grid_address(address_double,
-                                     grid_address[i],
+                                grid_address[i],
-                                     mesh,
+                                mesh,
-                                     is_shift);
+                                is_shift);
     for (j = 0; j < 3; j++) {
       long_address_double[j] = address_double[j] * divisor[j];
     }
@@ -372,7 +432,7 @@ get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
       }
       if (indivisible) {continue;}
       grid_point_rot =
-        kgd_get_dense_grid_point_double_mesh(address_double_rot, mesh);
+        rgd_get_double_grid_index(address_double_rot, mesh);
       if (grid_point_rot < ir_mapping_table[i]) {
 #ifdef _OPENMP
         ir_mapping_table[i] = grid_point_rot;
@@ -387,7 +447,7 @@ get_dense_ir_reciprocal_mesh_distortion(int grid_address[][3],
   return get_dense_num_ir(ir_mapping_table, mesh);
 }
 
-static long get_dense_num_ir(long ir_mapping_table[], const int mesh[3])
+static long get_dense_num_ir(long ir_mapping_table[], const long mesh[3])
 {
   long i, num_ir;
 
@@ -409,12 +469,12 @@ static long get_dense_num_ir(long ir_mapping_table[], const int mesh[3])
   return num_ir;
 }
 
-static int check_mesh_symmetry(const int mesh[3],
+static long check_mesh_symmetry(const long mesh[3],
-                               const int is_shift[3],
+                                const long is_shift[3],
-                               const MatINT *rot_reciprocal)
+                                const MatLONG *rot_reciprocal)
 {
-  int i, j, k, sum;
+  long i, j, k, sum;
-  int eq[3];
+  long eq[3];
 
   eq[0] = 0; /* a=b */
   eq[1] = 0; /* b=c */
@@ -425,7 +485,7 @@ static int check_mesh_symmetry(const int mesh[3],
     sum = 0;
     for (j = 0; j < 3; j++) {
       for (k = 0; k < 3; k++) {
-        sum += abs(rot_reciprocal->mat[i][j][k]);
+        sum += labs(rot_reciprocal->mat[i][j][k]);
       }
     }
     if (sum > 3) {
@@ -450,3 +510,93 @@ static int check_mesh_symmetry(const int mesh[3],
           ((eq[1] && mesh[1] == mesh[2] && is_shift[1] == is_shift[2]) || (!eq[1])) &&
           ((eq[2] && mesh[2] == mesh[0] && is_shift[2] == is_shift[0]) || (!eq[2])));
 }
+
+
+static long Nint(const double a)
+{
+  if (a < 0.0)
+    return (long) (a - 0.5);
+  else
+    return (long) (a + 0.5);
+}
+
+static double Dabs(const double a)
+{
+  if (a < 0.0)
+    return -a;
+  else
+    return a;
+}
+
+static void transpose_matrix_l3(long a[3][3], KPTCONST long b[3][3])
+{
+  long c[3][3];
+  c[0][0] = b[0][0];
+  c[0][1] = b[1][0];
+  c[0][2] = b[2][0];
+  c[1][0] = b[0][1];
+  c[1][1] = b[1][1];
+  c[1][2] = b[2][1];
+  c[2][0] = b[0][2];
+  c[2][1] = b[1][2];
+  c[2][2] = b[2][2];
+  kpt_copy_matrix_l3(a, c);
+}
+
+static void multiply_matrix_l3(long m[3][3],
+                               KPTCONST long a[3][3],
+                               KPTCONST long b[3][3])
+{
+  long i, j;                   /* a_ij */
+  long c[3][3];
+  for (i = 0; i < 3; i++) {
+    for (j = 0; j < 3; j++) {
+      c[i][j] =
+        a[i][0] * b[0][j] + a[i][1] * b[1][j] + a[i][2] * b[2][j];
+    }
+  }
+  kpt_copy_matrix_l3(m, c);
+}
+
+static long check_identity_matrix_l3(KPTCONST long a[3][3],
+                                     KPTCONST long b[3][3])
+{
+  if ( a[0][0] - b[0][0] ||
+       a[0][1] - b[0][1] ||
+       a[0][2] - b[0][2] ||
+       a[1][0] - b[1][0] ||
+       a[1][1] - b[1][1] ||
+       a[1][2] - b[1][2] ||
+       a[2][0] - b[2][0] ||
+       a[2][1] - b[2][1] ||
+       a[2][2] - b[2][2]) {
+    return 0;
+  }
+  else {
+    return 1;
+  }
+}
+
+static void multiply_matrix_vector_ld3(double v[3],
+                                       KPTCONST long a[3][3],
+                                       const double b[3])
+{
+  long i;
+  double c[3];
+  for (i = 0; i < 3; i++)
+    c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
+  for (i = 0; i < 3; i++)
+    v[i] = c[i];
+}
+
+static void multiply_matrix_vector_l3(long v[3],
+                                      KPTCONST long a[3][3],
+                                      const long b[3])
+{
+  long i;
+  long c[3];
+  for (i = 0; i < 3; i++)
+    c[i] = a[i][0] * b[0] + a[i][1] * b[1] + a[i][2] * b[2];
+  for (i = 0; i < 3; i++)
+    v[i] = c[i];
+}

c/kpoint.h

@@ -35,19 +35,28 @@
 #ifndef __kpoint_H__
 #define __kpoint_H__
 
-#include <stddef.h>
+#ifndef KPTCONST
-#include "mathfunc.h"
+#define KPTCONST
+#endif
 
-long kpt_get_dense_irreducible_reciprocal_mesh(int grid_address[][3],
+typedef struct {
+  long size;
+  long (*mat)[3][3];
+} MatLONG;
+
+long kpt_get_dense_irreducible_reciprocal_mesh(long grid_address[][3],
                                                long ir_mapping_table[],
-                                               const int mesh[3],
+                                               const long mesh[3],
-                                               const int is_shift[3],
+                                               const long is_shift[3],
-                                               const MatINT *rot_reciprocal);
+                                               const MatLONG *rot_reciprocal);
-MatINT *kpt_get_point_group_reciprocal(const MatINT * rotations,
+MatLONG *kpt_get_point_group_reciprocal(const MatLONG * rotations,
-                                       const int is_time_reversal);
+                                        const long is_time_reversal);
-MatINT *kpt_get_point_group_reciprocal_with_q(const MatINT * rot_reciprocal,
+MatLONG *kpt_get_point_group_reciprocal_with_q(const MatLONG * rot_reciprocal,
-                                              const double symprec,
+                                               const double symprec,
-                                              const long num_q,
+                                               const long num_q,
-                                              SPGCONST double qpoints[][3]);
+                                               KPTCONST double qpoints[][3]);
+void kpt_copy_matrix_l3(long a[3][3], KPTCONST long b[3][3]);
+MatLONG * kpt_alloc_MatLONG(const long size);
+void kpt_free_MatLONG(MatLONG * matlong);
 
 #endif

c/phono3py.c

@@ -47,8 +47,6 @@
 #include "triplet.h"
 #include "triplet_iw.h"
 
-#include "kgrid.h"
-
 #include <stdio.h>
 
 
@@ -58,18 +56,18 @@ void ph3py_get_interaction(Darray *fc3_normal_squared,
                            const lapack_complex_double *eigenvectors,
                            const long (*triplets)[3],
                            const long num_triplets,
-                           const int *grid_address,
+                           const long *grid_address,
-                           const int *mesh,
+                           const long *mesh,
                            const double *fc3,
-                           const int is_compact_fc3,
+                           const long is_compact_fc3,
                            const double *shortest_vectors,
-                           const int svecs_dims[3],
+                           const long svecs_dims[3],
-                           const int *multiplicity,
+                           const long *multiplicity,
                            const double *masses,
-                           const int *p2s_map,
+                           const long *p2s_map,
-                           const int *s2p_map,
+                           const long *s2p_map,
-                           const int *band_indices,
+                           const long *band_indices,
-                           const int symmetrize_fc3_q,
+                           const long symmetrize_fc3_q,
                            const double cutoff_frequency)
 {
   itr_get_interaction(fc3_normal_squared,
@@ -95,27 +93,27 @@ void ph3py_get_interaction(Darray *fc3_normal_squared,
 
 
 void ph3py_get_pp_collision(double *imag_self_energy,
-                            PHPYCONST int relative_grid_address[24][4][3], /* thm */
+                            PHPYCONST long relative_grid_address[24][4][3], /* thm */
                             const double *frequencies,
                             const lapack_complex_double *eigenvectors,
                             const long (*triplets)[3],
                             const long num_triplets,
                             const long *triplet_weights,
-                            const int *grid_address, /* thm */
+                            const long *grid_address, /* thm */
                             const long *bz_map, /* thm */
-                            const int *mesh, /* thm */
+                            const long *mesh, /* thm */
                             const double *fc3,
-                            const int is_compact_fc3,
+                            const long is_compact_fc3,
                             const double *shortest_vectors,
-                            const int svecs_dims[3],
+                            const long svecs_dims[3],
-                            const int *multiplicity,
+                            const long *multiplicity,
                             const double *masses,
-                            const int *p2s_map,
+                            const long *p2s_map,
-                            const int *s2p_map,
+                            const long *s2p_map,
-                            const Iarray *band_indices,
+                            const Larray *band_indices,
                             const Darray *temperatures,
-                            const int is_NU,
+                            const long is_NU,
-                            const int symmetrize_fc3_q,
+                            const long symmetrize_fc3_q,
                             const double cutoff_frequency)
 {
   ppc_get_pp_collision(imag_self_energy,
@@ -153,20 +151,20 @@ void ph3py_get_pp_collision_with_sigma(
   const long (*triplets)[3],
   const long num_triplets,
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
-  const int *mesh,
+  const long *mesh,
   const double *fc3,
-  const int is_compact_fc3,
+  const long is_compact_fc3,
   const double *shortest_vectors,
-  const int svecs_dims[3],
+  const long svecs_dims[3],
-  const int *multiplicity,
+  const long *multiplicity,
   const double *masses,
-  const int *p2s_map,
+  const long *p2s_map,
-  const int *s2p_map,
+  const long *s2p_map,
-  const Iarray *band_indices,
+  const Larray *band_indices,
   const Darray *temperatures,
-  const int is_NU,
+  const long is_NU,
-  const int symmetrize_fc3_q,
+  const long symmetrize_fc3_q,
   const double cutoff_frequency)
 {
   ppc_get_pp_collision_with_sigma(imag_self_energy,
@@ -205,8 +203,8 @@ void ph3py_get_imag_self_energy_at_bands_with_g(
   const char *g_zero,
   const double temperature,
   const double cutoff_frequency,
-  const int num_frequency_points,
+  const long num_frequency_points,
-  const int frequency_point_index)
+  const long frequency_point_index)
 {
   ise_get_imag_self_energy_at_bands_with_g(imag_self_energy,
                                            fc3_normal_squared,
@@ -230,7 +228,7 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
   const double *g,
   const char *g_zero,
   const double temperature,
@@ -253,7 +251,7 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
 
 void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
                                          const Darray *fc3_normal_squared,
-                                         const int *band_indices,
+                                         const long *band_indices,
                                          const double *frequencies,
                                          const long (*triplets)[3],
                                          const long *triplet_weights,
@@ -279,7 +277,7 @@ void ph3py_get_real_self_energy_at_frequency_point(
   double *real_self_energy,
   const double frequency_point,
   const Darray *fc3_normal_squared,
-  const int *band_indices,
+  const long *band_indices,
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,
@@ -369,7 +367,7 @@ void ph3py_get_isotope_scattering_strength(
   const double *frequencies,
   const lapack_complex_double *eigenvectors,
   const long num_grid_points,
-  const int *band_indices,
+  const long *band_indices,
   const long num_band,
   const long num_band0,
   const double sigma,
@@ -398,7 +396,7 @@ void ph3py_get_thm_isotope_scattering_strength
  const double *frequencies,
  const lapack_complex_double *eigenvectors,
  const long num_ir_grid_points,
- const int *band_indices,
+ const long *band_indices,
  const long num_band,
  const long num_band0,
  const double *integration_weights,
@@ -420,9 +418,9 @@ void ph3py_get_thm_isotope_scattering_strength
 }
 
 void ph3py_distribute_fc3(double *fc3,
-                          const int target,
+                          const long target,
-                          const int source,
+                          const long source,
-                          const int *atom_mapping,
+                          const long *atom_mapping,
                           const long num_atom,
                           const double *rot_cart)
 {
@@ -439,7 +437,7 @@ void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
                             PHPYCONST double (*delta_fc2s)[3][3],
                             const double *inv_U,
                             PHPYCONST double (*site_sym_cart)[3][3],
-                            const int *rot_map_syms,
+                            const long *rot_map_syms,
                             const long num_atom,
                             const long num_site_sym,
                             const long num_disp)
@@ -462,10 +460,10 @@ void ph3py_set_permutation_symmetry_fc3(double *fc3, const long num_atom)
 
 
 void ph3py_set_permutation_symmetry_compact_fc3(double * fc3,
-                                                const int p2s[],
+                                                const long p2s[],
-                                                const int s2pp[],
+                                                const long s2pp[],
-                                                const int nsym_list[],
+                                                const long nsym_list[],
-                                                const int perms[],
+                                                const long perms[],
                                                 const long n_satom,
                                                 const long n_patom)
 {
@@ -479,13 +477,13 @@ void ph3py_set_permutation_symmetry_compact_fc3(double * fc3,
 }
 
 void ph3py_transpose_compact_fc3(double * fc3,
-                                 const int p2s[],
+                                 const long p2s[],
-                                 const int s2pp[],
+                                 const long s2pp[],
-                                 const int nsym_list[],
+                                 const long nsym_list[],
-                                 const int perms[],
+                                 const long perms[],
                                  const long n_satom,
                                  const long n_patom,
-                                 const int t_type)
+                                 const long t_type)
 {
   fc3_transpose_compact_fc3(fc3,
                             p2s,
@@ -500,13 +498,13 @@ void ph3py_transpose_compact_fc3(double * fc3,
 
 long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                              long *map_q,
-                                             int (*grid_address)[3],
+                                             long (*grid_address)[3],
                                              const long grid_point,
-                                             const int mesh[3],
+                                             const long mesh[3],
-                                             const int is_time_reversal,
+                                             const long is_time_reversal,
                                              const long num_rot,
-                                             PHPYCONST int (*rotations)[3][3],
+                                             PHPYCONST long (*rotations)[3][3],
-                                             const int swappable)
+                                             const long swappable)
 {
   return tpl_get_triplets_reciprocal_mesh_at_q(map_triplets,
                                                map_q,
@@ -522,11 +520,11 @@ long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
 
 long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
                                 const long grid_point,
-                                PHPYCONST int (*bz_grid_address)[3],
+                                PHPYCONST long (*bz_grid_address)[3],
                                 const long *bz_map,
                                 const long *map_triplets,
                                 const long num_map_triplets,
-                                const int mesh[3])
+                                const long mesh[3])
 {
   return tpl_get_BZ_triplets_at_q(triplets,
                                   grid_point,
@@ -542,19 +540,19 @@ void ph3py_get_integration_weight(double *iw,
                                   char *iw_zero,
                                   const double *frequency_points,
                                   const long num_band0,
-                                  PHPYCONST int relative_grid_address[24][4][3],
+                                  PHPYCONST long relative_grid_address[24][4][3],
-                                  const int mesh[3],
+                                  const long mesh[3],
                                   PHPYCONST long (*triplets)[3],
                                   const long num_triplets,
-                                  PHPYCONST int (*bz_grid_address)[3],
+                                  PHPYCONST long (*bz_grid_address)[3],
                                   const long *bz_map,
                                   const double *frequencies1,
                                   const long num_band1,
                                   const double *frequencies2,
                                   const long num_band2,
                                   const long tp_type,
-                                  const int openmp_per_triplets,
+                                  const long openmp_per_triplets,
-                                  const int openmp_per_bands)
+                                  const long openmp_per_bands)
 {
   tpl_get_integration_weight(iw,
                              iw_zero,
@@ -700,9 +698,9 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
 
 void ph3py_get_neighboring_gird_points(long *relative_grid_points,
                                        const long *grid_points,
-                                       PHPYCONST int (*relative_grid_address)[3],
+                                       PHPYCONST long (*relative_grid_address)[3],
-                                       const int mesh[3],
+                                       const long mesh[3],
-                                       PHPYCONST int (*bz_grid_address)[3],
+                                       PHPYCONST long (*bz_grid_address)[3],
                                        const long *bz_map,
                                        const long num_grid_points,
                                        const long num_relative_grid_address)
@@ -728,10 +726,10 @@ void ph3py_set_integration_weights(double *iw,
                                    const long num_band0,
                                    const long num_band,
                                    const long num_gp,
-                                   PHPYCONST int (*relative_grid_address)[4][3],
+                                   PHPYCONST long (*relative_grid_address)[4][3],
-                                   const int mesh[3],
+                                   const long mesh[3],
                                    const long *grid_points,
-                                   PHPYCONST int (*bz_grid_address)[3],
+                                   PHPYCONST long (*bz_grid_address)[3],
                                    const long *bz_map,
                                    const double *frequencies)
 {

c/phono3py.h

@@ -48,41 +48,41 @@ void ph3py_get_interaction(Darray *fc3_normal_squared,
                            const lapack_complex_double *eigenvectors,
                            const long (*triplets)[3],
                            const long num_triplets,
-                           const int *grid_address,
+                           const long *grid_address,
-                           const int *mesh,
+                           const long *mesh,
                            const double *fc3,
-                           const int is_compact_fc3,
+                           const long is_compact_fc3,
                            const double *shortest_vectors,
-                           const int svecs_dims[3],
+                           const long svecs_dims[3],
-                           const int *multiplicity,
+                           const long *multiplicity,
                            const double *masses,
-                           const int *p2s_map,
+                           const long *p2s_map,
-                           const int *s2p_map,
+                           const long *s2p_map,
-                           const int *band_indices,
+                           const long *band_indices,
-                           const int symmetrize_fc3_q,
+                           const long symmetrize_fc3_q,
                            const double cutoff_frequency);
 void ph3py_get_pp_collision(double *imag_self_energy,
-                            PHPYCONST int relative_grid_address[24][4][3], /* thm */
+                            PHPYCONST long relative_grid_address[24][4][3], /* thm */
                             const double *frequencies,
                             const lapack_complex_double *eigenvectors,
                             const long (*triplets)[3],
                             const long num_triplets,
                             const long *triplet_weights,
-                            const int *grid_address, /* thm */
+                            const long *grid_address, /* thm */
                             const long *bz_map, /* thm */
-                            const int *mesh, /* thm */
+                            const long *mesh, /* thm */
                             const double *fc3,
-                            const int is_compact_fc3,
+                            const long is_compact_fc3,
                             const double *shortest_vectors,
-                            const int svecs_dims[3],
+                            const long svecs_dims[3],
-                            const int *multiplicity,
+                            const long *multiplicity,
                             const double *masses,
-                            const int *p2s_map,
+                            const long *p2s_map,
-                            const int *s2p_map,
+                            const long *s2p_map,
-                            const Iarray *band_indices,
+                            const Larray *band_indices,
                             const Darray *temperatures,
-                            const int is_NU,
+                            const long is_NU,
-                            const int symmetrize_fc3_q,
+                            const long symmetrize_fc3_q,
                             const double cutoff_frequency);
 void ph3py_get_pp_collision_with_sigma(
   double *imag_self_energy,
@@ -93,20 +93,20 @@ void ph3py_get_pp_collision_with_sigma(
   const long (*triplets)[3],
   const long num_triplets,
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
-  const int *mesh,
+  const long *mesh,
   const double *fc3,
-  const int is_compact_fc3,
+  const long is_compact_fc3,
   const double *shortest_vectors,
-  const int svecs_dims[3],
+  const long svecs_dims[3],
-  const int *multiplicity,
+  const long *multiplicity,
   const double *masses,
-  const int *p2s_map,
+  const long *p2s_map,
-  const int *s2p_map,
+  const long *s2p_map,
-  const Iarray *band_indices,
+  const Larray *band_indices,
   const Darray *temperatures,
-  const int is_NU,
+  const long is_NU,
-  const int symmetrize_fc3_q,
+  const long symmetrize_fc3_q,
   const double cutoff_frequency);
 void ph3py_get_imag_self_energy_at_bands_with_g(
   double *imag_self_energy,
@@ -118,8 +118,8 @@ void ph3py_get_imag_self_energy_at_bands_with_g(
   const char *g_zero,
   const double temperature,
   const double cutoff_frequency,
-  const int num_frequency_points,
+  const long num_frequency_points,
-  const int frequency_point_index);
+  const long frequency_point_index);
 void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
   double *detailed_imag_self_energy,
   double *imag_self_energy_N,
@@ -128,14 +128,14 @@ void ph3py_get_detailed_imag_self_energy_at_bands_with_g(
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
   const double *g,
   const char *g_zero,
   const double temperature,
   const double cutoff_frequency);
 void ph3py_get_real_self_energy_at_bands(double *real_self_energy,
                                          const Darray *fc3_normal_squared,
-                                         const int *band_indices,
+                                         const long *band_indices,
                                          const double *frequencies,
                                          const long (*triplets)[3],
                                          const long *triplet_weights,
@@ -147,7 +147,7 @@ void ph3py_get_real_self_energy_at_frequency_point(
   double *real_self_energy,
   const double frequency_point,
   const Darray *fc3_normal_squared,
-  const int *band_indices,
+  const long *band_indices,
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,
@@ -188,7 +188,7 @@ void ph3py_get_isotope_scattering_strength(
   const double *frequencies,
   const lapack_complex_double *eigenvectors,
   const long num_grid_points,
-  const int *band_indices,
+  const long *band_indices,
   const long num_band,
   const long num_band0,
   const double sigma,
@@ -202,74 +202,74 @@ void ph3py_get_thm_isotope_scattering_strength(
   const double *frequencies,
   const lapack_complex_double *eigenvectors,
   const long num_ir_grid_points,
-  const int *band_indices,
+  const long *band_indices,
   const long num_band,
   const long num_band0,
   const double *integration_weights,
   const double cutoff_frequency);
 void ph3py_distribute_fc3(double *fc3,
-                          const int target,
+                          const long target,
-                          const int source,
+                          const long source,
-                          const int *atom_mapping,
+                          const long *atom_mapping,
                           const long num_atom,
                           const double *rot_cart);
 void ph3py_rotate_delta_fc2(double (*fc3)[3][3][3],
                             PHPYCONST double (*delta_fc2s)[3][3],
                             const double *inv_U,
                             PHPYCONST double (*site_sym_cart)[3][3],
-                            const int *rot_map_syms,
+                            const long *rot_map_syms,
                             const long num_atom,
                             const long num_site_sym,
                             const long num_disp);
 void ph3py_set_permutation_symmetry_fc3(double *fc3, const long num_atom);
 void ph3py_set_permutation_symmetry_compact_fc3(double * fc3,
-                                                const int p2s[],
+                                                const long p2s[],
-                                                const int s2pp[],
+                                                const long s2pp[],
-                                                const int nsym_list[],
+                                                const long nsym_list[],
-                                                const int perms[],
+                                                const long perms[],
                                                 const long n_satom,
                                                 const long n_patom);
 void ph3py_transpose_compact_fc3(double * fc3,
-                                 const int p2s[],
+                                 const long p2s[],
-                                 const int s2pp[],
+                                 const long s2pp[],
-                                 const int nsym_list[],
+                                 const long nsym_list[],
-                                 const int perms[],
+                                 const long perms[],
                                  const long n_satom,
                                  const long n_patom,
-                                 const int t_type);
+                                 const long t_type);
 long ph3py_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                              long *map_q,
-                                             int (*grid_address)[3],
+                                             long (*grid_address)[3],
                                              const long grid_point,
-                                             const int mesh[3],
+                                             const long mesh[3],
-                                             const int is_time_reversal,
+                                             const long is_time_reversal,
                                              const long num_rot,
-                                             PHPYCONST int (*rotations)[3][3],
+                                             PHPYCONST long (*rotations)[3][3],
-                                             const int swappable);
+                                             const long swappable);
 long ph3py_get_BZ_triplets_at_q(long (*triplets)[3],
                                 const long grid_point,
-                                PHPYCONST int (*bz_grid_address)[3],
+                                PHPYCONST long (*bz_grid_address)[3],
                                 const long *bz_map,
                                 const long *map_triplets,
                                 const long num_map_triplets,
-                                const int mesh[3]);
+                                const long mesh[3]);
 void ph3py_get_integration_weight(double *iw,
                                   char *iw_zero,
                                   const double *frequency_points,
                                   const long num_band0,
-                                  PHPYCONST int relative_grid_address[24][4][3],
+                                  PHPYCONST long relative_grid_address[24][4][3],
-                                  const int mesh[3],
+                                  const long mesh[3],
                                   PHPYCONST long (*triplets)[3],
                                   const long num_triplets,
-                                  PHPYCONST int (*bz_grid_address)[3],
+                                  PHPYCONST long (*bz_grid_address)[3],
                                   const long *bz_map,
                                   const double *frequencies1,
                                   const long num_band1,
                                   const double *frequencies2,
                                   const long num_band2,
                                   const long tp_type,
-                                  const int openmp_per_triplets,
+                                  const long openmp_per_triplets,
-                                  const int openmp_per_bands);
+                                  const long openmp_per_bands);
 void ph3py_get_integration_weight_with_sigma(double *iw,
                                              char *iw_zero,
                                              const double sigma,
@@ -298,9 +298,9 @@ void ph3py_expand_collision_matrix(double *collision_matrix,
                                    const long num_band);
 void ph3py_get_neighboring_gird_points(long *relative_grid_points,
                                        const long *grid_points,
-                                       PHPYCONST int (*relative_grid_address)[3],
+                                       PHPYCONST long (*relative_grid_address)[3],
-                                       const int mesh[3],
+                                       const long mesh[3],
-                                       PHPYCONST int (*bz_grid_address)[3],
+                                       PHPYCONST long (*bz_grid_address)[3],
                                        const long *bz_map,
                                        const long num_grid_points,
                                        const long num_relative_grid_address);
@@ -309,10 +309,10 @@ void ph3py_set_integration_weights(double *iw,
                                    const long num_band0,
                                    const long num_band,
                                    const long num_gp,
-                                   PHPYCONST int (*relative_grid_address)[4][3],
+                                   PHPYCONST long (*relative_grid_address)[4][3],
-                                   const int mesh[3],
+                                   const long mesh[3],
                                    const long *grid_points,
-                                   PHPYCONST int (*bz_grid_address)[3],
+                                   PHPYCONST long (*bz_grid_address)[3],
                                    const long *bz_map,
                                    const double *frequencies);
 

c/phonoc_array.h

@@ -39,9 +39,9 @@
 
 /* It is assumed that number of dimensions is known for each array. */
 typedef struct {
-  int dims[MAX_NUM_DIM];
+  long dims[MAX_NUM_DIM];
-  int *data;
+  long *data;
-} Iarray;
+} Larray;
 
 typedef struct {
   int dims[MAX_NUM_DIM];

c/pp_collision.c

@@ -55,59 +55,59 @@ static void get_collision(double *ise,
                           const lapack_complex_double *eigenvectors,
                           const long triplet[3],
                           const long triplet_weight,
-                          const int *grid_address,
+                          const long *grid_address,
-                          const int *mesh,
+                          const long *mesh,
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
                           const double *masses,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map,
+                          const long *s2p_map,
-                          const int *band_indices,
+                          const long *band_indices,
-                          const int symmetrize_fc3_q,
+                          const long symmetrize_fc3_q,
                           const double cutoff_frequency,
-                          const int openmp_per_triplets);
+                          const long openmp_per_triplets);
 static void finalize_ise(double *imag_self_energy,
                          const double *ise,
-                         const int *grid_address,
+                         const long *grid_address,
                          const long (*triplets)[3],
                          const long num_triplets,
                          const long num_temps,
                          const long num_band0,
-                         const int is_NU);
+                         const long is_NU);
 
 void ppc_get_pp_collision(double *imag_self_energy,
-                          PHPYCONST int relative_grid_address[24][4][3], /* thm */
+                          PHPYCONST long relative_grid_address[24][4][3], /* thm */
                           const double *frequencies,
                           const lapack_complex_double *eigenvectors,
                           const long (*triplets)[3],
                           const long num_triplets,
                           const long *triplet_weights,
-                          const int *grid_address, /* thm */
+                          const long *grid_address, /* thm */
                           const long *bz_map, /* thm */
-                          const int *mesh, /* thm */
+                          const long *mesh, /* thm */
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
                           const double *masses,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map,
+                          const long *s2p_map,
-                          const Iarray *band_indices,
+                          const Larray *band_indices,
                           const Darray *temperatures,
-                          const int is_NU,
+                          const long is_NU,
-                          const int symmetrize_fc3_q,
+                          const long symmetrize_fc3_q,
                           const double cutoff_frequency)
 {
   long i;
   long num_band, num_band0, num_band_prod, num_temps;
-  int openmp_per_triplets;
+  long openmp_per_triplets;
   double *ise, *freqs_at_gp, *g;
   char *g_zero;
-  int tp_relative_grid_address[2][24][4][3];
+  long tp_relative_grid_address[2][24][4][3];
 
   ise = NULL;
   freqs_at_gp = NULL;
@@ -147,7 +147,7 @@ void ppc_get_pp_collision(double *imag_self_energy,
                                mesh,
                                triplets[i],
                                1,
-                               (int(*)[3])grid_address,
+                               (long(*)[3])grid_address,
                                bz_map,
                                frequencies,  /* used as f1 */
                                num_band,
@@ -212,25 +212,25 @@ void ppc_get_pp_collision_with_sigma(
   const long (*triplets)[3],
   const long num_triplets,
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
-  const int *mesh,
+  const long *mesh,
   const double *fc3,
-  const int is_compact_fc3,
+  const long is_compact_fc3,
   const double *shortest_vectors,
-  const int svecs_dims[3],
+  const long svecs_dims[3],
-  const int *multiplicity,
+  const long *multiplicity,
   const double *masses,
-  const int *p2s_map,
+  const long *p2s_map,
-  const int *s2p_map,
+  const long *s2p_map,
-  const Iarray *band_indices,
+  const Larray *band_indices,
   const Darray *temperatures,
-  const int is_NU,
+  const long is_NU,
-  const int symmetrize_fc3_q,
+  const long symmetrize_fc3_q,
   const double cutoff_frequency)
 {
   long i;
   long num_band, num_band0, num_band_prod, num_temps;
-  int openmp_per_triplets, const_adrs_shift;
+  long openmp_per_triplets, const_adrs_shift;
   double cutoff;
   double *ise, *freqs_at_gp, *g;
   char *g_zero;
@@ -336,32 +336,32 @@ static void get_collision(double *ise,
                           const lapack_complex_double *eigenvectors,
                           const long triplet[3],
                           const long triplet_weight,
-                          const int *grid_address,
+                          const long *grid_address,
-                          const int *mesh,
+                          const long *mesh,
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
                           const double *masses,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map,
+                          const long *s2p_map,
-                          const int *band_indices,
+                          const long *band_indices,
-                          const int symmetrize_fc3_q,
+                          const long symmetrize_fc3_q,
                           const double cutoff_frequency,
-                          const int openmp_per_triplets)
+                          const long openmp_per_triplets)
 {
   long i;
   long num_band_prod, num_g_pos;
   double *fc3_normal_squared;
-  int (*g_pos)[4];
+  long (*g_pos)[4];
 
   fc3_normal_squared = NULL;
   g_pos = NULL;
 
   num_band_prod = num_band0 * num_band * num_band;
   fc3_normal_squared = (double*)malloc(sizeof(double) * num_band_prod);
-  g_pos = (int(*)[4])malloc(sizeof(int[4]) * num_band_prod);
+  g_pos = (long(*)[4])malloc(sizeof(long[4]) * num_band_prod);
 
   for (i = 0; i < num_band_prod; i++) {
     fc3_normal_squared[i] = 0;
@@ -424,15 +424,15 @@ static void get_collision(double *ise,
 
 static void finalize_ise(double *imag_self_energy,
                          const double *ise,
-                         const int *grid_address,
+                         const long *grid_address,
                          const long (*triplets)[3],
                          const long num_triplets,
                          const long num_temps,
                          const long num_band0,
-                         const int is_NU)
+                         const long is_NU)
 {
   long i, j, k;
-  int is_N;
+  long is_N;
 
   if (is_NU) {
     for (i = 0; i < 2 * num_temps * num_band0; i++) {

c/pp_collision.h

@@ -40,27 +40,27 @@
 #include "lapack_wrapper.h"
 
 void ppc_get_pp_collision(double *imag_self_energy,
-                          PHPYCONST int relative_grid_address[24][4][3],
+                          PHPYCONST long relative_grid_address[24][4][3],
                           const double *frequencies,
                           const lapack_complex_double *eigenvectors,
                           const long (*triplets)[3],
                           const long num_triplets,
                           const long *triplet_weights,
-                          const int *grid_address,
+                          const long *grid_address,
                           const long *bz_map,
-                          const int *mesh,
+                          const long *mesh,
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
                           const double *masses,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map,
+                          const long *s2p_map,
-                          const Iarray *band_indices,
+                          const Larray *band_indices,
                           const Darray *temperatures,
-                          const int is_NU,
+                          const long is_NU,
-                          const int symmetrize_fc3_q,
+                          const long symmetrize_fc3_q,
                           const double cutoff_frequency);
 
 void ppc_get_pp_collision_with_sigma(
@@ -72,20 +72,20 @@ void ppc_get_pp_collision_with_sigma(
   const long (*triplets)[3],
   const long num_triplets,
   const long *triplet_weights,
-  const int *grid_address,
+  const long *grid_address,
-  const int *mesh,
+  const long *mesh,
   const double *fc3,
-  const int is_compact_fc3,
+  const long is_compact_fc3,
   const double *shortest_vectors,
-  const int svecs_dims[3],
+  const long svecs_dims[3],
-  const int *multiplicity,
+  const long *multiplicity,
   const double *masses,
-  const int *p2s_map,
+  const long *p2s_map,
-  const int *s2p_map,
+  const long *s2p_map,
-  const Iarray *band_indices,
+  const Larray *band_indices,
   const Darray *temperatures,
-  const int is_NU,
+  const long is_NU,
-  const int symmetrize_fc3_q,
+  const long symmetrize_fc3_q,
   const double cutoff_frequency);
 
 #endif

c/real_self_energy.c

@@ -39,7 +39,7 @@
 #include "real_self_energy.h"
 #include "real_to_reciprocal.h"
 
-static double get_real_self_energy_at_band(const int band_index,
+static double get_real_self_energy_at_band(const long band_index,
                                            const Darray *fc3_normal_squared,
                                            const double fpoint,
                                            const double *frequencies,
@@ -49,7 +49,7 @@ static double get_real_self_energy_at_band(const int band_index,
                                            const double temperature,
                                            const double unit_conversion_factor,
                                            const double cutoff_frequency);
-static double sum_real_self_energy_at_band(const int num_band,
+static double sum_real_self_energy_at_band(const long num_band,
                                            const double *fc3_normal_squared,
                                            const double fpoint,
                                            const double *freqs1,
@@ -57,7 +57,7 @@ static double sum_real_self_energy_at_band(const int num_band,
                                            const double epsilon,
                                            const double temperature,
                                            const double cutoff_frequency);
-static double sum_real_self_energy_at_band_0K(const int num_band,
+static double sum_real_self_energy_at_band_0K(const long num_band,
                                               const double *fc3_normal_squared,
                                               const double fpoint,
                                               const double *freqs1,
@@ -67,7 +67,7 @@ static double sum_real_self_energy_at_band_0K(const int num_band,
 
 void rse_get_real_self_energy_at_bands(double *real_self_energy,
                                        const Darray *fc3_normal_squared,
-                                       const int *band_indices,
+                                       const long *band_indices,
                                        const double *frequencies,
                                        const long (*triplets)[3],
                                        const long *triplet_weights,
@@ -76,7 +76,7 @@ void rse_get_real_self_energy_at_bands(double *real_self_energy,
                                        const double unit_conversion_factor,
                                        const double cutoff_frequency)
 {
-  int i, num_band0, num_band, gp0;
+  long i, num_band0, num_band, gp0;
   double fpoint;
 
   num_band0 = fc3_normal_squared->dims[1];
@@ -108,7 +108,7 @@ void rse_get_real_self_energy_at_frequency_point(
   double *real_self_energy,
   const double frequency_point,
   const Darray *fc3_normal_squared,
-  const int *band_indices,
+  const long *band_indices,
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,
@@ -117,7 +117,7 @@ void rse_get_real_self_energy_at_frequency_point(
   const double unit_conversion_factor,
   const double cutoff_frequency)
 {
-  int i, num_band0;
+  long i, num_band0;
 
   num_band0 = fc3_normal_squared->dims[1];
 
@@ -141,7 +141,7 @@ void rse_get_real_self_energy_at_frequency_point(
   }
 }
 
-static double get_real_self_energy_at_band(const int band_index,
+static double get_real_self_energy_at_band(const long band_index,
                                            const Darray *fc3_normal_squared,
                                            const double fpoint,
                                            const double *frequencies,
@@ -152,7 +152,7 @@ static double get_real_self_energy_at_band(const int band_index,
                                            const double unit_conversion_factor,
                                            const double cutoff_frequency)
 {
-  int i, num_triplets, num_band0, num_band, gp1, gp2;
+  long i, num_triplets, num_band0, num_band, gp1, gp2;
   double shift;
 
   num_triplets = fc3_normal_squared->dims[0];
@@ -194,7 +194,7 @@ static double get_real_self_energy_at_band(const int band_index,
   return shift;
 }
 
-static double sum_real_self_energy_at_band(const int num_band,
+static double sum_real_self_energy_at_band(const long num_band,
                                            const double *fc3_normal_squared,
                                            const double fpoint,
                                            const double *freqs1,
@@ -203,7 +203,7 @@ static double sum_real_self_energy_at_band(const int num_band,
                                            const double temperature,
                                            const double cutoff_frequency)
 {
-  int i, j;
+  long i, j;
   double n1, n2, f1, f2, f3, f4, shift;
   /* double sum; */
 
@@ -246,7 +246,7 @@ static double sum_real_self_energy_at_band(const int num_band,
   return shift;
 }
 
-static double sum_real_self_energy_at_band_0K(const int num_band,
+static double sum_real_self_energy_at_band_0K(const long num_band,
                                               const double *fc3_normal_squared,
                                               const double fpoint,
                                               const double *freqs1,
@@ -254,7 +254,7 @@ static double sum_real_self_energy_at_band_0K(const int num_band,
                                               const double epsilon,
                                               const double cutoff_frequency)
 {
-  int i, j;
+  long i, j;
   double f1, f2, shift;
 
   shift = 0;

c/real_self_energy.h

@@ -41,7 +41,7 @@
 
 void rse_get_real_self_energy_at_bands(double *real_self_energy,
                                        const Darray *fc3_normal_squared,
-                                       const int *band_indices,
+                                       const long *band_indices,
                                        const double *frequencies,
                                        const long (*triplets)[3],
                                        const long *triplet_weights,
@@ -53,7 +53,7 @@ void rse_get_real_self_energy_at_frequency_point(
   double *real_self_energy,
   const double frequency_point,
   const Darray *fc3_normal_squared,
-  const int *band_indices,
+  const long *band_indices,
   const double *frequencies,
   const long (*triplets)[3],
   const long *triplet_weights,

c/real_to_reciprocal.c

@@ -46,57 +46,57 @@ static void
 real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
                                  const double q[9],
                                  const double *fc3,
-                                 const int is_compact_fc3,
+                                 const long is_compact_fc3,
                                  const double *shortest_vectors,
-                                 const int svecs_dims[3],
+                                 const long svecs_dims[3],
-                                 const int *multiplicity,
+                                 const long *multiplicity,
-                                 const int *p2s_map,
+                                 const long *p2s_map,
-                                 const int *s2p_map);
+                                 const long *s2p_map);
 static void
 real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
                           const double q[9],
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map);
+                          const long *s2p_map);
 static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
                                         const double q[9],
                                         const double *fc3,
-                                        const int is_compact_fc3,
+                                        const long is_compact_fc3,
                                         const double *shortest_vectors,
-                                        const int svecs_dims[3],
+                                        const long svecs_dims[3],
-                                        const int *multiplicity,
+                                        const long *multiplicity,
-                                        const int *p2s,
+                                        const long *p2s,
-                                        const int *s2p,
+                                        const long *s2p,
                                         const long pi0,
                                         const long pi1,
                                         const long pi2);
 static lapack_complex_double get_phase_factor(const double q[],
-                                              const int qi,
+                                              const long qi,
                                               const double *shortest_vectors,
-                                              const int multi);
+                                              const long multi);
 static lapack_complex_double
 get_pre_phase_factor(const long i,
                      const double q[9],
                      const double *shortest_vectors,
-                     const int svecs_dims[3],
+                     const long svecs_dims[3],
-                     const int *multiplicity,
+                     const long *multiplicity,
-                     const int *p2s_map);
+                     const long *p2s_map);
 
 /* fc3_reciprocal[num_patom, num_patom, num_patom, 3, 3, 3] */
 void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
                             const double q[9],
                             const double *fc3,
-                            const int is_compact_fc3,
+                            const long is_compact_fc3,
                             const double *shortest_vectors,
-                            const int svecs_dims[3],
+                            const long svecs_dims[3],
-                            const int *multiplicity,
+                            const long *multiplicity,
-                            const int *p2s_map,
+                            const long *p2s_map,
-                            const int *s2p_map,
+                            const long *s2p_map,
-                            const int openmp_at_bands)
+                            const long openmp_at_bands)
 {
   if (openmp_at_bands) {
     real_to_reciprocal_openmp(fc3_reciprocal,
@@ -126,12 +126,12 @@ static void
 real_to_reciprocal_single_thread(lapack_complex_double *fc3_reciprocal,
                                  const double q[9],
                                  const double *fc3,
-                                 const int is_compact_fc3,
+                                 const long is_compact_fc3,
                                  const double *shortest_vectors,
-                                 const int svecs_dims[3],
+                                 const long svecs_dims[3],
-                                 const int *multiplicity,
+                                 const long *multiplicity,
-                                 const int *p2s_map,
+                                 const long *p2s_map,
-                                 const int *s2p_map)
+                                 const long *s2p_map)
 {
   long i, j, k;
   long num_patom, adrs_shift;
@@ -172,12 +172,12 @@ static void
 real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
                           const double q[9],
                           const double *fc3,
-                          const int is_compact_fc3,
+                          const long is_compact_fc3,
                           const double *shortest_vectors,
-                          const int svecs_dims[3],
+                          const long svecs_dims[3],
-                          const int *multiplicity,
+                          const long *multiplicity,
-                          const int *p2s_map,
+                          const long *p2s_map,
-                          const int *s2p_map)
+                          const long *s2p_map)
 {
   long i, j, k, jk;
   long num_patom, adrs_shift;
@@ -219,12 +219,12 @@ real_to_reciprocal_openmp(lapack_complex_double *fc3_reciprocal,
 static void real_to_reciprocal_elements(lapack_complex_double *fc3_rec_elem,
                                         const double q[9],
                                         const double *fc3,
-                                        const int is_compact_fc3,
+                                        const long is_compact_fc3,
                                         const double *shortest_vectors,
-                                        const int svecs_dims[3],
+                                        const long svecs_dims[3],
-                                        const int *multiplicity,
+                                        const long *multiplicity,
-                                        const int *p2s,
+                                        const long *p2s,
-                                        const int *s2p,
+                                        const long *s2p,
                                         const long pi0,
                                         const long pi1,
                                         const long pi2)
@@ -289,11 +289,11 @@ static lapack_complex_double
 get_pre_phase_factor(const long i,
                      const double q[9],
                      const double *shortest_vectors,
-                     const int svecs_dims[3],
+                     const long svecs_dims[3],
-                     const int *multiplicity,
+                     const long *multiplicity,
-                     const int *p2s_map)
+                     const long *p2s_map)
 {
-  int j;
+  long j;
   double pre_phase;
   lapack_complex_double pre_phase_factor;
 
@@ -313,11 +313,11 @@ get_pre_phase_factor(const long i,
 }
 
 static lapack_complex_double get_phase_factor(const double q[],
-                                              const int qi,
+                                              const long qi,
                                               const double *shortest_vectors,
-                                              const int multi)
+                                              const long multi)
 {
-  int i, j;
+  long i, j;
   double sum_real, sum_imag, phase;
 
   sum_real = 0;

c/real_to_reciprocal.h

@@ -41,11 +41,11 @@
 void r2r_real_to_reciprocal(lapack_complex_double *fc3_reciprocal,
                             const double q[9],
                             const double *fc3,
-                            const int is_compact_fc3,
+                            const long is_compact_fc3,
                             const double *shortest_vectors,
-                            const int svecs_dims[3],
+                            const long svecs_dims[3],
-                            const int *multiplicity,
+                            const long *multiplicity,
-                            const int *p2s_map,
+                            const long *p2s_map,
-                            const int *s2p_map,
+                            const long *s2p_map,
-                            const int openmp_at_bands);
+                            const long openmp_at_bands);
 #endif

c/reciprocal_to_normal.c

@@ -73,7 +73,7 @@ static double get_fc3_sum
 
 void reciprocal_to_normal_squared
 (double *fc3_normal_squared,
- PHPYCONST int (*g_pos)[4],
+ PHPYCONST long (*g_pos)[4],
  const long num_g_pos,
  const lapack_complex_double *fc3_reciprocal,
  const double *freqs0,
@@ -83,11 +83,11 @@ void reciprocal_to_normal_squared
  const lapack_complex_double *eigvecs1,
  const lapack_complex_double *eigvecs2,
  const double *masses,
- const int *band_indices,
+ const long *band_indices,
  const long num_band0,
  const long num_band,
  const double cutoff_frequency,
- const int openmp_at_bands)
+ const long openmp_at_bands)
 {
   long i, num_atom;
 
@@ -124,9 +124,9 @@ void reciprocal_to_normal_squared
   }
 
 #ifdef MEASURE_R2N
-      loopTotalCPUTime = (double)(clock() - loopStartCPUTime) / CLOCKS_PER_SEC;
+  loopTotalCPUTime = (double)(clock() - loopStartCPUTime) / CLOCKS_PER_SEC;
-      loopTotalWallTime = difftime(time(NULL), loopStartWallTime);
+  loopTotalWallTime = difftime(time(NULL), loopStartWallTime);
-      printf("  %1.3fs (%1.3fs CPU)\n", loopTotalWallTime, loopTotalCPUTime);
+  printf("  %1.3fs (%1.3fs CPU)\n", loopTotalWallTime, loopTotalCPUTime);
 #endif
 
 }

c/reciprocal_to_normal.h

@@ -41,7 +41,7 @@
 
 void reciprocal_to_normal_squared
 (double *fc3_normal_squared,
- PHPYCONST int (*g_pos)[4],
+ PHPYCONST long (*g_pos)[4],
  const long num_g_pos,
  const lapack_complex_double *fc3_reciprocal,
  const double *freqs0,
@@ -51,10 +51,10 @@ void reciprocal_to_normal_squared
  const lapack_complex_double *eigvecs1,
  const lapack_complex_double *eigvecs2,
  const double *masses,
- const int *band_indices,
+ const long *band_indices,
  const long num_band0,
  const long num_band,
  const double cutoff_frequency,
- const int openmp_at_bands);
+ const long openmp_at_bands);
 
 #endif

c/triplet.c

@@ -34,29 +34,28 @@
 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
 /* POSSIBILITY OF SUCH DAMAGE. */
 
-#include <stddef.h>
+#include "kpoint.h"
-#include <mathfunc.h>
 #include "triplet.h"
 #include "triplet_iw.h"
 #include "triplet_kpoint.h"
 
 static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                               long *map_q,
-                                              int (*grid_address)[3],
+                                              long (*grid_address)[3],
-                                              const int grid_point,
+                                              const long grid_point,
-                                              const int mesh[3],
+                                              const long mesh[3],
-                                              const int is_time_reversal,
+                                              const long is_time_reversal,
                                               const long num_rot,
-                                              TPLCONST int (*rotations)[3][3],
+                                              TPLCONST long (*rotations)[3][3],
-                                              const int swappable);
+                                              const long swappable);
 
 long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
                               const long grid_point,
-                              TPLCONST int (*bz_grid_address)[3],
+                              TPLCONST long (*bz_grid_address)[3],
                               const long *bz_map,
                               const long *map_triplets,
                               const long num_map_triplets,
-                              const int mesh[3])
+                              const long mesh[3])
 {
   return tpk_get_BZ_triplets_at_q(triplets,
                                   grid_point,
@@ -69,13 +68,13 @@ long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
 
 long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                            long *map_q,
-                                           int (*grid_address)[3],
+                                           long (*grid_address)[3],
                                            const long grid_point,
-                                           const int mesh[3],
+                                           const long mesh[3],
-                                           const int is_time_reversal,
+                                           const long is_time_reversal,
                                            const long num_rot,
-                                           TPLCONST int (*rotations)[3][3],
+                                           TPLCONST long (*rotations)[3][3],
-                                           const int swappable)
+                                           const long swappable)
 {
   return get_triplets_reciprocal_mesh_at_q(map_triplets,
                                            map_q,
@@ -92,22 +91,22 @@ void tpl_get_integration_weight(double *iw,
                                 char *iw_zero,
                                 const double *frequency_points,
                                 const long num_band0,
-                                TPLCONST int relative_grid_address[24][4][3],
+                                TPLCONST long relative_grid_address[24][4][3],
-                                const int mesh[3],
+                                const long mesh[3],
                                 TPLCONST long (*triplets)[3],
                                 const long num_triplets,
-                                TPLCONST int (*bz_grid_address)[3],
+                                TPLCONST long (*bz_grid_address)[3],
                                 const long *bz_map,
                                 const double *frequencies1,
                                 const long num_band1,
                                 const double *frequencies2,
                                 const long num_band2,
                                 const long tp_type,
-                                const int openmp_per_triplets,
+                                const long openmp_per_triplets,
-                                const int openmp_per_bands)
+                                const long openmp_per_bands)
 {
   long i, num_band_prod;
-  int tp_relative_grid_address[2][24][4][3];
+  long tp_relative_grid_address[2][24][4][3];
 
   tpl_set_relative_grid_address(tp_relative_grid_address,
                                 relative_grid_address,
@@ -174,9 +173,9 @@ void tpl_get_integration_weight_with_sigma(double *iw,
 }
 
 
-int tpl_is_N(const long triplet[3], const int *grid_address)
+long tpl_is_N(const long triplet[3], const long *grid_address)
 {
-  int i, j, sum_q, is_N;
+  long i, j, sum_q, is_N;
 
   is_N = 1;
   for (i = 0; i < 3; i++) {
@@ -193,12 +192,12 @@ int tpl_is_N(const long triplet[3], const int *grid_address)
 }
 
 void tpl_set_relative_grid_address(
-  int tp_relative_grid_address[2][24][4][3],
+  long tp_relative_grid_address[2][24][4][3],
-  TPLCONST int relative_grid_address[24][4][3],
+  TPLCONST long relative_grid_address[24][4][3],
   const long tp_type)
 {
   long i, j, k, l;
-  int signs[2];
+  long signs[2];
 
   signs[0] = 1;
   signs[1] = 1;
@@ -223,20 +222,20 @@ void tpl_set_relative_grid_address(
 
 static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                               long *map_q,
-                                              int (*grid_address)[3],
+                                              long (*grid_address)[3],
-                                              const int grid_point,
+                                              const long grid_point,
-                                              const int mesh[3],
+                                              const long mesh[3],
-                                              const int is_time_reversal,
+                                              const long is_time_reversal,
                                               const long num_rot,
-                                              TPLCONST int (*rotations)[3][3],
+                                              TPLCONST long (*rotations)[3][3],
-                                              const int swappable)
+                                              const long swappable)
 {
-  MatINT *rot_real;
+  MatLONG *rot_real;
   long i, num_ir;
 
-  rot_real = mat_alloc_MatINT(num_rot);
+  rot_real = kpt_alloc_MatLONG(num_rot);
   for (i = 0; i < num_rot; i++) {
-    mat_copy_matrix_i3(rot_real->mat[i], rotations[i]);
+    kpt_copy_matrix_l3(rot_real->mat[i], rotations[i]);
   }
 
   num_ir = tpk_get_ir_triplets_at_q(map_triplets,
@@ -248,7 +247,7 @@ static long get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                     rot_real,
                                     swappable);
 
-  mat_free_MatINT(rot_real);
+  kpt_free_MatLONG(rot_real);
 
   return num_ir;
 }

c/triplet.h

@@ -51,13 +51,13 @@
 /* in the input. */
 long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
                                            long *map_q,
-                                           int (*grid_address)[3],
+                                           long (*grid_address)[3],
                                            const long grid_point,
-                                           const int mesh[3],
+                                           const long mesh[3],
-                                           const int is_time_reversal,
+                                           const long is_time_reversal,
                                            const long num_rot,
-                                           TPLCONST int (*rotations)[3][3],
+                                           TPLCONST long (*rotations)[3][3],
-                                           const int swappable);
+                                           const long swappable);
 /* Irreducible grid-point-triplets in BZ are stored. */
 /* triplets are recovered from grid_point and triplet_weights. */
 /* BZ boundary is considered in this recovery. Therefore grid addresses */
@@ -66,28 +66,28 @@ long tpl_get_triplets_reciprocal_mesh_at_q(long *map_triplets,
 /* Number of ir-triplets is returned. */
 long tpl_get_BZ_triplets_at_q(long (*triplets)[3],
                               const long grid_point,
-                              TPLCONST int (*bz_grid_address)[3],
+                              TPLCONST long (*bz_grid_address)[3],
                               const long *bz_map,
                               const long *map_triplets,
                               const long num_map_triplets,
-                              const int mesh[3]);
+                              const long mesh[3]);
 void tpl_get_integration_weight(double *iw,
                                 char *iw_zero,
                                 const double *frequency_points,
                                 const long num_band0,
-                                TPLCONST int relative_grid_address[24][4][3],
+                                TPLCONST long relative_grid_address[24][4][3],
-                                const int mesh[3],
+                                const long mesh[3],
                                 TPLCONST long (*triplets)[3],
                                 const long num_triplets,
-                                TPLCONST int (*bz_grid_address)[3],
+                                TPLCONST long (*bz_grid_address)[3],
                                 const long *bz_map,
                                 const double *frequencies1,
                                 const long num_band1,
                                 const double *frequencies2,
                                 const long num_band2,
                                 const long tp_type,
-                                const int openmp_per_triplets,
+                                const long openmp_per_triplets,
-                                const int openmp_per_bands);
+                                const long openmp_per_bands);
 void tpl_get_integration_weight_with_sigma(double *iw,
                                            char *iw_zero,
                                            const double sigma,
@@ -100,10 +100,10 @@ void tpl_get_integration_weight_with_sigma(double *iw,
                                            const long num_band,
                                            const long tp_type);
 
-int tpl_is_N(const long triplet[3], const int *grid_address);
+long tpl_is_N(const long triplet[3], const long *grid_address);
 void tpl_set_relative_grid_address(
-  int tp_relative_grid_address[2][24][4][3],
+  long tp_relative_grid_address[2][24][4][3],
-  TPLCONST int relative_grid_address[24][4][3],
+  TPLCONST long relative_grid_address[24][4][3],
   const long tp_type);
 
 #endif

c/triplet_iw.c

@@ -32,9 +32,8 @@
 /* ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE */
 /* POSSIBILITY OF SUCH DAMAGE. */
 
-#include <stddef.h>
 #include <math.h>
-#include "kgrid.h"
+#include "rgrid.h"
 #include "phonoc_utils.h"
 #include "triplet.h"
 #include "triplet_iw.h"
@@ -44,22 +43,22 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
                               const double *frequencies1,
                               const double *frequencies2,
                               TPLCONST long vertices[2][24][4],
-                              const int num_band1,
+                              const long num_band1,
-                              const int num_band2,
+                              const long num_band2,
-                              const int b1,
+                              const long b1,
-                              const int b2,
+                              const long b2,
                               const long tp_type);
-static int set_g(double g[3],
+static long set_g(double g[3],
-                 const double f0,
+                  const double f0,
-                 TPLCONST double freq_vertices[3][24][4],
+                  TPLCONST double freq_vertices[3][24][4],
-                 const long max_i);
+                  const long max_i);
-static int in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4]);
+static long in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4]);
 static void get_triplet_tetrahedra_vertices(
   long vertices[2][24][4],
-  TPLCONST int tp_relative_grid_address[2][24][4][3],
+  TPLCONST long tp_relative_grid_address[2][24][4][3],
-  const int mesh[3],
+  const long mesh[3],
   const long triplet[3],
-  TPLCONST int (*bz_grid_address)[3],
+  TPLCONST long (*bz_grid_address)[3],
   const long *bz_map);
 
 void
@@ -67,18 +66,18 @@ tpi_get_integration_weight(double *iw,
                            char *iw_zero,
                            const double *frequency_points,
                            const long num_band0,
-                           TPLCONST int tp_relative_grid_address[2][24][4][3],
+                           TPLCONST long tp_relative_grid_address[2][24][4][3],
-                           const int mesh[3],
+                           const long mesh[3],
                            const long triplets[3],
                            const long num_triplets,
-                           TPLCONST int (*bz_grid_address)[3],
+                           TPLCONST long (*bz_grid_address)[3],
                            const long *bz_map,
                            const double *frequencies1,
                            const long num_band1,
                            const double *frequencies2,
                            const long num_band2,
                            const long tp_type,
-                           const int openmp_per_bands)
+                           const long openmp_per_bands)
 {
   long max_i, j, b1, b2, b12, num_band_prod, adrs_shift;
   long vertices[2][24][4];
@@ -151,7 +150,7 @@ void tpi_get_integration_weight_with_sigma(double *iw,
                                            const double *frequencies,
                                            const long num_band,
                                            const long tp_type,
-                                           const int openmp_per_bands)
+                                           const long openmp_per_bands)
 {
   long j, b12, b1, b2, adrs_shift;
   double f0, f1, f2, g0, g1, g2;
@@ -212,13 +211,13 @@ void tpi_get_integration_weight_with_sigma(double *iw,
 void
 tpi_get_dense_neighboring_grid_points(long neighboring_grid_points[],
                                       const long grid_point,
-                                      TPLCONST int relative_grid_address[][3],
+                                      TPLCONST long relative_grid_address[][3],
-                                      const int num_relative_grid_address,
+                                      const long num_relative_grid_address,
-                                      const int mesh[3],
+                                      const long mesh[3],
-                                      TPLCONST int bz_grid_address[][3],
+                                      TPLCONST long bz_grid_address[][3],
                                       const long bz_map[])
 {
-  int bzmesh[3], address_double[3], bz_address_double[3];
+  long bzmesh[3], address_double[3], bz_address_double[3];
   long i, j, bz_gp, prod_bz_mesh;
 
   prod_bz_mesh = 1;
@@ -232,10 +231,10 @@ tpi_get_dense_neighboring_grid_points(long neighboring_grid_points[],
                            relative_grid_address[i][j]) * 2;
       bz_address_double[j] = address_double[j];
     }
-    bz_gp = bz_map[kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh)];
+    bz_gp = bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
     if (bz_gp == prod_bz_mesh) {
       neighboring_grid_points[i] =
-        kgd_get_dense_grid_point_double_mesh(address_double, mesh);
+        rgd_get_double_grid_index(address_double, mesh);
     } else {
       neighboring_grid_points[i] = bz_gp;
     }
@@ -246,13 +245,13 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
                               const double *frequencies1,
                               const double *frequencies2,
                               TPLCONST long vertices[2][24][4],
-                              const int num_band1,
+                              const long num_band1,
-                              const int num_band2,
+                              const long num_band2,
-                              const int b1,
+                              const long b1,
-                              const int b2,
+                              const long b2,
                               const long tp_type)
 {
-  int i, j;
+  long i, j;
   double f1, f2;
 
   for (i = 0; i < 24; i++) {
@@ -281,12 +280,12 @@ static void set_freq_vertices(double freq_vertices[3][24][4],
 /* iw_zero=1 information can be used to omit to compute particles */
 /* interaction strength that is often heaviest part in throughout */
 /* calculation. */
-static int set_g(double g[3],
+static long set_g(double g[3],
-                 const double f0,
+                  const double f0,
-                 TPLCONST double freq_vertices[3][24][4],
+                  TPLCONST double freq_vertices[3][24][4],
-                 const long max_i)
+                  const long max_i)
 {
-  int i, iw_zero;
+  long i, iw_zero;
 
   iw_zero = 1;
 
@@ -302,9 +301,9 @@ static int set_g(double g[3],
   return iw_zero;
 }
 
-static int in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4])
+static long in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4])
 {
-  int i, j;
+  long i, j;
   double fmin, fmax;
 
   fmin = freq_vertices[0][0];
@@ -330,13 +329,13 @@ static int in_tetrahedra(const double f0, TPLCONST double freq_vertices[24][4])
 
 static void get_triplet_tetrahedra_vertices(
   long vertices[2][24][4],
-  TPLCONST int tp_relative_grid_address[2][24][4][3],
+  TPLCONST long tp_relative_grid_address[2][24][4][3],
-  const int mesh[3],
+  const long mesh[3],
   const long triplet[3],
-  TPLCONST int (*bz_grid_address)[3],
+  TPLCONST long (*bz_grid_address)[3],
   const long *bz_map)
 {
-  int i, j;
+  long i, j;
 
   for (i = 0; i < 2; i++) {
     for (j = 0; j < 24; j++) {

c/triplet_iw.h

@@ -35,7 +35,6 @@
 #ifndef __triplet_iw_H__
 #define __triplet_iw_H__
 
-#include <stddef.h>
 #include "triplet.h"
 
 void
@@ -43,18 +42,18 @@ tpi_get_integration_weight(double *iw,
                            char *iw_zero,
                            const double *frequency_points,
                            const long num_band0,
-                           TPLCONST int tp_relative_grid_address[2][24][4][3],
+                           TPLCONST long tp_relative_grid_address[2][24][4][3],
-                           const int mesh[3],
+                           const long mesh[3],
                            const long triplets[3],
                            const long num_triplets,
-                           TPLCONST int (*bz_grid_address)[3],
+                           TPLCONST long (*bz_grid_address)[3],
                            const long *bz_map,
                            const double *frequencies1,
                            const long num_band1,
                            const double *frequencies2,
                            const long num_band2,
                            const long tp_type,
-                           const int openmp_per_bands);
+                           const long openmp_per_bands);
 void tpi_get_integration_weight_with_sigma(double *iw,
                                            char *iw_zero,
                                            const double sigma,
@@ -66,14 +65,14 @@ void tpi_get_integration_weight_with_sigma(double *iw,
                                            const double *frequencies,
                                            const long num_band,
                                            const long tp_type,
-                                           const int openmp_per_bands);
+                                           const long openmp_per_bands);
 void
 tpi_get_dense_neighboring_grid_points(long neighboring_grid_points[],
                                       const long grid_point,
-                                      TPLCONST int relative_grid_address[][3],
+                                      TPLCONST long relative_grid_address[][3],
-                                      const int num_relative_grid_address,
+                                      const long num_relative_grid_address,
-                                      const int mesh[3],
+                                      const long mesh[3],
-                                      TPLCONST int bz_grid_address[][3],
+                                      TPLCONST long bz_grid_address[][3],
                                       const long bz_map[]);
 
 #endif

c/triplet_kpoint.c

@@ -36,14 +36,13 @@
 
 #include <stddef.h>
 #include <stdlib.h>
-#include "mathfunc.h"
 #include "kpoint.h"
-#include "kgrid.h"
+#include "rgrid.h"
 #include "triplet.h"
 #include "triplet_kpoint.h"
 
 #define KPT_NUM_BZ_SEARCH_SPACE 125
-static int bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
+static long bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
   { 0,  0,  0},
   { 0,  0,  1},
   { 0,  0,  2},
@@ -171,42 +170,42 @@ static int bz_search_space[KPT_NUM_BZ_SEARCH_SPACE][3] = {
   {-1, -1, -1}
 };
 
-static void grid_point_to_address_double(int address_double[3],
+static void grid_point_to_address_double(long address_double[3],
                                          const long grid_point,
-                                         const int mesh[3],
+                                         const long mesh[3],
-                                         const int is_shift[3]);
+                                         const long is_shift[3]);
 static long get_ir_triplets_at_q(long *map_triplets,
                                  long *map_q,
-                                 int (*grid_address)[3],
+                                 long (*grid_address)[3],
                                  const long grid_point,
-                                 const int mesh[3],
+                                 const long mesh[3],
-                                 const MatINT * rot_reciprocal,
+                                 const MatLONG * rot_reciprocal,
-                                 const int swappable);
+                                 const long swappable);
 static long get_BZ_triplets_at_q(long (*triplets)[3],
                                  const long grid_point,
-                                 TPLCONST int (*bz_grid_address)[3],
+                                 TPLCONST long (*bz_grid_address)[3],
                                  const long *bz_map,
                                  const long *map_triplets,
                                  const long num_map_triplets,
-                                 const int mesh[3]);
+                                 const long mesh[3]);
-static int get_third_q_of_triplets_at_q(int bz_address[3][3],
+static long get_third_q_of_triplets_at_q(long bz_address[3][3],
-                                        const int q_index,
+                                         const long q_index,
-                                        const long *bz_map,
+                                         const long *bz_map,
-                                        const int mesh[3],
+                                         const long mesh[3],
-                                        const int bzmesh[3]);
+                                         const long bzmesh[3]);
-static void modulo_i3(int v[3], const int m[3]);
+static void modulo_l3(long v[3], const long m[3]);
 
 long tpk_get_ir_triplets_at_q(long *map_triplets,
                               long *map_q,
-                              int (*grid_address)[3],
+                              long (*grid_address)[3],
-                              const int grid_point,
+                              const long grid_point,
-                              const int mesh[3],
+                              const long mesh[3],
-                              const int is_time_reversal,
+                              const long is_time_reversal,
-                              const MatINT * rotations,
+                              const MatLONG * rotations,
-                              const int swappable)
+                              const long swappable)
 {
-  int num_ir;
+  long num_ir;
-  MatINT *rot_reciprocal;
+  MatLONG *rot_reciprocal;
 
   rot_reciprocal = kpt_get_point_group_reciprocal(rotations, is_time_reversal);
   num_ir = get_ir_triplets_at_q(map_triplets,
@@ -216,17 +215,17 @@ long tpk_get_ir_triplets_at_q(long *map_triplets,
                                 mesh,
                                 rot_reciprocal,
                                 swappable);
-  mat_free_MatINT(rot_reciprocal);
+  kpt_free_MatLONG(rot_reciprocal);
   return num_ir;
 }
 
 long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
                               const long grid_point,
-                              TPLCONST int (*bz_grid_address)[3],
+                              TPLCONST long (*bz_grid_address)[3],
                               const long *bz_map,
                               const long *map_triplets,
                               const long num_map_triplets,
-                              const int mesh[3])
+                              const long mesh[3])
 {
   return get_BZ_triplets_at_q(triplets,
                               grid_point,
@@ -239,19 +238,19 @@ long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
 
 static long get_ir_triplets_at_q(long *map_triplets,
                                  long *map_q,
-                                 int (*grid_address)[3],
+                                 long (*grid_address)[3],
                                  const long grid_point,
-                                 const int mesh[3],
+                                 const long mesh[3],
-                                 const MatINT * rot_reciprocal,
+                                 const MatLONG * rot_reciprocal,
-                                 const int swappable)
+                                 const long swappable)
 {
   long i, j, num_grid, q_2, num_ir_q, num_ir_triplets, ir_grid_point;
-  int mesh_double[3], is_shift[3];
+  long mesh_double[3], is_shift[3];
-  int address_double0[3], address_double1[3], address_double2[3];
+  long address_double0[3], address_double1[3], address_double2[3];
   long *ir_grid_points, *third_q;
   double tolerance;
   double stabilizer_q[1][3];
-  MatINT *rot_reciprocal_q;
+  MatLONG *rot_reciprocal_q;
 
   ir_grid_points = NULL;
   third_q = NULL;
@@ -283,7 +282,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
                                                        mesh,
                                                        is_shift,
                                                        rot_reciprocal_q);
-  mat_free_MatINT(rot_reciprocal_q);
+  kpt_free_MatLONG(rot_reciprocal_q);
   rot_reciprocal_q = NULL;
 
   third_q = (long*) malloc(sizeof(long) * num_ir_q);
@@ -309,7 +308,7 @@ static long get_ir_triplets_at_q(long *map_triplets,
     for (j = 0; j < 3; j++) { /* q'' */
       address_double2[j] = - address_double0[j] - address_double1[j];
     }
-    third_q[i] = kgd_get_dense_grid_point_double_mesh(address_double2, mesh);
+    third_q[i] = rgd_get_double_grid_index(address_double2, mesh);
   }
 
   num_ir_triplets = 0;
@@ -348,15 +347,15 @@ static long get_ir_triplets_at_q(long *map_triplets,
 
 static long get_BZ_triplets_at_q(long (*triplets)[3],
                                  const long grid_point,
-                                 TPLCONST int (*bz_grid_address)[3],
+                                 TPLCONST long (*bz_grid_address)[3],
                                  const long *bz_map,
                                  const long *map_triplets,
                                  const long num_map_triplets,
-                                 const int mesh[3])
+                                 const long mesh[3])
 {
   long i, num_ir;
-  int j, k;
+  long j, k;
-  int bz_address[3][3], bz_address_double[3], bzmesh[3];
+  long bz_address[3][3], bz_address_double[3], bzmesh[3];
   long *ir_grid_points;
 
   ir_grid_points = NULL;
@@ -395,7 +394,7 @@ static long get_BZ_triplets_at_q(long (*triplets)[3],
         bz_address_double[k] = bz_address[j][k] * 2;
       }
       triplets[i][j] =
-        bz_map[kgd_get_dense_grid_point_double_mesh(bz_address_double, bzmesh)];
+        bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
     }
   }
 
@@ -405,20 +404,20 @@ static long get_BZ_triplets_at_q(long (*triplets)[3],
   return num_ir;
 }
 
-static int get_third_q_of_triplets_at_q(int bz_address[3][3],
+static long get_third_q_of_triplets_at_q(long bz_address[3][3],
-                                        const int q_index,
+                                         const long q_index,
-                                        const long *bz_map,
+                                         const long *bz_map,
-                                        const int mesh[3],
+                                         const long mesh[3],
-                                        const int bzmesh[3])
+                                         const long bzmesh[3])
 {
-  int i, j, smallest_g, smallest_index, sum_g, delta_g[3];
+  long i, j, smallest_g, smallest_index, sum_g, delta_g[3];
   long prod_bzmesh;
   long bzgp[KPT_NUM_BZ_SEARCH_SPACE];
-  int bz_address_double[3];
+  long bz_address_double[3];
 
   prod_bzmesh = (long)bzmesh[0] * bzmesh[1] * bzmesh[2];
 
-  modulo_i3(bz_address[q_index], mesh);
+  modulo_l3(bz_address[q_index], mesh);
   for (i = 0; i < 3; i++) {
     delta_g[i] = 0;
     for (j = 0; j < 3; j++) {
@@ -432,8 +431,7 @@ static int get_third_q_of_triplets_at_q(int bz_address[3][3],
       bz_address_double[j] = (bz_address[q_index][j] +
                               bz_search_space[i][j] * mesh[j]) * 2;
     }
-    bzgp[i] = bz_map[kgd_get_dense_grid_point_double_mesh(bz_address_double,
+    bzgp[i] = bz_map[rgd_get_double_grid_index(bz_address_double, bzmesh)];
-                                                          bzmesh)];
   }
 
   for (i = 0; i < KPT_NUM_BZ_SEARCH_SPACE; i++) {
@@ -449,9 +447,9 @@ escape:
 
   for (i = 0; i < KPT_NUM_BZ_SEARCH_SPACE; i++) {
     if (bzgp[i] < prod_bzmesh) { /* q'' is in BZ */
-      sum_g = (abs(delta_g[0] + bz_search_space[i][0]) +
+      sum_g = (labs(delta_g[0] + bz_search_space[i][0]) +
-               abs(delta_g[1] + bz_search_space[i][1]) +
+               labs(delta_g[1] + bz_search_space[i][1]) +
-               abs(delta_g[2] + bz_search_space[i][2]));
+               labs(delta_g[2] + bz_search_space[i][2]));
       if (sum_g < smallest_g) {
         smallest_index = i;
         smallest_g = sum_g;
@@ -466,13 +464,13 @@ escape:
   return smallest_g;
 }
 
-static void grid_point_to_address_double(int address_double[3],
+static void grid_point_to_address_double(long address_double[3],
                                          const long grid_point,
-                                         const int mesh[3],
+                                         const long mesh[3],
-                                         const int is_shift[3])
+                                         const long is_shift[3])
 {
-  int i;
+  long i;
-  int address[3];
+  long address[3];
 
 #ifndef GRID_ORDER_XYZ
   address[2] = grid_point / (mesh[0] * mesh[1]);
@@ -489,9 +487,9 @@ static void grid_point_to_address_double(int address_double[3],
   }
 }
 
-static void modulo_i3(int v[3], const int m[3])
+static void modulo_l3(long v[3], const long m[3])
 {
-  int i;
+  long i;
 
   for (i = 0; i < 3; i++) {
     v[i] = v[i] % m[i];

c/triplet_kpoint.h

@@ -37,24 +37,23 @@
 #ifndef __triplet_kpoint_H__
 #define __triplet_kpoint_H__
 
-#include <stddef.h>
+#include "kpoint.h"
-#include "mathfunc.h"
 #include "triplet.h"
 
 long tpk_get_ir_triplets_at_q(long *map_triplets,
                               long *map_q,
-                              int (*grid_address)[3],
+                              long (*grid_address)[3],
-                              const int grid_point,
+                              const long grid_point,
-                              const int mesh[3],
+                              const long mesh[3],
-                              const int is_time_reversal,
+                              const long is_time_reversal,
-                              const MatINT * rotations,
+                              const MatLONG * rotations,
-                              const int swappable);
+                              const long swappable);
 long tpk_get_BZ_triplets_at_q(long (*triplets)[3],
                               const long grid_point,
-                              TPLCONST int (*bz_grid_address)[3],
+                              TPLCONST long (*bz_grid_address)[3],
                               const long *bz_map,
                               const long *map_triplets,
                               const long num_map_triplets,
-                              const int mesh[3]);
+                              const long mesh[3]);
 
 #endif

phono3py/other/isotope.py

@@ -68,7 +68,7 @@ class Isotope(object):
                  symprec=1e-5,
                  cutoff_frequency=None,
                  lapack_zheev_uplo='L'):
-        self._mesh = np.array(mesh, dtype='intc')
+        self._mesh = np.array(mesh, dtype='int_')
 
         if mass_variances is None:
             self._mass_variances = get_mass_variances(primitive)
@@ -97,9 +97,9 @@ class Isotope(object):
 
         num_band = len(self._primitive) * 3
         if band_indices is None:
-            self._band_indices = np.arange(num_band, dtype='intc')
+            self._band_indices = np.arange(num_band, dtype='int_')
         else:
-            self._band_indices = np.array(band_indices, dtype='intc')
+            self._band_indices = np.array(band_indices, dtype='int_')
 
     def set_grid_point(self, grid_point):
         self._grid_point = grid_point
@@ -263,7 +263,7 @@ class Isotope(object):
         phono3c.integration_weights(
             self._integration_weights,
             freq_points,
-            np.array(thm.get_tetrahedra(), dtype='intc', order='C'),
+            thm.get_tetrahedra(),
             self._mesh,
             self._grid_points,
             self._frequencies,

phono3py/phonon/solver.py

@@ -83,16 +83,16 @@ def run_phonon_solver_c(dm,
         frequencies,
         eigenvectors,
         phonon_done,
-        np.array(grid_points, dtype='int_'),
+        grid_points,
-        np.array(grid_address, dtype='int_', order='C'),
+        grid_address,
         np.array(mesh, dtype='int_'),
         fc,
         svecs,
-        np.array(multiplicity, dtype='int_', order='C'),
+        multiplicity,
         positions,
         masses,
-        np.array(fc_p2s, dtype='int_'),
+        fc_p2s,
-        np.array(fc_s2p, dtype='int_'),
+        fc_s2p,
         frequency_conversion_factor,
         born,
         dielectric,
@@ -143,7 +143,7 @@ def _extract_params(dm):
         dielectric = None
 
     return (svecs,
-            multiplicity,
+            np.array(multiplicity, dtype='int_', order='C'),
             masses,
             rec_lattice,
             positions,
@@ -166,4 +166,4 @@ def _get_fc_elements_mapping(dm, fc):
         fc_p2s = np.arange(len(p2s_map), dtype='intc')
         fc_s2p = s2pp_map
 
-    return fc_p2s, fc_s2p
+    return np.array(fc_p2s, dtype='int_'), np.array(fc_s2p, dtype='int_')

phono3py/phonon3/collision_matrix.py

@@ -86,8 +86,8 @@ class CollisionMatrix(ImagSelfEnergy):
 
         if not self._is_reducible_collision_matrix:
             self._ir_grid_points = ir_grid_points
-            self._rot_grid_points = rot_grid_points
+            self._rot_grid_points = rot_grid_points  # dtype='int_'
-            self._point_operations = point_operations
+            self._point_operations = point_operations  # dtype='int_', order='C'
             self._primitive = self._pp.primitive
             rec_lat = np.linalg.inv(self._primitive.cell)
             self._rotations_cartesian = np.array(

phono3py/phonon3/conductivity.py

@@ -69,8 +69,8 @@ def write_pp(conductivity,
     sigma_cutoff = conductivity.get_sigma_cutoff_width()
     mesh = conductivity.get_mesh_numbers()
     triplets, weights, map_triplets, _ = pp.get_triplets_at_q()
-    grid_address = pp.get_grid_address()
+    grid_address = pp.grid_address
-    bz_map = pp.get_bz_map()
+    bz_map = pp.bz_map
     if map_triplets is None:
         all_triplets = None
     else:
@@ -138,11 +138,12 @@ class Conductivity(object):
         self._symmetry = symmetry
 
         if not self._is_kappa_star:
-            self._point_operations = np.array([np.eye(3, dtype='intc')],
+            self._point_operations = np.array([np.eye(3, dtype='int_')],
-                                              dtype='intc')
+                                              dtype='int_', order='C')
         else:
-            self._point_operations = symmetry.get_reciprocal_operations()
+            self._point_operations = np.array(symmetry.reciprocal_operations,
-        rec_lat = np.linalg.inv(self._primitive.get_cell())
+                                              dtype='int_', order='C')
+        rec_lat = np.linalg.inv(self._primitive.cell)
         self._rotations_cartesian = np.array(
             [similarity_transformation(rec_lat, r)
              for r in self._point_operations], dtype='double')
@@ -338,7 +339,8 @@ class Conductivity(object):
             (self._ir_grid_points,
              self._ir_grid_weights) = self._get_ir_grid_points()
         elif not self._is_kappa_star:  # All grid points
-            coarse_grid_address = get_grid_address(self._coarse_mesh)
+            coarse_grid_address = np.array(get_grid_address(self._coarse_mesh),
+                                           dtype='int_', order='C')
             coarse_grid_points = np.arange(np.prod(self._coarse_mesh),
                                            dtype='int_')
             self._grid_points = from_coarse_to_dense_grid_points(
@@ -347,7 +349,7 @@ class Conductivity(object):
                 coarse_grid_points,
                 coarse_grid_address,
                 coarse_mesh_shifts=self._coarse_mesh_shifts)
-            self._grid_weights = np.ones(len(self._grid_points), dtype='intc')
+            self._grid_weights = np.ones(len(self._grid_points), dtype='int_')
             self._ir_grid_points = self._grid_points
             self._ir_grid_weights = self._grid_weights
         else:  # Automatic sampling
@@ -394,7 +396,7 @@ class Conductivity(object):
         self._mesh = self._pp.mesh_numbers
 
         if mesh_divisors is None:
-            self._mesh_divisors = np.array([1, 1, 1], dtype='intc')
+            self._mesh_divisors = np.array([1, 1, 1], dtype='int_')
         else:
             self._mesh_divisors = []
             for i, (m, n) in enumerate(zip(self._mesh, mesh_divisors)):
@@ -405,7 +407,7 @@ class Conductivity(object):
                     print(("Mesh number %d for the " +
                            ["first", "second", "third"][i] +
                            " axis is not dividable by divisor %d.") % (m, n))
-            self._mesh_divisors = np.array(self._mesh_divisors, dtype='intc')
+            self._mesh_divisors = np.array(self._mesh_divisors, dtype='int_')
             if coarse_mesh_shifts is None:
                 self._coarse_mesh_shifts = [False, False, False]
             else:

phono3py/phonon3/conductivity_LBTE.py

@@ -1115,7 +1115,7 @@ class Conductivity_LBTE(Conductivity):
                 self._average_collision_matrix_by_degeneracy()
                 self._expand_collisions()
             self._combine_reducible_collisions()
-            weights = np.ones(np.prod(self._mesh), dtype='intc')
+            weights = np.ones(np.prod(self._mesh), dtype='int_')
             self._symmetrize_collision_matrix()
         else:
             self._combine_collisions()

phono3py/phonon3/conductivity_RTA.py

@@ -158,8 +158,8 @@ def _write_gamma_detail(br, interaction, i, compression="gzip", filename=None,
     sigmas = br.get_sigmas()
     sigma_cutoff = br.get_sigma_cutoff_width()
     triplets, weights, map_triplets, _ = interaction.get_triplets_at_q()
-    grid_address = interaction.get_grid_address()
+    grid_address = interaction.grid_address
-    bz_map = interaction.get_bz_map()
+    bz_map = interaction.bz_map
     if map_triplets is None:
         all_triplets = None
     else:
@@ -832,7 +832,7 @@ class Conductivity_RTA(Conductivity):
             if sigma is None:
                 phono3c.pp_collision(
                     collisions,
-                    np.array(thm.get_tetrahedra(), dtype='intc', order='C'),
+                    thm.get_tetrahedra(),
                     self._frequencies,
                     self._eigenvectors,
                     triplets_at_q,

phono3py/phonon3/fc3.py

@@ -156,7 +156,7 @@ def distribute_fc3(fc3,
             rot_indices = np.where(permutations[:, i_target] == i_done)[0]
             if len(rot_indices) > 0:
                 atom_mapping = np.array(permutations[rot_indices[0]],
-                                        dtype='intc')
+                                        dtype='int_')
                 rot = rotations[rot_indices[0]]
                 rot_cart_inv = np.array(
                     similarity_transformation(lattice, rot).T,
@@ -175,8 +175,8 @@ def distribute_fc3(fc3,
         try:
             import phono3py._phono3py as phono3c
             phono3c.distribute_fc3(fc3,
-                                   int(s2compact[i_target]),
+                                   s2compact[i_target],
-                                   int(s2compact[i_done]),
+                                   s2compact[i_done],
                                    atom_mapping,
                                    rot_cart_inv)
         except ImportError:
@@ -206,18 +206,19 @@ def set_permutation_symmetry_fc3(fc3):
 def set_permutation_symmetry_compact_fc3(fc3, primitive):
     try:
         import phono3py._phono3py as phono3c
-        s2p_map = primitive.get_supercell_to_primitive_map()
+        s2p_map = primitive.s2p_map
-        p2s_map = primitive.get_primitive_to_supercell_map()
+        p2s_map = primitive.p2s_map
-        p2p_map = primitive.get_primitive_to_primitive_map()
+        p2p_map = primitive.p2p_map
         permutations = primitive.get_atomic_permutations()
         s2pp_map, nsym_list = get_nsym_list_and_s2pp(s2p_map,
                                                      p2p_map,
                                                      permutations)
-        phono3c.permutation_symmetry_compact_fc3(fc3,
+        phono3c.permutation_symmetry_compact_fc3(
-                                                 permutations,
+            fc3,
-                                                 s2pp_map,
+            np.array(permutations, dtype='int_', order='C'),
-                                                 p2s_map,
+            np.array(s2pp_map, dtype='int_'),
-                                                 nsym_list)
+            np.array(p2s_map, dtype='int_'),
+            np.array(nsym_list, dtype='int_'))
     except ImportError:
         text = ("Import error at phono3c.permutation_symmetry_compact_fc3. "
                 "Corresponding python code is not implemented.")
@@ -254,26 +255,33 @@ def set_translational_invariance_fc3(fc3):
 def set_translational_invariance_compact_fc3(fc3, primitive):
     try:
         import phono3py._phono3py as phono3c
-        s2p_map = primitive.get_supercell_to_primitive_map()
+        s2p_map = primitive.s2p_map
-        p2s_map = primitive.get_primitive_to_supercell_map()
+        p2s_map = primitive.p2s_map
-        p2p_map = primitive.get_primitive_to_primitive_map()
+        p2p_map = primitive.p2p_map
         permutations = primitive.get_atomic_permutations()
         s2pp_map, nsym_list = get_nsym_list_and_s2pp(s2p_map,
                                                      p2p_map,
                                                      permutations)
-        phono3c.transpose_compact_fc3(fc3,
+
-                                      permutations,
+        permutations = np.array(permutations, dtype='int_', order='C')
-                                      s2pp_map,
+        s2pp_map = np.array(s2pp_map, dtype='int_')
-                                      p2s_map,
+        p2s_map = np.array(p2s_map, dtype='int_')
-                                      nsym_list,
+        nsym_list = np.array(nsym_list, dtype='int_')
-                                      0)  # dim[0] <--> dim[1]
+        phono3c.transpose_compact_fc3(
+            fc3,
+            permutations,
+            s2pp_map,
+            p2s_map,
+            nsym_list,
+            0)  # dim[0] <--> dim[1]
         set_translational_invariance_fc3_per_index(fc3, index=1)
-        phono3c.transpose_compact_fc3(fc3,
+        phono3c.transpose_compact_fc3(
-                                      permutations,
+            fc3,
-                                      s2pp_map,
+            permutations,
-                                      p2s_map,
+            s2pp_map,
-                                      nsym_list,
+            p2s_map,
-                                      0)  # dim[0] <--> dim[1]
+            nsym_list,
+            0)  # dim[0] <--> dim[1]
         set_translational_invariance_fc3_per_index(fc3, index=1)
         set_translational_invariance_fc3_per_index(fc3, index=2)
 
@@ -441,6 +449,7 @@ def solve_fc3(first_atom_num,
                                                  positions,
                                                  site_symmetry,
                                                  symprec)
+    rot_map_syms = np.array(rot_map_syms, dtype='int_', order='C')
     rot_disps = get_rotated_displacement(displacements_first, site_sym_cart)
 
     logger.debug("pinv")

phono3py/phonon3/imag_self_energy.py

@@ -588,7 +588,7 @@ class ImagSelfEnergy(object):
             self._pp_strength,
             self._triplets_at_q,
             self._weights_at_q,
-            self._pp.get_grid_address(),
+            self._pp.grid_address,
             self._frequencies,
             self._temperature,
             self._g,

phono3py/phonon3/interaction.py

@@ -62,7 +62,7 @@ class Interaction(object):
                  lapack_zheev_uplo='L'):
         self._supercell = supercell
         self._primitive = primitive
-        self._mesh = np.array(mesh, dtype='intc')
+        self._mesh = np.array(mesh, dtype='int_')
         self._symmetry = symmetry
 
         self._band_indices = None
@@ -114,10 +114,10 @@ class Interaction(object):
 
         svecs, multiplicity = self._primitive.get_smallest_vectors()
         self._smallest_vectors = svecs
-        self._multiplicity = multiplicity
+        self._multiplicity = np.array(multiplicity, dtype='int_', order='C')
         self._masses = np.array(self._primitive.masses, dtype='double')
-        self._p2s = self._primitive.p2s_map
+        self._p2s = np.array(self._primitive.p2s_map, dtype='int_')
-        self._s2p = self._primitive.s2p_map
+        self._s2p = np.array(self._primitive.s2p_map, dtype='int_')
 
         self._allocate_phonon()
 
@@ -468,9 +468,9 @@ class Interaction(object):
     def _set_band_indices(self, band_indices):
         num_band = len(self._primitive) * 3
         if band_indices is None:
-            self._band_indices = np.arange(num_band, dtype='intc')
+            self._band_indices = np.arange(num_band, dtype='int_')
         else:
-            self._band_indices = np.array(band_indices, dtype='intc')
+            self._band_indices = np.array(band_indices, dtype='int_')
 
     def _run_c(self, g_zero):
         import phono3py._phono3py as phono3c

phono3py/phonon3/triplets.py

@@ -53,7 +53,7 @@ def get_triplets_at_q(grid_point,
         A grid point
     mesh : array_like
         Mesh numbers
-        shape=(3,), dtype='intc'
+        shape=(3,), dtype='int_'
     point_group : array_like
         Rotation matrices in real space. Note that those in reciprocal space
         mean these matrices transposed (local terminology).
@@ -109,11 +109,12 @@ def get_triplets_at_q(grid_point,
         is_time_reversal=is_time_reversal,
         swappable=swappable)
     bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
-        grid_address,
+        np.array(grid_address, dtype='intc', order='C'),
         mesh,
         reciprocal_lattice,
         is_dense=True)
     bz_map = np.array(bz_map, dtype='int_')
+    bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
     triplets_at_q, weights = _get_BZ_triplets_at_q(
         grid_point,
         bz_grid_address,
@@ -151,14 +152,14 @@ def get_nosym_triplets_at_q(grid_point,
                             mesh,
                             reciprocal_lattice,
                             stores_triplets_map=False):
-    grid_address = get_grid_address(mesh)
     bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
-        grid_address,
+        get_grid_address(mesh),
         mesh,
         reciprocal_lattice,
         is_dense=True)
     bz_map = np.array(bz_map, dtype='int_')
-    map_triplets = np.arange(len(grid_address), dtype=bz_map.dtype)
+    bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
+    map_triplets = np.arange(np.prod(mesh), dtype=bz_map.dtype)
     triplets_at_q, weights = _get_BZ_triplets_at_q(
         grid_point,
         bz_grid_address,
@@ -176,6 +177,7 @@ def get_nosym_triplets_at_q(grid_point,
 
 
 def get_grid_address(mesh):
+    """Returns grid_address of dtype='intc'"""
     grid_mapping_table, grid_address = spglib.get_stabilized_reciprocal_mesh(
         mesh,
         [[[1, 0, 0], [0, 1, 0], [0, 0, 1]]],
@@ -192,6 +194,8 @@ def get_bz_grid_address(mesh, reciprocal_lattice, with_boundary=False):
         mesh,
         reciprocal_lattice,
         is_dense=True)
+    bz_map = np.array(bz_map, dtype='int_')
+    bz_grid_address = np.array(bz_grid_address, dtype='int_', order='C')
     if with_boundary:
         return bz_grid_address, bz_map
     else:
@@ -237,6 +241,8 @@ def get_ir_grid_points(mesh, rotations, mesh_shifts=None):
         rotations,
         is_shift=np.where(mesh_shifts, 1, 0),
         is_dense=True)
+    grid_mapping_table = np.array(grid_mapping_table, dtype='int_', order='C')
+    grid_address = np.array(grid_address, dtype='int_', order='C')
     (ir_grid_points,
      ir_grid_weights) = extract_ir_grid_points(grid_mapping_table)
 
@@ -262,11 +268,11 @@ def reduce_grid_points(mesh_divisors,
                        dense_grid_points,
                        dense_grid_weights=None,
                        coarse_mesh_shifts=None):
-    divisors = np.array(mesh_divisors, dtype='intc')
+    divisors = np.array(mesh_divisors, dtype='int_')
     if (divisors == 1).all():
         coarse_grid_points = np.array(dense_grid_points, dtype='int_')
         if dense_grid_weights is not None:
-            coarse_grid_weights = np.array(dense_grid_weights, dtype='intc')
+            coarse_grid_weights = np.array(dense_grid_weights, dtype='int_')
     else:
         if coarse_mesh_shifts is None:
             shift = [0, 0, 0]
@@ -306,10 +312,10 @@ def get_coarse_ir_grid_points(primitive,
                               coarse_mesh_shifts,
                               is_kappa_star=True,
                               symprec=1e-5):
-    mesh = np.array(mesh, dtype='intc')
+    mesh = np.array(mesh, dtype='int_')
 
     symmetry = Symmetry(primitive, symprec)
-    point_group = symmetry.get_pointgroup_operations()
+    point_group = symmetry.pointgroup_operations
 
     if mesh_divisors is None:
         (ir_grid_points,
@@ -317,13 +323,14 @@ def get_coarse_ir_grid_points(primitive,
          grid_address,
          grid_mapping_table) = get_ir_grid_points(mesh, point_group)
     else:
-        mesh_divs = np.array(mesh_divisors, dtype='intc')
+        mesh_divs = np.array(mesh_divisors, dtype='int_')
         coarse_mesh = mesh // mesh_divs
         if coarse_mesh_shifts is None:
             coarse_mesh_shifts = [False, False, False]
 
         if not is_kappa_star:
-            coarse_grid_address = get_grid_address(coarse_mesh)
+            coarse_grid_address = np.array(get_grid_address(coarse_mesh),
+                                           dtype='int_', order='C')
             coarse_grid_points = np.arange(np.prod(coarse_mesh), dtype='int_')
         else:
             (coarse_ir_grid_points,
@@ -342,9 +349,9 @@ def get_coarse_ir_grid_points(primitive,
         grid_address = get_grid_address(mesh)
         ir_grid_weights = ir_grid_weights
 
-    reciprocal_lattice = np.linalg.inv(primitive.get_cell())
+    reciprocal_lattice = np.linalg.inv(primitive.cell)
     bz_grid_address, bz_map = spglib.relocate_BZ_grid_address(
-        grid_address,
+        np.array(grid_address, dtype='intc', order='C'),
         mesh,
         reciprocal_lattice,
         is_dense=True)
@@ -360,10 +367,9 @@ def get_number_of_triplets(primitive,
                            grid_point,
                            swappable=True,
                            symprec=1e-5):
-    mesh = np.array(mesh, dtype='intc')
     symmetry = Symmetry(primitive, symprec)
-    point_group = symmetry.get_pointgroup_operations()
+    point_group = symmetry.pointgroup_operations
-    reciprocal_lattice = np.linalg.inv(primitive.get_cell())
+    reciprocal_lattice = np.linalg.inv(primitive.cell)
     triplets_at_q, _, _, _, _, _ = get_triplets_at_q(
         grid_point,
         mesh,
@@ -497,16 +503,16 @@ def _get_triplets_reciprocal_mesh_at_q(fixed_grid_number,
 
     map_triplets = np.zeros(np.prod(mesh), dtype='int_')
     map_q = np.zeros(np.prod(mesh), dtype='int_')
-    grid_address = np.zeros((np.prod(mesh), 3), dtype='intc')
+    grid_address = np.zeros((np.prod(mesh), 3), dtype='int_')
 
     phono3c.triplets_reciprocal_mesh_at_q(
         map_triplets,
         map_q,
         grid_address,
         fixed_grid_number,
-        np.array(mesh, dtype='intc'),
+        np.array(mesh, dtype='int_'),
         is_time_reversal * 1,
-        np.array(rotations, dtype='intc', order='C'),
+        np.array(rotations, dtype='int_', order='C'),
         swappable * 1)
 
     return map_triplets, map_q, grid_address
@@ -530,7 +536,7 @@ def _get_BZ_triplets_at_q(grid_point,
                                           bz_grid_address,
                                           bz_map,
                                           map_triplets,
-                                          np.array(mesh, dtype='intc'))
+                                          np.array(mesh, dtype='int_'))
     assert num_ir_ret == len(ir_weights)
 
     return triplets, np.array(ir_weights, dtype='int_')
@@ -558,18 +564,19 @@ def _set_triplets_integration_weights_c(g,
             phono3c.neighboring_grid_points(
                 neighboring_grid_points,
                 np.array(triplets_at_q[:, i], dtype='int_').ravel(),
-                np.array(j * unique_vertices, dtype='intc', order='C'),
+                np.array(j * unique_vertices, dtype='int_', order='C'),
                 mesh,
                 grid_address,
                 bz_map)
-            interaction.run_phonon_solver(np.unique(neighboring_grid_points))
+            interaction.run_phonon_solver(
+                np.array(np.unique(neighboring_grid_points), dtype='int_'))
 
     frequencies = interaction.get_phonons()[0]
     phono3c.triplets_integration_weights(
         g,
         g_zero,
         frequency_points,  # f0
-        np.array(thm.get_tetrahedra(), dtype='intc', order='C'),
+        thm.get_tetrahedra(),
         mesh,
         triplets_at_q,
         frequencies,  # f1
@@ -581,10 +588,10 @@ def _set_triplets_integration_weights_c(g,
 
 def _set_triplets_integration_weights_py(g, interaction, frequency_points):
     reciprocal_lattice = np.linalg.inv(interaction.get_primitive().get_cell())
-    mesh = interaction.get_mesh_numbers()
+    mesh = interaction.mesh_numbers
     thm = TetrahedronMethod(reciprocal_lattice, mesh=mesh)
-    grid_address = interaction.get_grid_address()
+    grid_address = interaction.grid_address
-    bz_map = interaction.get_bz_map()
+    bz_map = interaction.bz_map
     triplets_at_q = interaction.get_triplets_at_q()[0]
     tetrahedra_vertices = get_tetrahedra_vertices(
         thm.get_tetrahedra(),
@@ -592,7 +599,8 @@ def _set_triplets_integration_weights_py(g, interaction, frequency_points):
         triplets_at_q,
         grid_address,
         bz_map)
-    interaction.run_phonon_solver(np.unique(tetrahedra_vertices))
+    interaction.run_phonon_solver(
+        np.array(np.unique(tetrahedra_vertices), dtype='int_'))
     frequencies = interaction.get_phonons()[0]
     num_band = frequencies.shape[1]
     for i, vertices in enumerate(tetrahedra_vertices):