text: #include "spglib.h"
  #include "tetrahedron_method.h"
  #include <stdio.h>
  #include <stdlib.h>
  
  static void test_tetrahedron_method(void);
  static void mat_copy_matrix_d3(double a[3][3], double b[3][3]);
  static double mat_get_determinant_d3(double a[3][3]);
  static int mat_inverse_matrix_d3(double m[3][3],
  				 double a[3][3],
  				 const double precision);
  
  int main(void)
  {
    test_tetrahedron_method();
  
    return 0;
  }
  
  /* frequency.dat is in the example directory. */
  /* The values in this file are the phonon frequencies of NaCl */
  /* with 20x20x20 mesh. Calculation was done with reducing */
  /* k-points to the irreducible k-points using phonopy. */
  /* (http://phonopy.sf.net/) */
  static void test_tetrahedron_method(void)
  {
    printf("*** Example of tetrahedron method of NaCl to calculate DOS ***:\n");
    printf("Read data from frequency.dat and write DOS to dos.dat.\n");
  
    int i, j, k, l, q, r;
  
    /* NaCl 20x20x20 mesh */
    double lattice[3][3] = {
      {0.000000000000000, 2.845150738087836, 2.845150738087836},
      {2.845150738087836, 0.000000000000000, 2.845150738087836},
      {2.845150738087836, 2.845150738087836, 0.000000000000000}
    };
    double position[][3] =
      {{0, 0, 0},
       {0.5, 0.5, 0.5}};
    int types[] = {1, 2};
    int num_atom = 2;
    int m = 20;
    int mesh[] = {m, m, m};
    int num_gp = mesh[0] * mesh[1] * mesh[2];
    int is_shift[] = {0, 0, 0};
    int grid_address[num_gp][3];
    int grid_mapping_table[num_gp];
    int weights[num_gp];
    int num_ir = spg_get_ir_reciprocal_mesh(grid_address,
  					  grid_mapping_table,
  					  mesh,
  					  is_shift,
  					  1,
  					  lattice,
  					  position,
  					  types,
  					  num_atom,
  					  1e-5);
    int ir_gp[num_ir];
    int ir_weights[num_ir];
    int gp_ir_index[num_gp];
    int relative_grid_address[24][4][3];
    double rec_lat[3][3];
    FILE *fp;
    char * line = NULL;
    size_t len = 0;
    ssize_t read;
    double frequency[num_ir * num_atom * 3];
    double max_f, min_f;
    double t_omegas[24][4];
    int g_addr[3];
    int gp;
    int num_freqs = 201;
    double dos[num_freqs];
    double integral_dos[num_freqs];
    double omegas[num_freqs];
    double iw;
    
    for (i = 0; i < num_gp; i++) {
      weights[i] = 0;
    }
  
    for (i = 0; i < num_gp; i++) {
      weights[grid_mapping_table[i]]++;
    }
  
    j = 0;
    for (i = 0; i < num_gp; i++) {
      if (weights[i] != 0) {
        ir_gp[j] = i;
        ir_weights[j] = weights[i];
        gp_ir_index[i] = j;
        j++;
      } else {
        gp_ir_index[i] = gp_ir_index[grid_mapping_table[i]];
      }
    }
  
    printf("Number of irreducible k-points: %d\n", num_ir);
    
    mat_inverse_matrix_d3(rec_lat, lattice, 1e-5);
    thm_get_relative_grid_address(relative_grid_address, rec_lat);
  
    /* for (i = 0; i < 24; i++) { */
    /*   for (j = 0; j < 4; j++) { */
    /*     printf("[%2d %2d %2d] ", */
    /* 	     relative_grid_address[i][j][0], */
    /* 	     relative_grid_address[i][j][1], */
    /* 	     relative_grid_address[i][j][2]); */
    /*   } */
    /*   printf("\n"); */
    /* } */
  
    fp = fopen("frequency.dat", "r");
  
    for (i = 0; i < num_ir * num_atom * 3; i++) {
      read = getline(&line, &len, fp);
      if (read == -1) {
        break;
      }
      frequency[i] = strtod(line, NULL);
    }
  
    fclose(fp);
  
    max_f = frequency[0];
    min_f = frequency[0];
    for (i = 0; i < num_ir * num_atom * 3; i++) {
      if (max_f < frequency[i]) {
        max_f = frequency[i];
      }
      if (min_f > frequency[i]) {
        min_f = frequency[i];
      }
    }
    
    printf("Number of frequencies: %d\n", i);
    
  
  #pragma omp parallel for private(j, k, l, q, r, g_addr, gp, t_omegas, iw)
    for (i = 0; i < num_freqs; i++) {
      dos[i] = 0;
      integral_dos[i] = 0;
      omegas[i] = min_f + (max_f - min_f) / (num_freqs - 1) * i;
      for (j = 0; j < num_ir;  j++) {
        for (k = 0; k < num_atom * 3; k++) {
  	for (l = 0; l < 24; l++) {
  	  for (q = 0; q < 4; q++) {
  	    for (r = 0; r < 3; r++) {
  	      g_addr[r] = grid_address[ir_gp[j]][r] +
  		relative_grid_address[l][q][r];
  	    }
  	    gp = spg_get_grid_point_from_address(g_addr, mesh);
  	    t_omegas[l][q] = frequency[gp_ir_index[gp] * num_atom * 3 + k];
  	  }
  	}
  	iw = thm_get_integration_weight(omegas[i], t_omegas, 'J');
  	dos[i] += iw * ir_weights[j];
  	iw = thm_get_integration_weight(omegas[i], t_omegas, 'I');
  	integral_dos[i] += iw * ir_weights[j];
        }
      }
    }
  
    fp = fopen("dos.dat", "w");
  
    for (i = 0; i < num_freqs; i++) {
      fprintf(fp, "%f %f\n", omegas[i], dos[i] / num_gp);
    }
  
    fprintf(fp, "\n\n");
    
    for (i = 0; i < num_freqs; i++) {
      fprintf(fp, "%f %f\n", omegas[i], integral_dos[i] / num_gp);
    }
      
    fclose(fp);
  }
  
  static void mat_copy_matrix_d3(double a[3][3], double b[3][3])
  {
    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];
  }
  
  static double mat_get_determinant_d3(double a[3][3])
  {
    return a[0][0] * (a[1][1] * a[2][2] - a[1][2] * a[2][1])
      + a[0][1] * (a[1][2] * a[2][0] - a[1][0] * a[2][2])
      + a[0][2] * (a[1][0] * a[2][1] - a[1][1] * a[2][0]);
  }
  
  static int mat_inverse_matrix_d3(double m[3][3],
  				 double a[3][3],
  				 const double precision)
  {
    double det;
    double c[3][3];
    det = mat_get_determinant_d3(a);
  
    c[0][0] = (a[1][1] * a[2][2] - a[1][2] * a[2][1]) / det;
    c[1][0] = (a[1][2] * a[2][0] - a[1][0] * a[2][2]) / det;
    c[2][0] = (a[1][0] * a[2][1] - a[1][1] * a[2][0]) / det;
    c[0][1] = (a[2][1] * a[0][2] - a[2][2] * a[0][1]) / det;
    c[1][1] = (a[2][2] * a[0][0] - a[2][0] * a[0][2]) / det;
    c[2][1] = (a[2][0] * a[0][1] - a[2][1] * a[0][0]) / det;
    c[0][2] = (a[0][1] * a[1][2] - a[0][2] * a[1][1]) / det;
    c[1][2] = (a[0][2] * a[1][0] - a[0][0] * a[1][2]) / det;
    c[2][2] = (a[0][0] * a[1][1] - a[0][1] * a[1][0]) / det;
    mat_copy_matrix_d3(m, c);
    return 1;
  }
  
  